WIRE Image Co-Addition Software: wkoadda.f

Preface

The FORTRAN software wkoadda.f co-adds a user-specified number of WIRE (Wide-Field Infrared Explorer) images of the same field, which may be from either the WIRE simulator or the actual mission itself. The software is commonly refered to as the "coadder." This document describes the WIRE coadder and how to use it.

Overview

The main purpose of image co-addition is to enhance detector sensitivity through noise migitation. Co-addition of sub-pixel dithered images also greatly improves the quality of undersampled images, which is the case for the WIRE 12-µm band.

For the simplest case of a number of different identically-registered, defect-free images of the same field, where each image includes inherent photon noise as its dominant noise source, image co-addition is simply a sum of the image data divided by the number of images in the sum on a pixel-by-pixel basis. Since photon noise is additive, zero-mean, sample independent, and has a variance that is proportional to the underlying image intensity, it can only be reduced by image co-addition. The resulting co-add image will have a reduced photon noise variance that is lower by one over the number of input images in the co-addition. Other types of additive, zero-mean, sample-independent noise present in the data will be similarly reduced.

The coadder performs a number of processing steps, including image scaling and corrections for optical distortion, image registration (translations and rotations), interpolation, up-sampling, bad-pixel handling, tracking of sample distributions, image-statistics calculations, outlier rejection, co-addition, and normalization.

The interpolation algorithm was designed to preserve image photometry. An area-weighted interpolation is employed, which assures that all data number (DN) counts in the input image are co-added into the output co-add image, to a high degree of accuracy.

The same type of interpolation is used to produce an output "coverage" image, which stores the number of input images contributing to each pixel in the output co-add image. The coverage image data are used in the normalization step.

Normalization consists of multiplying each coadd pixel interpolate by the ratio of the area of a WIRE pixel to the area of a coad pixel, which is nominally 0.25, and dividing by the corresponding coverage image interpolate.

Image Scaling and Distortion Corrections

Input image data are transformed from scaled/distorted space to undistorted space with a scale consistent with f=1000 mm:

(xDIS, yDIS) = coordinates in scaled/undistorted space

(xUND, yUND) = coordinates in undistorted, f=1000-mm space

(xC, yC) = coordinates of image rotation center

(xDC, yDC) = coordinates of image scaling/distortion center

(xDIS, yDIS) = nonlinear distortion corrections

(xFPA25, yFPA25) = shifts applied to 25-µm data to correct for misalignments between 12-µm and 25-µm FPAs

(fx, fy) = focal lengths of the telescope optics, in mm

Distortion corrections are found by least-squares fits to ray-tracing model data, where fx and fy are the fit parameters and (xDIS, yDIS) are the residuals:

, , (,) = ray field angles

Distortion corrections, xDIS and yDIS, are each in general functions of (xDIS, yDIS) and the infrared band.

Interpolation, Co-Addition, and Normalization

Interpolation step:

Interpolate coadd-image intensity, IOk(X, Y) , and coverage (weight), wOk(X, Y), from at most four pixels in input image k:

Ii = input-image intensity (background-subtracted) at pixel i

wi = input-image coverage (weight) at i

Ai = input-image pixel area (after scaling/distortion correction) at pixel i

fOi = fractional overlap area of coadd pixel at (X, Y) onto input-image pixel i

AO = area of coadd pixel (same for all k)

Co-addition step:

Normalization step:

Outlier Rejection

Radiation in the form of energetic (MeV) protons are expected to have possibly three different kinds of unwanted effects on the WIRE focal plane arrays (FPAs):

  1. charge deposition
  2. dark-current elevation
  3. changes in detector responsivity

The first effect produces typically produces between 10,000 and 100,000 electrons of charge on a detector per proton. The second effect persists over a few readout intervals. Both transient and possibly permanent radiation damage are expected from the third effect. Synchrotron radiation tests conducted by Dr. Terry Herter have found that only the first two effects are significant. Prolonged radiation exposure of the FPAs may show the third effect to also be important, as with IRAS.

Radiation sources are energetic charged particles trapped in the Earth's magnetic field and cosmic rays. Radiation effects during the WIRE mission will be highest when the spacecraft passes through the South Atlantic Anomaly (SAA), a region in the Earth's southern hemisphere of high charged particle concentration. The trapped particles are generally modeled with an isotropic velocity distribution. Passage time through the SAA varys with geographic position, covering as much as a 60 orbital angle range at 320 longitude (this is for IRAS altitudes; for WIRE altitudes the angular range will be smaller). Detailed numerical modeling of the WIRE spacecraft trajectories through the SAA for various survey targets is forthcoming (Henderson, private communication); for now it is assumed that for surveys at SAA orbital angles, 10% of the survey will include strong SAA effects. It has been estimated that about 10% of the pixels in a 56-second WIRE-image exposure taken while passing through the SAA will be affected by rad hits, assuming a flux of approximately 25 protons/cm2/s.

The detectors in the WIRE FPAs have a larger extent in pixel dimension (75 m on a side) than in thickness (10-15 m), and therefore geometry effects are expected to be minimal, with most rad hits isolated in individual detectors. This is in contrast to the much thicker detectors used in the HST, in which contiguous lines of several detectors are affected by a single event.

Because the detector gain settings are 80 electrons per DN for the 12-m band and 40 electrons per DN for the 25-m band, the 12-bit ADC (analog-to-digital converter), which digitizes the detector output signals and has a 0-4095 DN range at its output, saturates at about 330,000 and 160,000 electrons per respective band. A single WIRE image consists of 14 co-added sub-frames (14 separate ADC readouts are hardware co-added), and the probability of the same pixel being hit more than once during the data acquisition of the 14 sub-frames is low. Therefore, when the effect of a single rad hit is spread out over 14 sub-frames, the net relative change in DN of WIRE image due to a rad hit is small when relatively bright sources are imaged, and since the detector gains are set for detection of these bright sources in the linear region of the ADCs response, detector saturation is unlikely. On the other hand, when dim sources are imaged, the effect of a rad hit on a detector is relatively more important.

Thus, rad hits will cause mostly measurements of larger values than would be otherwise obtained in the absence of this radiation. Detection and correction of these outlier data are necessary to remove the positive bias in the data caused by rad hits. It is natural that this processing be done in the coadder, since it involves using statistics that are derived from the same input data included in the image co-addition.

Outlier rejection is not just limited to positive outliers caused by rad hits. Positive outliers caused by either photon or rad-hit noise are indistinguishable; both types will be rejected. There are also negative outliers due to large negative photon noise contributions, which can also be rejected. Throwing away both positive and negative outliers (that pass the exceedance test) may minimize the overshoot in the correcting bias caused by truncation of the nonsymmetric distribution of photon plus rad-hit noise.

Based on the above considerations, the following algorithm for outlier rejection has been implemented in the coadder:

  1. Read from command-line input or namelist file the average background level (in units of DN) over all input images to be included in the coadd image. This constant is assumed to be both spatially and temporally constant.


  2. Compute the up-sampled coadd and coverage images from the background-subtracted samples, with no outlier rejection. These images are pixel-to-pixel variations of the mean and "number" of the background-subtracted samples used to compute the mean. Also compute the sample-noise (standard deviation) and sample-count images. During this process store the 10 largest and 10 smallest values of these background-subtracted samples and also their associated coverage, sample-noise and sample-count values, for each pixel; these are the candidate outliers for the rejection algorithm. This allows for up to 3.6x106 positive outliers to be rejected (and also an equal number of negative outliers), assuming 600x600-pixel output images. An ultra-deep survey for the no evolution assumption requires the co-addition of 3200 images; assuming 5% of the time is spent in the SAA, about 3.6x106 rad hits equally distributed over the FPAs are expected for the flux value quoted above (radiation shielding issues aside). The available computer memory limits the number of candidate outliers that can be stored and tested (see next section for a workaround to this limitation).


  3. Compute a map of the photon-noise standard deviation, This is given by



    where is the intensity at coadd-image pixel ij plus the background, and



  4. Reject those outliers that meet the criteria listed below, which depend on the type of outlier rejection method specified by the user (as defined below). Regardless of the outlier rejection method chosen, outlier rejection will not occur at coadd pixel locations where the coverage is insufficient; the required minimum coverage amount is user-specified by the command-line option -ollim or the namelist variable nmincov (see Table 1 below). The type of outlier rejection method is user-specified by the variable rsigma, which is overloaded to handle threshold values in different units depending on whether it is negative, lies in the [0,1) range, or is greater than 1. Thus, the value of rsigma determines which of the three different outlier rejection methods available is used by the code. The three outlier rejection methods are as follows:


Each outlier meeting the rejection criteria is removed from the coadd intensity image, and its coverage amount is removed from the coadd coverage image.

Positive outliers are rejected first, before negative outliers, until a user-specifiable minimum coverage criterion is met (see entry for -ollim and nmincov in Table 1 below).

Sofware Execution

A Makefile and located in the same directory as the source code, for building an executable file called wkoadda from the source code.

The image coadder is typically executed in the WIRE data analysis (DA) pipeline by a perl script. DA pipeline software executed prior to wkoadda creates all necessary input files.

The required format for coadder input images is FITS (Flexible Image Transport System); the output images created by the coadder are also in this format.

The coadder is executed as follows:

%wkoadda [command-line options] offsetfile

The offset file contains information about which input images are to be coadded along with relevant image registration data.

The coadder can also be executed with a namelist file as an alternative means to specifying command-line options:

%wkoadda namelistfile offsetfile

Required and Optional Inputs

Table 1. Coadder command-line options. Note that the corresponding namelist variable name is given in parantheses if different from the command-line option name.

Command-line option (namelist variable) Definition Default value(s)
maximagesize Maximum size of output coadd image, in terms of number of pixels on a side. Default is 1024 pixels. 1
band Integer flag indicating the WIRE infrared band of the data to be co-added. Set to 1 for the 12-µm band and 2 for the 25-µm band. 1
nrep The width (or height) of a WIRE pixel in units of coadd pixels. 2; in terms of areal coverage, this value is used to specify 4 coad pixels per WIRE pixel
OPlines, OPcols (nlino, ncolo) The number of pixels per row and per column in the coadder's output images. Used to size output coadd and coverage images. If either value is 0, then these quantities and xdel and ydel (see below) are recalculated. 600, 600; these values are used to size the memory allocation for each output image.
xOff, yOff (xdel, ydel) Relative x and y shifts, in coadd pixels, of the first WIRE image relative to the lower left-hand corner of the coadd image (positive values mean shift right and up). Used to center WIRE images in coadd image. If either value is -9999, then these quantities and nlino and ncolo are recalculated -9999.9,-9999.9
Nozap (Use complement zap)Logical variable to control blanking of "empty" coadd pixels. If Nozap is set to .false., then the blanking value is NaN; otherwise the blanking value is zero. .false.
noNormalize (Use complement norm)If this switch is set to .true., then do not normalize the coadd intensity image by the coadd coverage image. .false.
minCov (ncmin)If the coverage value of a given coadd pixel is ncmin, then that pixel is declared "empty" and blanked out. 0.0
readCov (readn)Logical flag to indicate that input coverage image data are to be used. If readn=.true., then nifile (see below) must be specified. .false.
(nifile) Path and filename of the input coverage image(s). The path and filename of the corresponding input image with the string nprefix appended to the filename's beginning.
covPrefix (nprefix) Filename prefix of the input coverage image(s). "n"
noOP (noout)Logical variable to suppress creation of output coadd intensity and coverage images, and rejected-outlier images if set to .true. .false.
ncoutLogical variable to allow output of coadd coverage image if set to .true. (noout must also be set to .false.). .true.
noiLogical variable to allow output of coadd noise image (standard deviations of coadd samples) if set to .true. (noout must also be set to .false.). .false.
ofilePath and filename of the coadd image. The path and filename of the first input image with a "c" appended to the filename's beginning.
nfilePath and filename of the output coverage image. The path and filename of the coadd image with an "n" appended to the filename's beginning.
focallenx, focalleny The x and y focal lengths (in mm) of the WIRE telescope that determine the x and y image scales. Use 1000.0 mm for simulated WIRE images. focallenx=999.752 mm, focalleny=994.438 mm for the 12-µm band; and focallenx=999.899 mm, focalleny=994.572 mm for the 25-µm band (as derived from Roy Esplin's ray-tracing model data for the WIRE system).
distortflag Logical variable to indicate whether or not corrections for distortion (nonlinear image scaling) are to be applied. Distortion corrections are computed from distortion model fit parameters in user-supplied subroutine distortf() in file distortion.f. The fit parameters are potentially band-dependent, and are derived either from measurements (such as during IOC) or from ray-tracing model data. .false. (currently Roy Esplin's ray-tracing model data show that distortion is negligible for the WIRE design).
xdc, ydc The x and y coordinates of the "distortion center" of the input images in WIRE pixel units (these are also potentially band-dependent). The distortion center is simply the origin of the coordinate system in which the scaling/distortion corrections are defined. The default coordinate values are (64.5, 64.5), which is the FPA center.
fpadtheta, fpadx, fpady The rotation (degrees) and x and y shifts (WIRE pixels) required to co-register a 25-µm image to a corresponding 12-µm image (positive fpadtheta,fpadx,and fpady mean rotate counter-clockwise, and shift right and up the 25 µm image). 0., 0., 0.
ebugLogical flag to cause activation of all IEEE signal handlers. .false.
forcfo Logical flag to force output of floating-point images. .true.
nmissingok Logical flag to indicate whether to continue the processing with the value 1 for input coverage in the event that the input coverage image does not exist. .false.
survfrtime Exposure time of a WIRE input image, in seconds. 56 seconds
background Average background level over all input images included in coadd image, assumed to be spatially and temporally constant. Read in from the input image FITS header if not given here. 0 DN
rsigma Outlier-rejection control parameter, overloaded as follows:
  • If rsigma is negative, then candidate outliers are rejected if they lie outside of the N-sigma envelope of the photon-noise distribution with its mean given by the value of the associated background-subtracted sample, where
    N =|rsigma|.
  • If rsigma=0, then no outlier rejection will be done, and substantially less RAM will be required.
  • If 0<rsigma<1, then the number of candidate positive outliers that are rejected is equal to the fraction rsigma of total number of samples in the coadd image at the pixel location of interest, up to the available number of candidate positive outliers; candidate negative outliers are handled similarly.
  • If rsigma>=1, then rsigma is the number of candidate positive outliers that are rejected at the coadd-image pixel location of interest, up to the available number of candidate positive outliers; candidate negative outliers are handled similarly.
0
staredgthresh Star edge threshold. A bright source is considered to be present when the ratio of sample noise to photon noise exceeds this threshold. 2
ollim (nmincov) The minimum coverage amount required at a given coadd-image pixel location for outlier rejection to occur. 20 samples or units of exposure time relative to the frame time
rejprefix Filename prefix of the rejected-outlier images. "ro"
rejimages A logical array of 4 elements to indicate whether (.true.) or not (.false.) the following types of rejected-outlier images are to be generated:
  • input coverage (ic)
  • input intensity and frame (ii and if)
  • coadd coverage (cc)
  • coadd intensity image (ci)
.false., .false., .false., .false.

Namelist file

The namelist file contains a list of namelist variables for parameter settings that override the default values that are set in the coadder. The available namelist variables are defined in Table 1 above. A sample namelist file is given as follows:

 $in
    nrep = 4,   ncout = t,
    zap = t, readn = f,
    ncmin = 0.0, nono = f,
    xdel = -9999, ydel = -9999, 
    ofile = "cws224700103.fits",
    nprefix = "n",
    fcg = f,   bik = f,    coff = f,
    focallenx = 1000., focalleny = 1000., distortflag = f,
    survfrtime = 56., rejprefix = "ro", rejimages = t,t,t,t,
    background = 8350., rsigma = 2., mincov = 1,
 $

The namelist variables that are omitted from the namelist file will be set to their default values.

Offset file

The offset file is generated by other codes in the WIRE data analysis (DA) pipeline, and contains, among other things, a list of the filenames of the input images, their sizes (x-pixels by y-pixels), image-intensity offset (in DN or data number), and image registration data (x-shift and y-shift in pixels, rotation angle in degrees, pixels coordinates of the rotation center). A sample offset file is given as follows:

|frame                                                                          |ncol|nlin|nmat|  xcen|  ycen|    zoff|   theta|    xdel|    ydel| xmin| xmax| ymin| ymax| xsig| ysig|
| c										|  i |  i |  i |    r |    r |      r |      r |      r |      r |   r |   r |   r |   r |   r |   r |
/stage/wire-pit/tim/pipe/reduce//sws301500103                                     128  128   18   0.00   0.00    0.000    0.000    0.000    0.000     0   129     0   129 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500104                                     128  128   18   0.00   0.00    0.000    0.040   10.430    3.740    10   140     4   133 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500105                                     128  128   18   0.00   0.00    0.000    0.080    7.260   16.470     7   136    16   145 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500106                                     128  128   17   0.00   0.00    0.000    0.100  -16.870   -5.910   -17   112    -6   123 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500107                                     128  128   17   0.00   0.00    0.000    0.080   -3.020   14.970    -3   126    15   144 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500108                                     128  128   19   0.00   0.00    0.000    0.180   -2.760    6.360    -3   127     6   135 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500109                                     128  128   19   0.00   0.00    0.000    0.020   -1.210   17.390    -1   128    17   146 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500110                                     128  128   19   0.00   0.00    0.000    0.050   -0.580   -1.900    -1   129    -2   127 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500203                                     128  128   24   0.00   0.00    0.000   -0.865   14.753  -16.336    13   144   -16   115 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500204                                     128  128   24   0.00   0.00    0.000   -0.795    9.167  -16.650     7   138   -17   114 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500205                                     128  128   25   0.00   0.00    0.000   -0.855    5.562    6.868     4   135     7   138 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500206                                     128  128   21   0.00   0.00    0.000   -0.925  -15.780    9.466   -18   113     9   141 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500207                                     128  128   22   0.00   0.00    0.000   -0.855  -22.463  -18.228   -24   107   -18   113 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500208                                     128  128   27   0.00   0.00    0.000   -0.905   10.943  -11.723     9   140   -12   119 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500209                                     128  128   27   0.00   0.00    0.000   -0.845   14.369   -7.460    12   143    -7   123 0.000 0.000
/stage/wire-pit/tim/pipe/reduce//sws301500210                                     128  128   20   0.00   0.00    0.000   -0.895  -11.685   16.439   -14   117    16   147 0.000 0.000
\#SUMRY: #frames= 160  master= 1    maxdx= -24.1  maxdy= -24.2
maxdtheta=-2.2 

Outputs

The outputs from the coadder are primarily in the form of FITS image files. Below are the types of FITS files that can be generated, along with example filenames.

The coadd intensity, coverage, and sample noise images are corrected for any outliers rejected (their effects are removed).

The generation of one or more output images depends on the command-line or namelist variable settings, and coadder outputs as well. Rejected-outlier images are generated if four conditions are satisfied: 1) noout=.true., 2) rsigma is nonzero, 3) the number of rejected outliers found is nonzero, and 4) rejimages()=.true.

Description of rejected-outlier images

The coadd photon noise image is given for comparison with the coadd sample noise image. The ratio of the latter to former is used to detect the presence of bright sources, by looking for threshold exceedances (see command-line option staredgethresh). Outlier candidates found on bright sources are subject to a more stringent rejection test; this is necessary to prevent throwing away bright-source samples that have higher variability due to large pixelization relative to PSF width.

The rejected-outlier coadd intensity image contains the sum of all the outliers rejected, both positive and negative. This image is normalized by the rejected-outlier coadd coverage image if nono is set to .false.

The rejected-outlier coadd coverage image contains the coverage of all the outliers rejected, both positive and negative.

The rejected-outlier input images store only information about positive rejected outliers. The rejected-outlier input image contains the difference between the biggest outlier rejected and the coadd intensity (after removal of rejected outliers) at a give pixel location; this quantity is normalized by the sample noise if nono is set to .false.

Examples

To run the coadder in the two-pass outlier rejection mode, first do the first pass:

/proj/wire/dev/bin/wcoad -v -xp /proj/wire/russ/pipe/coad/ -of cws450000100.of -co test1stpass-a_s.fits -nr 2 -b 1 -:wkoadda '-noi -nrep 2'

This produces the following coadd intensity, coverage, sample-noise, and sample-count images before outlier rejection:

test1stpass-a_s.fits
ntest1stpass-a_s.fits
cnoi-test1stpass-a_s.fits
nsco-test1stpass-a_s.fits

Now do the second pass:

/proj/wire/russ/pipe/coad/wkoadda2ndpass -df cws450000100.of -o test2ndpass-a_s.fits -p n -band 1 -fb -nr 2 -xdel -9999.9 -ydel -9999.9 -noi -stedg 3 -mis 400 -ollim 4 -sig -8 -bsig -10 -bg 7500 -reji t,t,t,t -debug -fpci test1stpass-a_s.fits -fpcc ntest1stpass-a_s.fits -fpcs cnoi-test1stpass-a_s.fits -fpcn nsco-test1stpass-a_s.fits

This produces the following coadd intensity, coverage, sample-noise, and sample-count images after outlier rejection, as well as the other coadder image products that describe the rad hits found and rejected:

test2ndpass-a_s.fits
ntest2ndpass-a_s.fits
cnoi-test2ndpass-a_s.fits
nsco-test2ndpass-a_s.fits
pnoi-test2ndpass-a_s.fits
roci-test2ndpass-a_s.fits
rocn-test2ndpass-a_s.fits
roii-test2ndpass-a_s.fits
roin-test2ndpass-a_s.fits
roif-test2ndpass-a_s.fits
fpcc-test2ndpass-a_s.fits
fpci-test2ndpass-a_s.fits
fpcn-test2ndpass-a_s.fits
fpcs-test2ndpass-a_s.fits
(The last four images in the list are made by enabling the -debug switch in wkoadda2ndpass.)

References

Wells, D. C., E. W. Greisen, and R. H. Harten, FITS: A Flexible Images Transport System, A&AS 44, p. 363, 1981.

Last revised: July 22, 1998

By: Russ Laher (laher@ipac.caltech.edu)

URL: http://spider.ipac.caltech.edu/staff/laher/wkoaddadir/wkoadda3.html

Return to Previous Page