Changeset 24557 for branches/eam_branches/20090522/magic/remove/src/streaksremove.c

Timestamp:

Jun 25, 2009, 2:00:56 PM (17 years ago)

Author:

eugene

Message:

merging changes from head

Location:

branches/eam_branches/20090522

Files:

• . (modified) (1 prop)
• magic/remove/src (modified) (1 prop)
• magic/remove/src/streaksremove.c (modified) (19 diffs)

branches/eam_branches/20090522/magic/remove/src

@@ -3 +3 @@
 streakscompare
 streaksrelease
+makefile

@@ -3 +3 @@
 streakscompare
 streaksrelease
+makefile

branches/eam_branches/20090522/magic/remove/src/streaksremove.c

@@ +1 @@
+/*
+ * streaksremove
+ *
+ * Convert satellite streak detctions into masks and remove the covered pixels from the
+ * input images.
+ * Optionally swap the inputs with the outputs so that subsequent references to the original
+ * images access the destreaked versions.
+ */
+
 #include "streaksremove.h"
 
@@ -23 +32 @@
     }
 
-    psU32 maskStreak = 0;
-    psU32 maskMask = 0;
+    // Values to set for masked pixels
+    psU32 maskStreak = 0;           // for the image and weight (usually NAN, MAXINT for integer images)
+    psU32 maskMask = 0;             // value looked up for MASK.STREAK 
 
     psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS");
 
+    // call Paul Sydney's code to parse the streaks file that DetectStreaks produced
     Streaks *streaks = readStreaksFile(streaksFileName);
     if (!streaks) {
-        psErrorStackPrint(stderr, "failed to read streaks file: %s", streaksFileName);
+        psError(PS_ERR_UNKNOWN, "failed to read streaks file: %s", streaksFileName);
         streaksExit("", PS_EXIT_PROG_ERROR);
     }
 
+    // open all of the input and output files, save their descriptions in the streakFiles struct
     streakFiles *sfiles = openFiles(config, true, argv[0]);
     setupAstrometry(sfiles);
 
+    // Optionally we can set pixels that are masked to NAN since they couldn't have been
+    // examined for streaks. Usually this is done by the distribution system just prior
+    // to release
     bool nanForRelease = psMetadataLookupBool(&status, config->arguments, "NAN_FOR_RELEASE");
     if (nanForRelease && (sfiles->inMask == NULL)) {
@@ -52 +67 @@
 
     if (checkNonWarpedPixels ) {
-        // From ICD:
+        // From magic ICD:
         // In the raw and detrended images, the pixels which were not
         // included in any of the streak-processed warps must also be masked.
@@ -69 +84 @@
     
     if (sfiles->stage == IPP_STAGE_RAW) {
-        // copy PHU to output files
+        // Except for raw stage, all of our (GPC1) files have one image extension.
+        // Raw files have a phu and multiple extensions, one per chip
+        // Since this is a raw file, copy it's PHU to output files
         copyPHU(sfiles, true);
 
@@ -82 +99 @@
     int totalStreakPixels = 0;
 
-    // Iterate through each component of the input (there is only one except for raw images)
+    // Iterate through each component of the input (except for raw images there is only one)
     do {
         bool exciseImageCube = false;
@@ -110 +127 @@
             psTimerStart("GET_STREAK_PIXELS");
 
-            StreakPixels *pixels = streak_on_component (streaks, sfiles->astrom,
-                                      sfiles->inImage->numCols, sfiles->inImage->numRows);
-
+            // call Paul Sydney's code to compute the set of pixels that are covered by the detected streaks
+            StreakPixels *pixels = streak_on_component(streaks, sfiles->astrom, sfiles->inImage->numCols,
+                                                        sfiles->inImage->numRows);
             psLogMsg("streaksremove", PS_LOG_INFO, "time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS"));
-
             
+            // if this extension contained an image, excise the streaked pixels.
+            // otherwise it contained an image cube (video cell) which is handled in the if block
             if (sfiles->inImage->image) {
                 if (checkNonWarpedPixels) {
@@ -126 +144 @@
                     psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED"));
                 }
+
                 totalStreakPixels +=  psArrayLength(pixels);
+
                 psTimerStart("REMOVE_STREAKS");
+
+                // for each pixel covered by the streak
                 for (int i = 0; i < psArrayLength (pixels); ++i) {
                     PixelPos *pixelPos = psArrayGet (pixels, i);
 
+                    // if this pixel was not part of the warp, skip because it has already been
+                    // excised (unless we weren't asked to check)
                     if (!checkNonWarpedPixels || warpedPixel(sfiles, pixelPos)) {
 
@@ -140 +164 @@
                     }
                 }
+
                 psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS"));
 
                 if (nanForRelease) {
+                    // set any pixels that were masked, to NAN (unless they are already NAN)
                     setMaskedToNAN(sfiles, maskMask, true);
                 }
 
             } else { 
-                // this component contains an image cube, excise it completely
+                // this component contains an image cube
+                // For now excise it completely
                 exciseImageCube = true;
             }
@@ -155 +182 @@
 
         if (sfiles->stage == IPP_STAGE_CHIP) {
+            // as a convience to the user of the output, replace the bogus WCS transform in the
+            // chip processed files with the data calcuated by psastro at the camera stage
+            // (actually we use a linear approximation)
             updateAstrometry(sfiles);
         }
 
-        // write out the destreaked temporary images and the recovery images
+        // write the destreaked "temporary" images and the recovery images
         writeImages(sfiles, exciseImageCube);
 
@@ -168 +198 @@
     psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, 100. * totalStreakPixels / totalPixels);
 
+    // all done close the files. This is where the files are written so it can take a long time.
+
     psTimerStart("CLOSE_IMAGES");
-    // close all files
+
     closeImages(sfiles);
+
     psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", psTimerClear("CLOSE_IMAGES"));
+
+
+    if (!replicateOutputs(sfiles)) {
+        psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
+        deleteTemps(sfiles);
+        psErrorStackPrint(stderr, "");
+        exit(PS_EXIT_UNKNOWN_ERROR);
+    }
 
     // NOTE: from here on we can't just quit if something goes wrong.
     // especially if we're working at the raw stage
-
-    if (!replicateOutputs(sfiles)) {
-        psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
-        psErrorStackPrint(stderr, "");
-        exit(PS_EXIT_UNKNOWN_ERROR);
-    }
+    // turn off automatic deletion of output files by streaksExit
+    setStreakFiles(NULL);
 
     if (psMetadataLookupBool(&status, config->arguments, "REPLACE")) {
         //     swap the instances for the input and output
-        //     Note this is a database operation. No file I/O is performed
+        //     Note this is a nebulous database operation. No file I/O is performed
         if (!swapOutputsToInputs(sfiles)) {
-            psError(PS_ERR_UNKNOWN, false, "failed to swap files");
-
-            // XXX: Now what? I guess swapOutputsToInputs will need to undo anything that
-            // it has done and give a detailed report of what happened
-
-            psErrorStackPrint(stderr, "");
+            // XXX: Now what?
+            // It is up to the program that reverts failed destreak runs to insure that
+            // any input files that have been swapped are restored
+
+            psErrorStackPrint(stderr, "failed to swap files");
+
+            // XXX: pick a specific error code for this failure
             exit(PS_EXIT_UNKNOWN_ERROR);
         }
-
-        if (psMetadataLookupBool(&status, config->arguments, "REMOVE")) {
-            //      delete the temporary storage objects (which now points to the original image(s)
-            if (!deleteTemps(sfiles)) {
-                psError(PS_ERR_UNKNOWN, false, "failed to delete temporary files");
-                // XXX: Now what? At this point the output files have been swapped, so we can't
-                // repeat the operation.
-
-                // Returning error status here is problematic. The inputs have been streak removed
-                // but they're still lying around
-                // Maybe just print an error message and
-                // let other system tools clean up
-                psErrorStackPrint(stderr, "");
-                exit(PS_EXIT_UNKNOWN_ERROR);
-            }
-        }
-    }
+    }
+    // all done. Clean up to look for memory leaks.
 
     psFree(sfiles);
@@ -403 +426 @@
             usage();
         }
-        psString dir = pathToDirectory(argv[argnum]);
-        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "directory for temporary files",
-            dir);
-        psFree(dir);
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "path for (temporary if replae) output files",
+            argv[argnum]);
         psArgumentRemove(argnum, &argc, argv);
     } else {
@@ -415 +436 @@
     if ((argnum = psArgumentGet(argc, argv, "-recovery"))) {
         psArgumentRemove(argnum, &argc, argv);
-        psString dir = pathToDirectory(argv[argnum]);
-        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "RECOVERY", 0, "directory for recovery files",
-            dir);
-        psFree(dir);
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "RECOVERY", 0, "path for recovery files",
+            argv[argnum]);
         psArgumentRemove(argnum, &argc, argv);
     } else if ((stage == IPP_STAGE_RAW) && gotReplace) {
         psError(PS_ERR_UNKNOWN, true, "-recovery is required for -stage raw with -replace\n");
         usage();
-    }
-
-    if ((argnum = psArgumentGet(argc, argv, "-remove"))) {
-        if (!gotReplace) {
-            psError(PS_ERR_UNKNOWN, true, "-replace is required with -remove\n");
-            usage();
-        }
-        psArgumentRemove(argnum, &argc, argv);
-        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "REMOVE", 0, "remove original files",
-            true);
     }
 
@@ -451 +460 @@
 updateAstrometry(streakFiles *sf)
 {
+    // XXX: why do I check this here? Shouldn't it be just around the call to linearizeTransforms?
     if (sf->bilevelAstrometry) {
 
@@ -501 +511 @@
         }
     }
-    sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header);
+    sf->outImage->header =  psMemIncrRefCounter(sf->inImage->header);
     if (sf->recImage) {
-        sf->recImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header);
-    }
+        sf->recImage->header = psMetadataCopy(NULL, sf->inImage->header);
+        addRecoveryKeyword(sf->recImage->header);
+    }
+    addDestreakKeyword(sf->outImage->header);
 
     if (!SFILE_IS_IMAGE(sf->inImage)) {
@@ -520 +532 @@
 
     // set up the compression parameters
-#ifdef STREAKS_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);
-    if (sf->recImage) {
-        psFitsSetCompression(sf->recImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-    }
-#endif
+    copyFitsOptions(sf->outImage, sf->recImage, sf->inImage, sf->tiles);
 
     if (sf->inMask) {
@@ -539 +539 @@
             sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
             if (sf->recMask) {
-                sf->recMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
-            }
+                sf->recMask->header = psMetadataCopy(NULL, sf->outMask->header);
+                addRecoveryKeyword(sf->recMask->header);
+            }
+            addDestreakKeyword(sf->outMask->header);
             if (updateAstrometry) {
                 pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
@@ -554 +556 @@
             }
 
-#ifdef STREAKS_COMPRESS_OUTPUT
-            // XXX: see note above
-            copyFitsOptions(sf->outMask, sf->recMask, sf->inMask);
-            psFitsSetCompression(sf->outMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-            if (sf->recMask) {
-                psFitsSetCompression(sf->recMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-            }
+            copyFitsOptions(sf->outMask, sf->recMask, sf->inMask, sf->tiles);
             if (sf->outChMask) {
-                copyFitsOptions(sf->outChMask, sf->recChMask, sf->inMask);
-                psFitsSetCompression(sf->outChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-                if (sf->recChMask) {
-                    psFitsSetCompression(sf->recChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-                }
-            }
-#endif
+                copyFitsOptions(sf->outChMask, sf->recChMask, sf->inMask, sf->tiles);
+            }
         }
     }
@@ -576 +567 @@
         sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);
         if (sf->recWeight) {
-            sf->recWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);
-        }
+            sf->recWeight->header = psMetadataCopy(NULL, sf->outWeight->header);
+            addRecoveryKeyword(sf->recWeight->header);
+        }
+        addDestreakKeyword(sf->outWeight->header);
         if (updateAstrometry) {
             pmAstromWriteWCS(sf->inWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
@@ -583 +576 @@
         setupImageRefs(sf->outWeight, sf->recWeight, sf->inWeight, sf->extnum, exciseAll);
 
-#ifdef STREAKS_COMPRESS_OUTPUT
-        copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight);
-        // XXX: see note above
-        psFitsSetCompression(sf->outWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-        if (sf->recWeight) {
-            psFitsSetCompression(sf->recWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-        }
-#endif
-    }
-    // and for raw images, create sub images that represent the actual image
-    // area (no overscan)
+        copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight, sf->tiles);
+    }
 
     return true;
 }
-
-
 
 static void

@@ -3 +3 @@
 streakscompare
 streaksrelease
+makefile