Changeset 20780 for trunk/magic/remove/src

Timestamp:

Nov 17, 2008, 2:44:52 PM (17 years ago)

Author:

bills

Message:

if no skycells files are provided, excise the entire image unless operating
at diff stage.
Don't force creation of recovery image unless stage is raw and -replace
was specified. Require recovery file if we're being asked to overwrite
a raw file.

File:

• trunk/magic/remove/src/streaksremove.c (modified) (29 diffs)

trunk/magic/remove/src/streaksremove.c

@@ -290 +290 @@
     char *inputBasename = basename(sf->inImage->name);
     sf->outImage = sFileOpen(config, stage, "OUTPUT", inputBasename, true);
-    sf->recImage = sFileOpen(config, stage, "RECOVERY", inputBasename, true);
+    sf->recImage = sFileOpen(config, stage, "RECOVERY", inputBasename, false);
 
     sf->inMask = sFileOpen(config, stage, "INPUT.MASK", NULL, false);
@@ -490 +490 @@
     }
     
-    
-    if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist")) { ;
-        if (CHIP_LEVEL_INPUT(stage)) {
-            psError(PS_ERR_UNKNOWN, true, "-skycelllist is required for raw and chip stages\n");
-            return NULL;
-        }
-    }
+    if ((argnum = psArgumentGet(argc, argv, "-keepnonwarped"))) {
+        psArgumentRemove(argnum, &argc, argv);
+        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "KEEP_NON_WARPED", 0,
+            "skip excising of non warped pixels", true);
+    }
+        
+    // if skycells are not provided then we have to execise all pixels  unless -keepnonwarped
+    pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist");
 
     if (!pmConfigFileSetsMD(config->arguments, &argc, argv, "ASTROM", "-astrom", NULL)) { ;
@@ -504 +505 @@
         }
     }
+
     if ((argnum = psArgumentGet(argc, argv, "-mask"))) {
         psArgumentRemove(argnum, &argc, argv);
@@ -516 +518 @@
         psArgumentRemove(argnum, &argc, argv);
     }
+
     if ((argnum = psArgumentGet(argc, argv, "-weight"))) {
         psArgumentRemove(argnum, &argc, argv);
@@ -528 +531 @@
         psArgumentRemove(argnum, &argc, argv);
     }
+
     if ((argnum = psArgumentGet(argc, argv, "-tmproot"))) {
         psArgumentRemove(argnum, &argc, argv);
@@ -544 +548 @@
         return NULL;
     }
+
     if ((argnum = psArgumentGet(argc, argv, "-recovery"))) {
         psArgumentRemove(argnum, &argc, argv);
@@ -550 +555 @@
             dir);
         psArgumentRemove(argnum, &argc, argv);
-    } else if (stage == IPP_STAGE_RAW) {
-        psError(PS_ERR_UNKNOWN, true, "-recovery is required for -stage raw\n");
+    } else if ((stage == IPP_STAGE_RAW) && gotReplace) {
+        psError(PS_ERR_UNKNOWN, true, "-recovery is required for -stage raw with -replace\n");
         return NULL;
     }
+
     if ((argnum = psArgumentGet(argc, argv, "-remove"))) {
         if (!gotReplace) {
@@ -589 +595 @@
     }
 
-    
     // TODO: add keyword indicating that streaks have been removed
     if (!psFitsWriteBlank(sfiles->outImage->fits, imageHeader, NULL)) {
@@ -597 +602 @@
     }
     // TODO: add keyword indicating that this is the recovery image
-    if (!psFitsWriteBlank(sfiles->recImage->fits, imageHeader, NULL)) {
+    if (sfiles->recImage && !psFitsWriteBlank(sfiles->recImage->fits, imageHeader, NULL)) {
         psError(PS_ERR_IO, false, "failed to write primary header to %s",
             sfiles->recImage->resolved_name);
@@ -619 +624 @@
         }
         // TODO: add keyword indicating that this is the recovery image
-        if (!psFitsWriteBlank(sfiles->recMask->fits, maskHeader, NULL)) {
+        if (sfiles->recMask && !psFitsWriteBlank(sfiles->recMask->fits, maskHeader, NULL)) {
             psError(PS_ERR_IO, false, "failed to write primary header to %s",
                 sfiles->recMask->resolved_name);
@@ -641 +646 @@
         }
         // TODO: add keyword indicating that this is a recovery image
-        if (!psFitsWriteBlank(sfiles->recWeight->fits, weightHeader, NULL)) {
+        if (sfiles->recWeight && !psFitsWriteBlank(sfiles->recWeight->fits, weightHeader, NULL)) {
             psError(PS_ERR_IO, false, "failed to write primary header to %s",
                 sfiles->recWeight->resolved_name);
@@ -764 +769 @@
             streaksremoveExit("", PS_EXIT_DATA_ERROR);
         }
-        in->numRows = in->image->numRows;
     }  else {
         if (stage != IPP_STAGE_RAW) {
@@ -784 +788 @@
 setFitsOptions(sFile *sfile, int bitpix, float bscale, float bzero)
 {
+    if (!sfile) {
+        return;
+    }
     if (sfile->fits->options) {
         psFree(sfile->fits->options);
@@ -822 +829 @@
     }
 
-//    setFitsOptions(out, bitpix, bscale, bzero);
-//    setFitsOptions(rec, bitpix, bscale, bzero);
+#ifdef COMPRESS_OUTPUT
+    setFitsOptions(out, bitpix, bscale, bzero);
+    setFitsOptions(rec, bitpix, bscale, bzero);
+#endif
 }
 
@@ -859 +868 @@
 
 static void
-createRecoveryImage(sFile *rec, sFile *in, int extnum)
-{
-    rec->image = psImageRecycle(rec->image, in->image->numCols, in->image->numRows, in->image->type.type);
-    if (!rec->image) {
-        psError(PS_ERR_UNKNOWN, false, "failed to allocate recoveyr image for extnsion %d", extnum);
+createNewImage(sFile *out, sFile *in, int extnum, double initValue)
+{
+    out->image = psImageRecycle(out->image, in->image->numCols, in->image->numRows, in->image->type.type);
+    if (!out->image) {
+        psError(PS_ERR_UNKNOWN, false, "failed to allocate image for extnsion %d", extnum);
         streaksremoveExit("", PS_EXIT_UNKNOWN_ERROR);
     }
-    psImageInit(rec->image, 0);
+    psImageInit(out->image, initValue);
+}
+
+static void
+setupImageRefs(sFile *out, sFile *recovery, sFile *in, int extnum, bool exciseAll)
+{
+    if (exciseAll) {
+        // set the output to an image all NANs
+        createNewImage(out, in, extnum, NAN);
+        if (recovery) {
+            // save the whole input as the recovery image
+            recovery->image = (psImage*) psMemIncrRefCounter(in->image);
+        }
+    } else {
+        // output is the input image
+        out->image = (psImage*) psMemIncrRefCounter(in->image);
+        if (recovery) {
+            // create a recovery image initialized to zero
+            createNewImage(recovery, in, extnum, 0);
+        }
+    }
 }
 
 // set the image for the output files to the contents of the corresponding input file.
 static bool
-readAndCopyToOutput(streakFiles *sf)
+readAndCopyToOutput(streakFiles *sf, bool exciseAll)
 {
     bool    updateAstrometry = false;
@@ -893 +922 @@
     if (sf->stage == IPP_STAGE_CHIP) {
         // For the chip level files, copy the WCS from the astrometry file to the header
-        updateAstrometry = true;
+        // updateAstrometry = true;
         if (!pmAstromWriteWCS(sf->inImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
             psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
@@ -900 +929 @@
     }
     sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header);
-    sf->recImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header);
+    if (sf->recImage) {
+        sf->recImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header);
+    }
 
     if (!SFILE_IS_IMAGE(sf->inImage)) {
@@ -910 +941 @@
 
     if (sf->inImage->image) {
-        sf->outImage->image = (psImage*) psMemIncrRefCounter(sf->inImage->image);
-        createRecoveryImage(sf->recImage, sf->inImage, sf->extnum);
+        setupImageRefs(sf->outImage, sf->recImage, sf->inImage, sf->extnum, exciseAll);
     }  else {
-        // we don't do a reference if we have an imagecube
+        // Image cubes are handeled specially
     }
 
     // set up the compression parameters
-    //  copyFitsOptions(sf->outImage, sf->recImage, sf->inImage);
+#ifdef COMPRESS_OUTPUT
+    // compression of the image pixels is disabled for now. Some consortium members
+    // have problems reading them
+    copyFitsOptions(sf->outImage, sf->recImage, sf->inImage);
 
     // XXX: TODO: can we derive these values from the input header?
     // psFitsCompressionGet(sf->inImage->image) gives compression none
     // perhaps we should just use the definition of COMP_IMG in the configuration 
-    //    psFitsSetCompression(sf->outImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-
-    psFitsSetCompression(sf->recImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
+    psFitsSetCompression(sf->outImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
+#endif
+    if (sf->recImage) {
+        psFitsSetCompression(sf->recImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
+    }
 
     if (sf->inMask) {
         readImage(sf->inMask, sf->extnum, sf->stage);
-        sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
         if (updateAstrometry) {
             pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
         }
-        sf->outMask->image = (psImage*) psMemIncrRefCounter(sf->inMask->image);
-        sf->recMask->image = (psImage*) psMemIncrRefCounter(sf->inMask->image);
-
-        createRecoveryImage(sf->recMask, sf->inMask, sf->extnum);
+        sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
+        if (sf->recMask) {
+            sf->recMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
+        }
+
+        setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll);
+
+        // XXX: see note above
         copyFitsOptions(sf->outMask, sf->recMask, sf->inMask);
-
-        // XXX: see note above
         psFitsSetCompression(sf->outMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
         psFitsSetCompression(sf->recMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
@@ -945 +981 @@
     if (sf->inWeight) {
         readImage(sf->inWeight, sf->extnum, sf->stage);
+        if (updateAstrometry) {
+            pmAstromWriteWCS(sf->inWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
+        }
         sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);
-        if (updateAstrometry) {
-            pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
-        }
-        sf->outWeight->image = (psImage*) psMemIncrRefCounter(sf->inWeight->image);
-
-        sf->recWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);
-        createRecoveryImage(sf->recWeight, sf->inWeight, sf->extnum);
+        if (sf->recWeight) {
+            sf->recWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);
+        }
+        setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll);
 
         copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight);
@@ -970 +1006 @@
 writeImage(sFile *sfile, psString extname, int extnum)
 {
-    if (!sfile->image) {
+    if (!sfile || !sfile->image) {
         return;
     }
@@ -987 +1023 @@
 writeImageCube(sFile *sfile, psArray *imagecube, psString extname, int extnum)
 {
+    if (!imagecube) {
+        return;
+    }
     if (!psFitsWriteImageCube(sfile->fits, sfile->header, imagecube, extname)) {
         psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d", 
@@ -997 +1036 @@
 
 static void
-writeImages(streakFiles *sf, bool exciseAll)
+writeImages(streakFiles *sf, bool exciseImageCube)
 {
     psString extname = NULL;
@@ -1005 +1044 @@
     }
     if (sf->inImage->numZPlanes == 0)  {
+        // note exciseing complete images is handled in readAndCopyToOutput
         writeImage(sf->outImage, extname, sf->extnum);
         writeImage(sf->recImage, extname, sf->extnum);
     } else {
+        // we have an image cube
         double initValue;
-        if (exciseAll) {
+        if (exciseImageCube) {
+            // copy the entire input image to the recovery image
             writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum);
             initValue = NAN;
         } else {
+            // otherwise write it to the output 
             writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
             initValue = 0;
         }
 
-        // borrow one of the images from the imagecube and set it to NAN
+        // borrow one of the images from the imagecube and set it to init value
         psImage *image = psArrayGet (sf->inImage->imagecube, 0);
         psMemIncrRefCounter(image);
         psImageInit(image, initValue);
-        if (exciseAll) {
+        if (exciseImageCube) {
             sf->outImage->image = image;
             writeImage(sf->outImage, extname, sf->extnum);
         } else {
-            sf->recImage->image = image;
-            writeImage(sf->recImage, extname, sf->extnum);
-        }
-        // Assumption there are no mask and weight images for video cells
+            // write zero valued image to reccovery
+            if (sf->recImage) {
+                sf->recImage->image = image;
+                writeImage(sf->recImage, extname, sf->extnum);
+            }
+        }
+        // Assumption: there are no mask and weight images for video cells
         return;
     }
@@ -1044 +1090 @@
 closeImage(sFile *sfile)
 {
-    if (!sfile->fits) {
+    if (!sfile) {
         return;
     }
-    if (!psFitsClose(sfile->fits)) {
+    if (sfile->fits && !psFitsClose(sfile->fits)) {
         psError(PS_ERR_IO, false, "failed to close image to %s", 
             sfile->resolved_name);
@@ -1253 +1299 @@
 
     double imageValue  = psImageGet (sfiles->inImage->image,  x, y);
-    psImageSet (sfiles->recImage->image,  x, y, imageValue);
+    if (sfiles->recImage) {
+        psImageSet (sfiles->recImage->image,  x, y, imageValue);
+    }
 
 #define ACTUALLY_MASK_PIXEL
@@ -1407 +1455 @@
 }
 
+static StreakPixels * 
+getStreakPixels(streakFiles *sfiles, Streaks *streaks)
+{
+    StreakPixels *pixels;
+#ifdef notyet
+    Streaks *componentStreaks = NULL;
+    Streaks *streaksToProcess;
+
+    bool only_component_streaks = false; // XXX flesh this out and make it an option
+    if (only_component_streaks) {
+        componentStreaks = restrictStreaksToComponent(sfiles, streaks);
+        streaksToProcess = componentStreaks;
+    } else { 
+        streaksToProcess = streaks;
+    }
+
+    pixels = streak_on_component (streaksToProcess, sfiles->astrom,
+                                  sfiles->inImage->numCols, sfiles->inImage->numRows);
+    if (componentStreaks) {
+        psFree(componentStreaks);
+    }
+#else
+        pixels = streak_on_component (streaks, sfiles->astrom,
+                                      sfiles->inImage->numCols, sfiles->inImage->numRows);
+#endif
+    return pixels;
+}
 int
 main(int argc, char *argv[])
@@ -1433 +1508 @@
     streakFiles *sfiles = openFiles(config);
 
+    bool exciseAll = false;
+    // --keepnonwarped is a test and debug mode
     bool status;
-    bool checkWarpedPixels = ! psMetadataLookupBool(&status, config->arguments, "INCLUDE_NON_WARPED");
-    if (checkWarpedPixels) {
+    bool keepNonWarpedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_WARPED");
+
+    // we need to check for non warped pixels unless we've been asked not to or the stage is diff
+    // (By definition pixels in diff images were included in a diff)
+    bool checkNonWarpedPixels = ! (keepNonWarpedPixels || (sfiles->stage == IPP_STAGE_DIFF) );
+
+    if (checkNonWarpedPixels ) {
         // From ICD:
         // In the raw and detrended images, the pixels which were not
@@ -1441 +1523 @@
         // 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.
-
-        computeWarpedPixels(sfiles);
+        // if no skycells are provided sfiles->exciseAll is set to true
+
+        if (! computeWarpedPixels(sfiles) ) {
+            // we have no choice to excise all pixels
+            exciseAll = true;
+        }
     }
     
@@ -1457 +1543 @@
     }
 
-    // Iterate through components of image
+    // Iterate through components of input files
     do {
-        bool exciseAll = false;
-
-        // read the images and copy the pixels image from the inputs to the outputs
-        // This also sets up the astrometry and initializes the recovery images
-        if (!readAndCopyToOutput(sfiles)) {
-            // this extension is not an image skip (video descriptor?) 
-            // XXX is there anything else that could be in an input file that we need to handle?
+        bool exciseImageCube = false;
+
+        // read the images and copy the data from the inputs to the outputs
+        // This also sets up the astrometry and initializes the recovery image(s)
+
+        if (!readAndCopyToOutput(sfiles, exciseAll)) {
+            // false return value indicates that this extension is not an image and
+            // the contents has already been copied to the output file.
+            // we don't need to process this extesion for streaks.
             continue;
         }
 
-        // Identify pixels to mask because of streaks
-
-#ifdef notyet
-        Streaks *componentStreaks = NULL;
-        Streaks *streaksToProcess;
-
-        bool only_component_streaks = false; // XXX flesh this out and make it an option
-        if (only_component_streaks) {
-            componentStreaks = restrictStreaksToComponent(sfiles, streaks);
-            streaksToProcess = componentStreaks;
-        } else { 
-            streaksToProcess = streaks;
-        }
-
-        pixels = streak_on_component (streaksToProcess, sfiles->astrom,
-                                      sfiles->inImage->numCols, sfiles->inImage->numRows);
-        if (componentStreaks) {
-            psFree(componentStreaks);
-        }
-#else
-        pixels = streak_on_component (streaks, sfiles->astrom,
-                                      sfiles->inImage->numCols, sfiles->inImage->numRows);
-#endif
-        // XXX: 
-        // 
-        // PFS: Need to get mask weight for PM_MASK_DETECTOR.  Is this stored
-        // in a config somewhere?  Not sure how to properly set the maskValue.
-        // WES just pick some numbers
-
-        if (checkWarpedPixels) {
-            exciseNonWarpedPixels(sfiles, MASK_BAD_WARP);
-        }
-
-        if (sfiles->inImage->image) {
-            for (i = 0; i < psArrayLength (pixels); ++i) {
-                pixelPos = psArrayGet (pixels, i);
-
-                if (!checkWarpedPixels || warpedPixel(sfiles, pixelPos)) {
-
-                    excisePixel(sfiles, pixelPos->x, pixelPos->y, MASK_STREAK);
-
-                } else {
-                    // This pixel was not included in any warp and has thus already excised above
+        if (!exciseAll) {
+            // Identify pixels to mask because of streaks
+
+            pixels = getStreakPixels(sfiles, streaks);
+
+            // XXX: 
+            // 
+            // PFS: Need to get mask weight for PM_MASK_DETECTOR.  Is this stored
+            // in a config somewhere?  Not sure how to properly set the maskValue.
+
+
+            if (checkNonWarpedPixels) {
+                exciseNonWarpedPixels(sfiles, MASK_BAD_WARP);
+            }
+            
+            if (sfiles->inImage->image) {
+                for (i = 0; i < psArrayLength (pixels); ++i) {
+                    pixelPos = psArrayGet (pixels, i);
+
+                    if (!checkNonWarpedPixels || warpedPixel(sfiles, pixelPos)) {
+
+                        excisePixel(sfiles, pixelPos->x, pixelPos->y, MASK_STREAK);
+
+                    } else {
+                        // This pixel was not included in any warp and has thus already excised above
+                    }
                 }
+            } else { 
+                // this component contains an image cube, excise it completely
+                exciseImageCube = true;
             }
-        } else { 
-            // This video cell is touched by a streak excise
-            // For now we excise all video cells.
-            // if (psArrayLength(pixels))
-            exciseAll = true;
-        }
-        psArrayElementsFree (pixels);
+            psArrayElementsFree (pixels);
+        }
 
         // write out the destreaked temporary images and the recovery images
-        writeImages(sfiles, exciseAll);
+        writeImages(sfiles, exciseImageCube);
 
     } while (streakFilesNextExtension(sfiles));