The description of the codes for removal of cosmic ray signatures on Deep Impact images is presented in section 2 of our paper on cosmic rays: http://www.astro.umd.edu/~ipatov/adsr-crs-2007.pdf
Below we present this section with the added links to the codes:
2.
Codes for removal of cosmic ray signatures
The most reliable way to recognize CR events on
astronomical images is to compare different images of the same region of the
sky, but it is not always possible to do this. The previous codes for removal of CR signatures in a single image were
worked out for studies of typical sky images. (Examples are imgclean http://www.astro.umd.edu/~ipatov/imclean.pro and crfind http://www.astro.umd.edu/~ipatov/crfind.pro, written by E. Deutsch and R. White, respectively. Imgclean can be found at several
websites, including http://pdssbn.astro.umd.edu/volume/didoc_0001/document/calibration_software/dical_v5/imgclean.pro.)
Sometimes the above codes do not work
well with visual images of a comet and with IR and raw (before calibration)
visual images. The code di_crrej http://www.astro.umd.edu/~ipatov/di_crrej.pro
written by D. Lindler in 2005 for analysis of DI images
has similar problems. All codes used can be found at http://www.astro.umd.edu/~ipatov/CRcodes.html.
These codes sometimes detect false CR signatures near the edge of a comet
nucleus or on its coma and may have problems with long (oblique entry) CR
signatures. Crfind
and di_crrej
only identify pixels corresponding to CR signatures, but do not replace these
pixels. We studied the above codes because they were readily available to us.
There are other codes (e.g., those described by van Dokkum
(2001) and Pych (2004)) that reject cosmic rays in
single CCD exposures; these have not yet been studied for use with DI images.
Imgclean uses the code starchck http://www.astro.umd.edu/~ipatov/starchck.pro
in order to make a decision whether a cluster corresponds to a star or to a
cosmic ray. The latter code compares the brightness of the brightest pixel of a
cluster with the brightness of several squares of close background and studies
the location of the pixel in the cluster. It is recommended to run imgclean twice on
the same image to achieve best cleaning results. In order to use imgclean for
analysis of DI images, we made a few changes. Our code imgcleanf http://www.astro.umd.edu/~ipatov/imcleanf.pro
(http://pdssbn.astro.umd.edu/volume/didoc_0001/document/calibration_software/dical_v5/imgcleanf.pro)
uses imgclean
as a subroutine and: (1) does not search for CRs at
the edges of Deep Impact MRI and HRI VIS full-frame images; (2) replaces the
brightness of clusters detected as CRs in a different
way.
Imgcleanf http://www.astro.umd.edu/~ipatov/imcleanf.pro
does not consider pixels located close to the edges of the frame as CR events
because these pixels do not contain an image itself. The sizes of these edges
are presented in Table IV in (Hampton et al., 2005), but usually we considered
the edges wider by 1-3 pixels (depending on the mode of the frame) than the
edges from the table. For example, for mode 1 (1024×1024 pixels) Hampton et al.
(2005) considered the edge of 8 pixels, but our analysis of images showed that
sometimes even the 11th line can
contain some bright pixels that are not CR signatures.
New replacement of the brightnesses
of pixels belonging to signatures of CRs was needed
because near the edge of a comet nucleus, imgclean considers too many
pixels as CRs and replaces their brightness with a
wrong brightness, so the edge of an image of the nucleus after imgclean can be
become ‘torn’ (though it must be smooth), part of the nucleus becomes 'eaten',
and a lot of bright pixels appear near the edge of the nucleus. With imgcleanf,
replaced pixels get the brightness of the median value of the adjacent
background. At the beginning of the IDL code we calculate sky=median(ima, medima), where ima is the
initial image and the default value of medima equals 7 (medima=7 was chosen based on analysis of corrected images
obtained at different values of medima). The
brightness of the pixels belonging to CR events is replaced by sky values. If this brightness is
greater than the initial brightness, then the initial brightness is left. This
approach allows pixels deleted as CRs (including
false CRs) near the edge of the nucleus to be reset
to practically the same level as the nearest background. So the deleted pixels
do not spoil the final image of a comet the way they do with imgclean.
Crfind http://www.astro.umd.edu/~ipatov/crfind.pro
removes signatures of cosmic rays on a single image using shape parameters. The
description of the code presented below is based on that presented at the
beginning of the code. Basically it rejects pixels that are too bright compared
with the median of neighboring pixel values. Iteratively it recomputes
the median after each rejection cycle so that large CRs
are gradually eliminated from the image. It calculates kcrf=[(IMAGE-SKY)/PSFRAT-(mfi-SKY)]/sigma, where sigma is the
noise expected in the difference, mfi is the median filtered image, PSFRAT is the ratio of central pixel
value to the median of surrounding pixels for a point-spread function centered
on pixel (default value equals 2.0, appropriate for DI PSFs),
SKY is an estimate of sky background
(default is to compute 15×15 pixels) median. If the value of kcrf
is greater than THRESH1×sigma, the pixel is flagged as bad.
Initially the C routine CMEDMASK was used for the initial median
computation, but then D. Lindler modified the code so
that it no longer needs a C routine at
the cost of a decrease in execution speed. It does a 3×3 median ignoring the
central pixel and any masked pixels. The median is recomputed with bad pixels
removed and then tested again with a lower threshold THRESH2, but only pixels that are adjacent to CRs
are allowed to be flagged. This algorithm is iterated as long as any new bad
pixels are found. The result is that CRs are allowed
to "grow" to include neighboring pixels that are likely to be bad.
Noise in a CCD pixel with D mean
counts (in digital numbers DN) is assumed to be sigma=((D+convar)/gain+(flin×D)2)1/2,
where convar=readnoise2×gain, readnoise
is CCD readout noise in DN (default value is 0.77), and gain is CCD gain in electrons/DN
(default value is 28.5). This is just the normal CCD noise model with the noise
becoming linear in D for large D. The default value of flin equals 0 and
means that there is no linear term. To customize for Deep Impact, D. Linder set
default THRESH1 to 3.5 and THRESH2 to 1.5.
The routine di_crrej http://www.astro.umd.edu/~ipatov/di_crrej.pro
uses the changes in image values from one pixel to the next along with a filter
to determine if the pixel is on a "sharp" edge of the image. It finds
both cosmic rays and small "holes" in the images.
We wrote the code rmcr for removal
of CR signatures. In contrast to the above codes, it can work also with raw
(not only with calibrated) images, but the present version is slow if an image
of a comet consists of a large number of pixels. The code rmcr analyses those pixels for
which raw digital numbers (DNs) or calibrated
radiances are greater than some limit lim. Only such
pixels are considered to be caused by CRs. This limit
can be an input parameter (e.g., for raw images) or it can be calculated as lim=limit0*klim (e.g., klim=10), where limit0 is the median value of all pixels on an image. For
calibrated images, one may not know lim in advance, so it is better to use the latter
calculation of lim.
In a ‘fixed’ version of rmcr http://www.astro.umd.edu/~ipatov/gjfrmcr.pro,
the value of lim
is the same for an entire image. In a ‘variable’ version http://www.astro.umd.edu/~ipatov/gjrmcr.pro,
limij=limitij*klim is different for different pixels with coordinates
i and j.
First we calculate limitij=median(ima,medima). If this
value is smaller than lim,
then limitij=lim (note
that if we do not use the second check, then sometimes we get many false CR
signatures). The ‘variable’ version of rmcr runs faster than the ‘fixed’ version if an image
contains many bright pixels, but usually it detects less CR signatures and had
more problems to run to completion. For some (not all) images considered with
the variable version, the number of detected CR signatures at medima=16 is
greater than at medima=7,
and smaller than at medima=32,
but these differences are not considerable. The same value of medima was used
for replacement of the brightness of pixels considered as CR signatures (note
that for imgcleanf
we used medima=7).
The code rmcr finds groups of pixels brighter than lim that are
located close to each other. Pixels are considered to belong to one group if (dx+1)2+(dy+1)2<ddilim, where dx and dy are differences in coordinates x
and y of two pixels (the width of a
pixel is equal to 1). The default value of ddilim is equal to 8.1, i.e., two
pixels separated by not more than one pixel belong to the same group.
The code makes statistics for sizes of all groups.
Sometimes rmcr (usually its simple
version trmcr
http://www.astro.umd.edu/~ipatov/trmcr.pro
intended only for statistics) was applied to pixels considered by other codes
as belonging to CR events in order to get the distribution of the events over
the number of pixels in one event.
Rmcr
removes all clusters with kp<0.17, where kp=npix/(dxx2+dyy2), dxx and dyy are the maximum differences of coordinates x and y
in a charge cluster each increased by 1, and npix is the total
number of pixels in the cluster. These clusters correspond to ‘long’ CR
signatures. Charge clusters consisting of not more than nlimit pixels are also considered
as CR signatures. Depending on a considered image and problem, the input
parameter nlimit
can be chosen to take any value (e.g., 10). In one version of the code, the
clusters (except for ‘long’ clusters) that are closer than dss (e.g., dss=10 pixels) to a defined
rectangle that includes the comet and its coma are not considered as CR
signatures. For small values of nlimit, it may be useful to run rmcr for calibrated images two
times - first with a greater lim, and then with a smaller lim. The code runs slowly when
there are a lot of pixels in all clusters (e.g., a comet occupies a
considerable part of an image) because the code analyzes the entire image at
once, rather than by small portions of the image at a time as do the other
codes.
The rmcr code replaces detected CR signatures with values of
their neighboring pixels in the same way as imgcleanf, and it also does not
consider pixels located close to the edges of the frame as CR events.
On ITS and IR images there are many pixels with high
brightness that appeared at the same places on different images. The
information about such ‘warm’ pixels can be extracted from a calibrated dark or
flat field image. The ‘warm’ pixels were found based on analysis of many
images. For removal of signatures of CRs on ITS and
IR calibrated images (see sections 5-6), we modified imgcleanf (the name of the modified code
is imgcleanw http://www.astro.umd.edu/~ipatov/imgcleanwitq.pro)
in such a way that if the brightness of a ‘warm’ pixel is greater than the
median value of close pixels, then it is replaced by this median value (in our
studies we also used the code rmwarm http://www.astro.umd.edu/~ipatov/rmwarm.pro
that only makes the above replacement without searching for CR signatures). The
search of CR signatures was made only for such modified ITS and IR images.
In our studies discussed in section 6, we used also
our code checkcr http://www.astro.umd.edu/~ipatov/checkcr.pro
that finds CR signatures based on two dark images and on bad pixel maps (these
maps were obtained by the DI team and can be found at http://pdssbn.astro.umd.edu/holdings/,
e.g., that for ITS images is located at
http://pdssbn.astro.umd.edu/holdings/dii-c-its-3_4-9p-encounter-v2.0/dataset.html). A cluster of pixels is considered by checkcr to be a
signature of a CR if it is present in only one image in succession and is
located far from ‘bad’ pixels. We considered that this cluster must be
separated by at least one pixel from any cluster on a subsequent image and from
bad pixels. The distance criteria (similar to ddilim) is
an input parameter and can be changed. Before using checkcr, we applied our code trmcr http://www.astro.umd.edu/~ipatov/trmcr.pro
(a simple version of rmcr) to the .fit files that contained
information about CR signatures detected by imgcleanf and to the file with
bad pixels in order to find groups of pixels that are separated by a distance less
than that corresponding to ddilim. Coordinates
of all pixels in such groups (including the number of the group to which pixels
belong) are written by rmcr
in the file delete.txt, and the main statistics are in the file groupk.txt. Note that a cluster can consist of a
different number of pixels on two images, but it can be considered as the same
cluster if it satisfies the distance criterion. The interactive code imr http://www.astro.umd.edu/~ipatov/imr.pro
for removal of CR signatures is discussed in section 7 (sub-version imrr http://www.astro.umd.edu/~ipatov/imrr.pro
produces not only an output array, as imr, but also a few .fit files) .
The codes described above
can be found at our website, starting search from http://www.astro.umd.edu/~ipatov/CRcodes.html.
Details of how to use the codes are at the beginning of the codes.
7. Interactive code for
removal of cosmic ray signatures
We developed
the interactive code imr http://www.astro.umd.edu/~ipatov/imr.pro, which allows one to control manually
the removal of CR signatures on images (sub-version imrr http://www.astro.umd.edu/~ipatov/imrr.pro
produces not only an output array, as imr, but also a few .fit files). The first step of the
interactive code is to replace the brightness of ‘warm’ pixels by the local
median brightness (this step can be deleted for MRI and HRI images). For ITS
images, ‘warm’ pixels can provide a great number of false detected CR
signatures. The code can show the clusters that imgclean detects as CR events (imgclean can be
replaced by any other code). It restores the edges of an image after imgclean (see
section 2). Based on the exposure time and size of the image, the code
calculates the expected number of CR signatures.
The code imr allows a user
to choose the rectangles “A” on a considered image where glitches recognized by
imgclean as
CR signatures are ignored. In other rectangles “B,” which can be chosen by the
user, the brightness of some pixels is replaced by the local
median brightness if the brightness of these pixels is greater by some
factor kdif
(the default value of kdif
is 2) than the median brightness.
The rectangles are
marked by clicking a mouse at two apices. If one knows in advance the
coordinates of corrected regions, he can print the coordinates. In the latter
case it is possible to work with images without seen them (i.e., for remote
access to an IDL computer).
The result of the work
of the code is the following: (1) inside rectangles “A”, an initial image
(i.e., before imgclean)
remains unchanged; (2) inside rectangles “B”, the brightness is changed for
clusters detected by imgclean
as CR signatures and also for other pixels whose brightness differes
by a factor of kdif
from the local median brightness; (3) outside rectangles “A” and “B”, the
brightness is changed for all the clusters detected by imgclean as CR signatures.
Imgclean often detects false CR
signatures near the edge of a comet nucleus. So, one can choose several
rectangles “A” which include such false CR signatures. These rectangles can
occupy the edges of the nucleus, or one can take one larger rectangle that
includes the entire image of the nucleus and its edges. On some ITS images,
there are regions where the brightness of nearby background pixels varies substantially,
and therefore imgclean
detects a great number of false CR signatures in these regions. Rectangles “A”
can also include such regions. We recommend using rectangles “B” in order to
delete those glitches of CRs that were not entirely
deleted by imgclean
(often imgclean
can not delete all pixels of a long CR event). Also one may want to delete
stars. In the above 'long' way in imr, one can use any number of rectangles and can make many
corrections of the image. This version is recommended if one wants to edit an
image carefully.
A 'quick' way in the
code allows a user to choose a rectangle “C” inside which clusters recognized
by imgclean
as CR signatures are not considered as CR events (i.e., inside this rectangle,
there are no changes compared with the original image). Outside this rectangle,
the brightness of all pixels which differed from the median value of local
brightness by a factor of more than kdif is replaced by the median value. Inside rectangle “C”
we make the same calculations as for the rectangles “A”, and outside the
rectangle “C”, we do the same as for rectangles “B”. It is suggested to choose
a rectangle “C” so that it will include the entire image of a nucleus and some
region near the edges of the nucleus (where usually imgclean detects false CR signatures; for images with outbursts this region
can be greater).
The 'quick' way is
recommended when one needs to consider quickly a large number of images. For
example, it can be used for small frames with a small image of a nucleus. Outside
the rectangle “C”, all stars are also deleted. For frames without a comet,
there is no sense to use a 'quick' way in imr.
The
interactive code allows a user to delete those long CR signatures that are not
recognized by imgclean
and to prevent consideration of the pixels close to the edge of a nucleus as CR
events, i.e., to solve two main problems that prevent obtaining good final MRI
and HRI images with imgclean.
Imr can be used for removal of long oblique CR signatures (vertical and horizontal
segments are not caused by CRs) on ITS and IR images. At small exposure times, the probability of long CR
signatures is small. It may be possible to check only ITS and IR images with
exposure time of at least a few seconds and to find those images where long CR
signatures are present. For such images one can delete only long CR signatures
(first take a rectangle “A” to be equal to the whole image in order to not
consider the CR signatures found by imgclean; then put rectangles “B” around long CR signatures
in order to delete these signatures).
//////////////////////////////////////////////////
Examples how to use codes:
raw = readfits(file,h)
di_crrej,raw,out1,/mask,header=h,num=num
dn
= raw/sxpar(h,'radcalv')*sxpar(h,'inttime')/1000.
out2 = crfind(dn)
d = raw
imgclean,d
out3 = d ne ra
grmcr,ima,jma,header=h,mode=2,klim=30
grmcr,ima,jma,header=h,klim=30,n246=181
gjrmcr,ima,jma,mode=2,klim=10.,nlim=10,medima=16
gjfrmcr,ima,jma,mode=2,klim=10.,nlim=10
imgcleanf,ima,jma,header=h
out3=ima
ne jma
fits_write,’h901-3.fit’,out3,h