
Streak Removal Utility ('Magic')

This tool analyzes a series of images and detects and remove streaks.
The detection process relies on the difference image from the mean
static sky image compared with the variance static sky image.  The
IPP provided images are reference by a base skycell name and sequence
number.  Three files are important to processing with nominal names:

- <skycell>.dif.<seqNum>.fits - difference image
- <skycell>.dif.<seqNum>.mk.fits - binary mask image
- <skycell>.dif.<seqNum>.wt.fits - variance image

A common base filename to use is <skycell>.dif.<seqNum>.  The command line
interface allows arbitrary file paths and names to be specified in the
--detect stage using:

-i <difference image> or --image <difference image>
-m <8bit mask image> or --mask <8bit mask image>
-w <weight image> or --weight <weight image>
--outroot <output root name used for detection files>

The difference and variances images will be encoded using Rice
compression in FITS format.  The binary mask image will be encoded
using PLIO compression in FITS format.  Currently, the difference
images have the World Coordinate System (WCS) settings stored in a
separate file which can be loaded using --wcs <wcs file>.  (TBR: the
wcs elements for the difference images will be supplied in the image
headers).  The weight or variance images will require scaling by the
CELL.VARFACTOR scale factor so that the variance of a pixel_{ij} in
the weight image is given as:

variance_{ij} = CELL.VARFACTOR * weight_{ij}

For the mask images, the header contains a list of the mask bit values
corresponding to masking concepts.  These are listed using the following keywords:

MSKNUM : an integer containing the number of defined mask bits
MSKNAMnn : a string containing the name of a bit. 'nn' has the value 
  0 - MSKNUM (%02d format)

MSKVALnn : an integer specifying the bit value(s).  this need not be a
  value of 2^n as a single name may represent more than a single set
  bit. 'nn' has the value 0 - MSKNUM (%02d format)

Examples:

MSKNAM01= 'FLAT    '           / Bitmask bit name                              
MSKVAL01=                    4 / Bitmask bit value                             
MSKNAM02= 'DARK    '           / Bitmask bit name                              
MSKVAL02=                   64 / Bitmask bit value                             
MSKNUM  =                   13 / Bitmask bit count                             

(other bits not listed)

To identify the bad pixels (pixels to ignore), the named bits listed
below should be selected and their bit values OR-ed together.  The
resulting value should then be AND-ed with the mask image pixel value.
If the resulting value is zero, the pixel is valid.

Named bits to ignore:

DETECTOR, FLAT, DARK, BLANK, RANGE, SAT, BAD, BAD.WARP, CR, GHOST

Detection is based on @f$ I_{diff}^2 > \tau^2 \sigma_{static}^2 @f$ where
@f$ \tau @f$ is the threshold in sigma.  The mask image is used to define
pixels to ignore.  This static sky mean image, @f$\overline{I_static}@f$
and variance image, @f$\sigma_static^2@f$ are scaled to match the PSF of
the observation image, @f$I_{obs}@f$.  The difference image is defined as
@f$I_{diff} = I_{obs} - \overline{I_static}@f$.

The streak removal code is run in 3 stages, '--detect', '--merge', and
'--remove'.  The detect stage should be run first on each skycell image
separately.  The detect stage generates 3 files for each skycell.

- <outroot>.clusters is a binary file containing all pixel clusters detected
- <outroot>.streaks is a text file containing all streaks detected at the
                     current stage of processing
- <outroot>_hough.fits is a FITS file containing the hough transform
                        computed from the pixel clusters detected

The merge stage should be processed as multiple steps between adjacent
skycells, combining skycells until the entire FOV is processed. The merge
stage updates all of the <outroot>.streaks files with any new streaks
detected.  Iterating from the first stage of 2x2 set of skycells up to the
final stage of full array of skycells will require ~3 merge stages.
The merge stage requires a list of outroot names specified in the detect
stage for the skycells currently merged. Although not currently used in
the merge stage, the list of difference, mask, and weight images along with
optional WCS list can be specified using:

--inputs <input list file>
-I <difference image list file> or --images <difference image list file>
-M  <binary mask image list file> or --masks <mask image list file>
-W <weight image list file> or --weights <weight image list file>
--wcsList <WCS list file>

Each input list file is a text file with a single input file path per line
terminated by '\n'.  The image in the first line of image.list should
correspond to the mask and weight files found in mask.list and weight.list
text files.  The number of images (-I | --images) specified must match
the number of mask (-M | --masks) and weight (-W | --weights) images
specified as well as optional wcs files (--wcsList) if specified.

The remove stage is the final processing stage which takes all streaks
detected and removes them from all source observation images,
including difference images.  For a given difference image on which
streaks are detected, the streaks are removed from the following
associated images: raw images, detrend images, warped images,
difference images.  For a given difference image, there may be
multiple versions of the detrended or warped images which must be
masked.  The remove stage will be called separately for each of the
image files which must be masked.

The remove stage must extend streaks to the sensor Field of View (FOV)
boundaries over all skycells in the image as well as map to the raw
images.  In the raw and detrended images, the pixels which were not
included in any of the streak-processed warps must also be masked.
Note that the warp and difference skycells are only generated if more
than a small fraction of the pixels are lit by the input image.

The remove stage will rely on the astrometry information for a given
type of image to perform the conversion between (ra,dec) and (x,y) as
needed.  For the warp and difference image, the coordinate
transformation information will be stored in the header using
traditional WCS keywords:

CTYPE1,2
CRVAL1,2
CRPIX1,2
CDELT1,2
PC00i00j (i & j = 1 or 2)

For the raw and detrended images, the coordinate transformations are
specified in the associated astrometry file, and use features beyond
the traditional linear WCS transformations to represent non-linear
distortions.  Conversion routines between pixels coordinates and
celestial coordinates (and vice-versa) will be performed using
IPP-supplied functions. 

The IPP higher-order WCS terms use two concepts beyond the traditional
WCS.  First, the additional keywords are defined to specify
the higher order coefficients using the following form:

PCAnXiYj 

where n = 1 or 2 for the output X or output Y direction, and i,j = 0..N
order, representing the coefficient for a term of the form X^i Y^j.
The maximum order is defined by the keyword NPLYTERM.  If higher-order
corrections are applied, the CTYPE1,2 values RA---PLY, DEC--PLY are
used in place of RA---TAN,DEC--TAN.  (NOTE: in future, we will adjust
the code to allow these keywords for any of the traditional WCS
transformation types, testing for their existence.  The current
implementation is somewhat rigid).

In addition, two new CTYPE values are defined which can be used with
higher order terms as well: WRP and DIS (ie RA---WRP,DEC--WRP;
RA---DIS,DEC--DIS).  These CTYPEs are used for two-level astrometry in
which the chip coordinates are transformed to a common focal-plane
coordinate system, which is in turn transformed to a celestial
coordinate system.  The WRP CTYPE defines a transformation which does
not include a projection: only the polynomial coordinate
transformation is applied.  The DIS behaves exactly like the
traditional TAN (or the new PLY projections) but implies the presence
of a correspondnig WRP transformation.  Note that in implementation,
there will be a single DIS transformation for a mosaic-camera
exposure, and one WRP transormation for each chip in the mosaic.

All removed pixels will be stored in a compressed FITS files based on
output_root, image, mask, and weight names and appended with
_archive.fits.  The archive files will be stored in a <archive>
directory accessible only by the IPP.  This will allow restoration of
raw and processed pixels as required.

An example streak removal script for 60 skycells would be:

# Detect stage which can be run at up to processors equal to the number
# of skycells (~60 processors)
RemoveStreak --detect --image image00.fits --mask mask00.fits \
             --weight weight00.fits --outroot output_root00
...
RemoveStreak --detect --image image76.fits --mask mask76.fits \
             --weight weight76.fits --outroot output_root76

# Merge stage (executed in 3 stages up to the full frame of skycells)
RemoveStreak --merge --inputs input.list --images image0_1.list \
             --masks mask0_1.list --weights weights0_1.list
...
RemoveStreak --merge --inputs input.list --images image6_7.list \
             --masks mask6_7.list --weights weights6_7.list
# Second stage merge
RemoveStreak --merge --inputs input.list --images image0_3.list \
             --masks mask0_3.list --weights weights0_3.list
...
RemoveStreak --merge --inputs input.list --images image4_7.list \
             --masks mask4_7.list --weights weights4_7.list

# Final merge stage
RemoveStreak --merge --inputs input.list --images image0_7.list \
             --masks mask0_7.list --weights weights0_7.list

# Remove stage (executed for each of the 60 raw chip image files, the
  60 detrend chip image files, and each of the 60 - 120 warp and
  difference skycell files).

(XXX NEED to clarify this)

RemoveStreak --remove --inputs input.list --image image.fits \
             --mask mask.fits --weight weight.fits  \
             --outroot output_root [--astrom astrom.fits]

All procecessing output should be stored in the log files in the
<archive> directory.
