Changeset 24283 for trunk/magic/remove/src/streaksremove.c

Timestamp:

May 30, 2009, 11:06:34 AM (17 years ago)

Author:

bills

Message:

add keywords to the headers indicating that an output image is
destreaked or a recovery image

File:

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

trunk/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) {
@@ -34 +45 @@
     }
 
+    // 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"));
 
-    // 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)) {
@@ -182 +213 @@
     }
 
+    // NOTE: from here on we can't just quit if something goes wrong.
+    // especially if we're working at the raw stage
+
     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");
@@ -192 +226 @@
 
             psErrorStackPrint(stderr, "");
+
+            // XXX: pick a specific error code for this failure
             exit(PS_EXIT_UNKNOWN_ERROR);
         }
 
+#ifdef notdef
+        // XXX: we've decided to not do the remove step here
+        // Instead we leave the backup images (the original images) in place and remove them
+        // later on
         if (psMetadataLookupBool(&status, config->arguments, "REMOVE")) {
             //      delete the temporary storage objects (which now points to the original image(s)
@@ -210 +250 @@
             }
         }
-    }
+#endif
+
+    }
+    // all done. Clean up to look for memory leaks.
 
     psFree(sfiles);
@@ -425 +468 @@
     }
 
+#ifdef notdef
+    // This parameter is no longer supported
     if ((argnum = psArgumentGet(argc, argv, "-remove"))) {
         if (!gotReplace) {
@@ -434 +479 @@
             true);
     }
+#endif
 
     if (argc != 1) {
@@ -451 +497 @@
 updateAstrometry(streakFiles *sf)
 {
+    // XXX: why do I check this here? Shouldn't it be just around the call to linearizeTransforms?
     if (sf->bilevelAstrometry) {
 
@@ -501 +548 @@
         }
     }
-    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)) {
@@ -539 +588 @@
             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);
@@ -576 +627 @@
         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);
@@ -597 +650 @@
     return true;
 }
-
-
 
 static void