Index: trunk/doc/design/ipp-magic-icd.txt =================================================================== --- trunk/doc/design/ipp-magic-icd.txt (revision 19842) +++ trunk/doc/design/ipp-magic-icd.txt (revision 19842) @@ -0,0 +1,213 @@ + +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: + +- .dif..fits - difference image +- .dif..mk.fits - binary mask image +- .dif..wt.fits - variance image + +A common base filename to use is .dif.. The command line +interface allows arbitrary file paths and names to be specified in the +--detect stage using: + +-i or --image +-m <8bit mask image> or --mask <8bit mask image> +-w or --weight +--outroot + +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 . (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. + +- .clusters is a binary file containing all pixel clusters detected +- .streaks is a text file containing all streaks detected at the + current stage of processing +- _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 .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 +-I or --images +-M or --masks +-W or --weights +--wcsList + +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 +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 + directory. Index: trunk/doc/design/remove-streaks-remove.c =================================================================== --- trunk/doc/design/remove-streaks-remove.c (revision 19842) +++ trunk/doc/design/remove-streaks-remove.c (revision 19842) @@ -0,0 +1,100 @@ +streaks_remove(inImage, // Input image + inMask, // Input mask (optional) + inWeight, // Input weight map (optional) + inAstrom, // Input astrometry (optional) + maskValue, // Mask value to set + recoveryRoot, // Root name for recovery files + streaks, // Streaks table + warps // Warps table + ) +{ + // Lock Nebulous files to prevent access while changing in-place + neb_lock(inImage); + if (inMask) neb_lock(inMask); + if (inWeight) neb_lock(inWeight); + + // Optionally load astrometry + if (inAstrom) astrom = astrometry_read(inAstrom); + + // Iterate through components of image + foreach extension(inImage) { + // Read header and extract astrometry + header = header_read(inImage); + if (!inAstrom) astrom = astrometry_from_header(header); + + // Identify pixels to mask because of streaks + streak_pixels = pixel_list_initialise(NONE); // List of pixels to be excised because of streak + foreach streak (streaks) { + // Calculate streak parameters on image + image_streak = streak_on_component(streak, astrom, num_cols, num_rows); + // Only bother with streaks that touch the image + next if image_streak not defined; + + foreach pixel specified by local_streak { + add_pixel(streak_pixels, pixel); + } + } + + // Identify pixels to mask because not in a warp + warp_pixels = pixel_list_initialise(ALL); // List of pixels to be excised because not in a warp + foreach warp (warps) { + // Calculate warp region on image + image_warp = warp_on_component(warp, astrom, num_cols, num_rows); + foreach pixel specified by image_warp { + remove_pixel(warp_pixels, pixel); + } + } + + // Mark the header as "destreaked" + mark_header(inImage); + if (inMask) mark_header(inImage); + if (inWeight) mark_header(inImage); + + if (size(streak_pixels) > 0 || size(warp_pixels) > 0) { + // Read the pixels + if (num_streaks == 0) { + image = image_read(inImage); + if (inMask) mask = image_read(inMask); + if (inWeight) weight = image_read(inWeight); + + recovery_image = initialise_image(0); + if (inMask) recovery_mask = initialise_image(0); + if (inWeight) recovery_weight = initialise_image(0); + } + + // Mask the pixels + foreach masked_pixel (streak_pixels, warp_pixels) { + recovery_image[x,y] = image[x,y]; + image[x,y] = NAN; + if (inMask) { + recovery_mask[x,y] = mask[x,y]; + mask[x,y] = maskValue; + } + if (inWeight) { + recovery_weight[x,y] = weight[x,y]; + weight[x,y] = NAN; + } + } + + // Save the outputs + image_save(image); + image_write_compressed(recoveryRoot + ".im.fits", recovery_image); + + if (inMask) { + image_save(mask); + image_write_compressed(recoveryRoot + ".mk.fits", recovery_mask); + } + if (inWeight) { + image_save(mask); + image_write_compressed(recoveryRoot + ".wt.fits", recovery_weight); + } + } + } + + // Unlock Nebulous + neb_unlock(inImage); + if (inMask) neb_unlock(inMask); + if (inWeight) neb_unlock(inWeight); + + return SUCCESS; +} Index: trunk/doc/design/remove-streaks-remove.txt =================================================================== --- trunk/doc/design/remove-streaks-remove.txt (revision 19842) +++ trunk/doc/design/remove-streaks-remove.txt (revision 19842) @@ -0,0 +1,50 @@ + +basic outline for remove-streaks-remove: + +inputs: + inImage, // Input image + inMask, // Input mask (optional) + inWeight, // Input weight map (optional) + inAstrom, // Input astrometry (optional) + maskValue, // Mask value to set (determined from header) + recoveryRoot, // Root name for recovery files + streaks, // Streaks table + warps // Warps table (list of overlapping warps) + +neb_lock images (inImage, inMask, inWeight) + +load astrometry (if supplied) + +select streaks touching inImage + +create 'recovery images' with 0.0 valued pixels + +loop over image components (64 cells for raw images; 1 cell for detrended, warp, and diff images) + + load image, mask, weight + + loop over selected streak + + skip streak if not on component + + determine streak pixels + (this function will require calls to the (ra,dec) <-> (x,y) transformations) + + copy streak pixel values to corresponding pixels in recovery image + + set streak pixel in input image and weight to 0.0, set mask bit + + save image, mask, weight (to temporary file) + + save recovery image, mask weight + +replace input image, mask, weight with temporary files + +neb_unlock images + +--- + +We propose to supply the APIs for everything except the code to +determine the streak pixels, though we will supply APIs to convert ra, +dec to/from x,y + Index: trunk/doc/design/remove-streaks-remove.v2.c =================================================================== --- trunk/doc/design/remove-streaks-remove.v2.c (revision 19842) +++ trunk/doc/design/remove-streaks-remove.v2.c (revision 19842) @@ -0,0 +1,150 @@ +streaks_remove(inImage, // Input image ?: char*? + inMask, // Input mask (optional) + inWeight, // Input weight map (optional) + inAstrom, // Input astrometry (optional) + maskValue, // Mask value to set + recoveryRoot, // Root name for recovery files + streaks, // Streaks table ?: table of RA/Dec endpoints? + warps // Warps table ?: what is this format? + ) +{ + // Lock Nebulous files to prevent access while changing in-place + neb_lock(inImage); + if (inMask) neb_lock(inMask); // ?: Is the call format, neb_lock (const char* path)? + if (inWeight) neb_lock(inWeight); + + // Optionally load astrometry + if (inAstrom) astrom = astrometry_read(inAstrom); + // ?: Is this a UH provided routine? EAM : yes + // ?: What is the structure for astrom? EAM : TBD + // ?: Definition of higher order WCS keywords? : EAM : added to ICD + + // ?: For detrend and raw files, are the higher order terms stored in + // a separate file? EAM : yes + + // Is the unique higher order WCS per extension or is it based on an image. EAM : see ICD + + // Iterate through components of image + foreach extension(inImage) { + // Read header and extract astrometry + header = header_read(inImage); + if (!inAstrom) astrom = astrometry_from_header(header); + // ?: Is this an affine transform only? + + // EAM : no. the optical distortion is third order. However, the non-linear terms are + // perturbations on the linear terms, and assuming there are no caustics in the optics + // (should not be!), there is a one-to-one relationship between input and output + // coordinates. This gives us two simplifications: 1) we can use a linear approximation + // for localized patches (we use this as an optimization in the warping, and I recommed + // we use it here as well). 2) the inverse transformation can be closely represented + // by a polynomial of the same order. Or, to be more accurate, an approximate inverse + // transformation of the same order only makes small errors. in our analysis, we + // define one direction (chip->focal plane->tangent plane->sky) to be the primary + // transformation, and the inverse is understood to be lower accuracy and used only for + // finding relationships. (ie, match reference star to detection). + + // Identify pixels to mask because of streaks + streak_pixels = pixel_list_initialise(NONE); // List of pixels to be excised because of streak + foreach streak (streaks) { + // Calculate streak parameters on image + image_streak = streak_on_component(streak, astrom, num_cols, num_rows); + // Only bother with streaks that touch the image + next if image_streak not defined; + + foreach pixel specified by local_streak { + add_pixel(streak_pixels, pixel); + } + } + + // Identify pixels to mask because not in a warp + warp_pixels = pixel_list_initialise(ALL); // List of pixels to be excised because not in a warp + foreach warp (warps) { + // Calculate warp region on image + image_warp = warp_on_component(warp, astrom, num_cols, num_rows); + foreach pixel specified by image_warp { + remove_pixel(warp_pixels, pixel); + } + } + + // Mark the header as "destreaked" + mark_header(inImage); + if (inMask) mark_header(inMask); + if (inWeight) mark_header(inWeight); + + // ?: mark_header(inMask); // Record masking in main image or mask? + // ?: mark_header(inWeight);// Record masking in weight image or mask? + + // EAM : these operations are mean to describe the addition of metadata to the headers + // to record that magic has been applied (and any interesting stats, like number of + // pixels removed -- NOTE: that means this should come after the following block). we + // would want to record the metadata in each of image, weight, and mask so that they + // are all self-documenting. + + if (size(streak_pixels) > 0 || size(warp_pixels) > 0) { + // Read the pixels + if (num_streaks == 0) { + image = image_read(inImage); + if (inMask) mask = image_read(inMask); + if (inWeight) weight = image_read(inWeight); + + recovery_image = initialise_image(0); + if (inMask) recovery_mask = initialise_image(0); + if (inWeight) recovery_weight = initialise_image(0); + } + + // Mask the pixels + foreach masked_pixel (streak_pixels, warp_pixels) { + recovery_image[x,y] = image[x,y]; + image[x,y] = NAN; + // ?: Really record it as Not A Number (NAN)? + // EAM : that is our preference: we use NAN to mark pixels with invalid data elsewhere + + if (inMask) { + recovery_mask[x,y] = mask[x,y]; + mask[x,y] = maskValue; + } + if (inWeight) { + recovery_weight[x,y] = weight[x,y]; + weight[x,y] = NAN; + } + } + + // ?: No temporary files? + + // EAM : temporary files are probably required, just ignored in this description. + // We are assuming that the process loops over the cells (extensions) in the image + // to minimize memory footprint, so we can't write the result to the same output + // file. + + // Save the outputs (write for current component/image extension) + image_save(image); + image_write_compressed(recoveryRoot + ".im.fits", recovery_image); + + if (inMask) { + image_save(mask); + image_write_compressed(recoveryRoot + ".mk.fits", recovery_mask); + } + if (inWeight) { + image_save(mask); + image_write_compressed(recoveryRoot + ".wt.fits", recovery_weight); + } + + // ?: Add writing of streak_pixels and warp_pixels list for image? + + // EAM : Paul Price has made an alternate proposal for the storage of the 'recovery + // image', which we used in defining this: make that image consist of only the + // masked pixels, with everything else set to 0.0. His point is that such an + // image, if written RICE compressed, would compress extremely well. It would also + // make recovery extremely simple: just read it in the two image and replace into + // the 'magicked image' any non-zero pixels from the recover image. + + } // if (size(streak_pixels) > 0 || size(warp_pixels) > 0) + } // foreach extension(inImage) + + // Unlock Nebulous + neb_unlock(inImage); + if (inMask) neb_unlock(inMask); + if (inWeight) neb_unlock(inWeight); + + return SUCCESS; +}