Combining Multiple Exposures

Using MSCSTACK to Generate Your Initial Stack

WARNING: THIS STEP IN OUR REDUCTIONS IS CURRENTLY DISCUSSED ONLY IN THE VERY SPECIFIC EXAMPLE OF OUR NDWFS REDUCTIONS AND RELIES ON SOFTWARE STILL UNDER ACTIVE DEVELOPMENT AND CHANGE. THESE NOTES WILL BE EXPANDED AND UPDATED ONCE THE RATE OF CHANGE IN THESE PROCEDURES SLOWS DOWN

We use the task mscstack to combine the images for a field into one deep exposure. Note that at the present time we do not make any effort to match the DIQ of the input images, although some groups using MOSAIC have successfully done this. We do, however, take some pains to enable the deepest possible final image, free of defects like faint residual images from satellite trails in an individual image. This can be accomplished if all such defects are masked out of the individual input images. In order to include as many defects in the input masks as possible, we use software written by Frank Valdes, which is still in active development, to compare each of the individual images to an initial average image that is created using a harsh rejection criteria. The software used to make the comparison, part of the ACE package Frank Valdes is writing, also produces a new bad-pixel mask that can be used in our "final" stacking of the images.

To make the "first pass" combined image that will be used to identify defects in the individual exposures, we use mscstack with a low threshold for rejection.


PACKAGE = mscred
   TASK = mscstack

input   =          @infile.list  List of images to combine
output  =       Stackharsh.fits  Output image
(plfile =        Stackharsh.pl) List of output pixel list files (optional)

(combine=              average) Type of combine operation (median|average)
(reject =              ccdclip) Type of rejection
(masktyp=            goodvalue) Mask type
(maskval=                   0.) Mask value
(blank  =               50000.) Value if there are no pixels

(scale  =            !mscscale) Image scaling
(zero   =             !msczero) Image zero point offset
(weight =                 none) Image weights
(statsec=                     ) Image section for computing statistics

(lthresh=                   1.) Lower threshold
(hthresh=                INDEF) Upper threshold
(nlow   =                    1) minmax: Number of low pixels to reject
(nhigh  =                    1) minmax: Number of high pixels to reject
(nkeep  =                    1) Minimum to keep (pos) or maximum to reject (neg)
(mclip  =                  yes) Use median in sigma clipping algorithms?
(lsigma =                   5.) Lower sigma clipping factor
(hsigma =                  2.9) Upper sigma clipping factor
(rdnoise=              rdnoise) ccdclip: CCD readout noise (electrons)
(gain   =                 gain) ccdclip: CCD gain (electrons/DN)
(snoise =                   0.) ccdclip: Sensitivity noise (fraction)
(sigscal=                  0.1) Tolerance for sigma clipping scaling corrections
(pclip  =                 -0.5) pclip: Percentile clipping parameter
(grow   =                   0.) Radius (pixels) for neighbor rejection
(mode   = 

Before diffdetect, part of ACE is run to identify differences between one of the individual exposures and the average or stacked image, mask output by mscstack must be modified. The mask that was created by mscstack is just a record of the number of input exposures that did not contribute to each pixel in the output image. The value for each pixel in the mask is an integer from 0, if that given pixel was contributed to by every input frame, to n, where n is the number of images in the stack whose data were excluded from the being used in the determination of the average. The task diffdetect needs to know which pixels have little or no information (like the positions of saturated stars, which in general will have been saturated in all of the input images). These positions should not be flagged further. For this to work, it is necessary to modify the bpm using imexpr. Here is the command.

 
PACKAGE = imutil 
TASK = imexpr

	expr    =    (a != 20 ? 0 : a)  expression
	output  =           junkmask.pl  output image
	(dims   =                 auto) output image dimensions
	(intype =                 auto) minimum type for input operands
	(outtype=                 auto) output image pixel datatype
	(refim  =                 auto) reference image for wcs etc
	(bwidth =                    0) boundary extension width
	(btype  =              nearest) boundary extension type
	(bpixval=                   0.) boundary pixel value
	(rangech=                  yes) perform range checking
	(verbose=                  yes) print informative messages
	(exprdb =                 none) expression database
	(lastout=   _blmask0_obj169s_8) last output image
	a       =        Stackharsh.pl  operand a
	z       =                       operand z
	(mode   =                   ql)
Note in this case, 20 images were combined to make the final stack.

We next run diffdetect on each individual image with this new mask. diffdetect compares each image with the stack and notices what is different taking into account the noise properties of both the single exposure and the combined stack. An example parameter file.


PACKAGE = ace
   TASK = diffdetect

images  =              mos068S  List of images
(masks  =                  BPM) List of bad pixel masks
(skys   =               SKYFIT) List of skys
(sigmas =               SKYSIG) List of sky sigmas
(scales =            !mscscale) List of image intensity scale factors

                        # Reference Image(s)
(rimages=      Stackharsh.fits) List of reference images
(rmasks =         junkmask4.pl) List of reference bad pixel masks
(rskys  =               SKYFIT) List of reference skys
(rsigmas=               SKYSIG) List of reference sky sigmas
(rscales=                   1.) List of reference intensity scale factors

                        # Output
objmasks=             diff068S  List of output object masks
catalogs=             diff068S  List of output catalogs
(logfile=               STDOUT) List of log files

                        # Sky
(newsky =                  yes) Determine new sky fit if one already exists?
(nskylin=                  100) Number of sky sample lines
(skyblk1=                   10) Sky block size for 1D sky estimation
(skyhcli=                   3.) High sky clipping during 1D sky estimation
(skylcli=                   3.) Low sky clippling during 1D sky estimation
(skyxord=                    4) Sky fitting x order
(skyyord=                    4) Sky fitting y order
(skyxter=                 half) Sky fitting y order

                        # Detection
(convolv=            block 3 3) Convolution kernel
(hsigma =                   3.) Sigma threshold above sky
(lsigma =                  10.) Sigma threshold below sky
(minpix =                   10) Minimum number of pixels in detected objects
(sigavg =                   6.) Sigma of mean flux cutoff
(sigmax =                   6.) Sigma of maximum pixel
(rfrac  =                  0.5) Minimum fraction of reference flux
(bpval  =                    1) Output bad pixel value

                        # Growing
(ngrow  =                    2) Number of grow rings
(mode   =                   ql)
This will generate two files, a list of objects in a STScI tables format, and a bad pixel mask.

We use acedisplay and reviewproto to examine the image and make sure that 1) no good areas of the image are included in the new mask by some error, and 2) that everything we want in the mask has been caught. First look at the stack and make sure there are no surviving satallite trails or scattered light. If there are, you might have to go through the individual images and mask these out by hand in the new ace-generated bad pixel mask using satzap for satallite trails and blkzap for bad areas in general (and tasks to "undo" corrections, unsatzap and unblkzap. The source for these scripts is available at these links: satzap.cl by James Rhoads and Arjun Dey
blkzap.cl by Arjun Dey, based on script by James Rhoads
and (unblkzap.cl and (unsatzap.cl).

Here are parameter files for acedisplay and reviewproto.


PACKAGE = ace
   TASK = acedisplay

image   =              mos068S  image to be displayed
frame   =                    3  frame to be written into
(bpmask =                  BPM) bad pixel mask
(bpdispl=                 none) bad pixel display (none|overlay|interpolate)
(bpcolor=                  red) bad pixel colors
(overlay=             !objmask) overlay mask
(ocolors=           1=red,+203) overlay colors
(erase  =                  yes) erase frame
(border_=                   no) erase unfilled area of window
(select_=                  yes) display frame being loaded
(repeat =                   no) repeat previous display parameters
(fill   =                  yes) scale image to fit display window
(zscale =                  yes) display range of greylevels near median
(contras=                 0.25) contrast adjustment for zscale algorithm
(zrange =                  yes) display full image intensity range
(zmask  =                     ) sample mask
(nsample=                 1000) maximum number of sample pixels to use
(xcenter=                  0.5) display window horizontal center
(ycenter=                  0.5) display window vertical center
(xsize  =                   1.) display window horizontal size
(ysize  =                   1.) display window vertical size
(xmag   =                   1.) display window horizontal magnification
(ymag   =                   1.) display window vertical magnification
(order  =                    0) spatial interpolator order (0=replicate, 1=linear)
(z1     =                     ) minimum greylevel to be displayed
(z2     =                     ) maximum greylevel to be displayed
(ztrans =               linear) greylevel transformation (linear|log|none|user)
(lutfile=                     ) file containing user defined look up table
(mode   =                   ql)

PACKAGE = ace
   TASK = reviewproto

catalog =         diff051S.tab  Catalog
(nooverl=                  yes) Display image without overlays
(overlay=                  yes) Display image with overlays
(compari=                  yes) Display comparison image
(compima=    11q1I-bandV7.fits) Comparison image
(box    =                  200) Box size (pixels)
(ocolors=                  204) Object mask color
(lcolor =                  205) Label color
(fd     =                     )
(mode   =                   ql)

PACKAGE = user
   TASK = blkzap

object  =              mos065S  Input image
x1      =                1247.  lower left corner: location 1 - X 
y1      =                3705.  lower left corner: location 1 - Y 
x2      =                1851.  upper right corner: location 2 - X 
y2      =                3787.  upper right corner: location 2 - Y 
(mode   =                   al)

PACKAGE = user
   TASK = satzap

object  =              mos061S  Input image
x1      =                1151.  Satellite trail: location 1 - X 
y1      =                3968.  Satellite trail: location 1 - Y 
x2      =                1596.  Satellite trail: location 2 - X 
y2      =                4011.  Satellite trail: location 2 - Y 
(hwidth =                  10.) half-width of flagged track?
(mode   =                   al)

Now the final stack can be made! This is done using mscstack again, but with a much less harsh rejection setting, counting on our improved bad-pixel masks to reject the contribution of unwanted pixels.


PACKAGE = mscred
   TASK = mscstack

input   =         @infile.list  List of images to combine
output  =           11q1I-band  Output image
(plfile =    11q1I-band_bpm.pl) List of output pixel list files (optional)

(combine=              average) Type of combine operation (median|average)
(reject =              ccdclip) Type of rejection
(masktyp=            goodvalue) Mask type
(maskval=                   0.) Mask value
(blank  =               50000.) Value if there are no pixels

(scale  =            !mscscale) Image scaling
(zero   =             !msczero) Image zero point offset
(weight =                 none) Image weights
(statsec=                     ) Image section for computing statistics

(lthresh=                   1.) Lower threshold
(hthresh=                INDEF) Upper threshold
(nlow   =                    1) minmax: Number of low pixels to reject
(nhigh  =                    1) minmax: Number of high pixels to reject
(nkeep  =                    1) Minimum to keep (pos) or maximum to reject (neg)
(mclip  =                  yes) Use median in sigma clipping algorithms?
(lsigma =                   5.) Lower sigma clipping factor
(hsigma =                   5.) Upper sigma clipping factor
(rdnoise=              rdnoise) ccdclip: CCD readout noise (electrons)
(gain   =                 gain) ccdclip: CCD gain (electrons/DN)
(snoise =                   0.) ccdclip: Sensitivity noise (fraction)
(sigscal=                  0.1) Tolerance for sigma clipping scaling corrections
(pclip  =                 -0.5) pclip: Percentile clipping parameter
(grow   =                   0.) Radius (pixels) for neighbor rejection
(mode   =                   ql)

Thus ends the "basic" reduction of the images. Next it is time to do the photometric calibrations and generate catalogues from these images.