Changeset 27838 for branches/tap_branches/magic/remove/src/streaksremove.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

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

branches/tap_branches/magic

@@ -1 +1 @@
 magic
 ssa-core-cpp
-Makefile
 Makefile.bak

@@ -1 +1 @@
 magic
 ssa-core-cpp
-Makefile
 Makefile.bak

branches/tap_branches/magic/remove/src

@@ -4 +4 @@
 streakscompare
 streaksrelease
-makefile
+Makefile
+Makefile.in
+config.h
+.deps
+streaksVersionDefinitions.h
+config.h.in
+stamp-h1

@@ -4 +4 @@
 streakscompare
 streaksrelease
-makefile
+Makefile
+Makefile.in
+config.h
+.deps
+streaksVersionDefinitions.h
+config.h.in
+stamp-h1

branches/tap_branches/magic/remove/src/streaksremove.c

@@ -12 +12 @@
 static pmConfig *parseArguments(int argc, char **argv);
 static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll);
-static void exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue);
-static bool warpedPixel(streakFiles *sfiles, int x, int y);
+static long exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue);
+static bool diffedPixel(streakFiles *sfiles, int x, int y);
 static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, psImageMaskType newMaskValue);
 static void writeImages(streakFiles *sf, bool exciseImageCube);
 static void updateAstrometry(streakFiles *sfiles);
 static void censorSources(streakFiles *sfiles, psImageMaskType maskStreak);
-static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonWarpedPixels, psU16 maskStreak);
-
+static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonDiffedPixels, psU16 maskStreak);
+
+char *streaksProgram = "streaksremove";
+
+// Note: For clarity the flow of this program is in main().
+// There is not a lot of error checking is done in main.
+// Until the end, where we might be doing Nebulous operations, called functions exit when an error
+// is encountered.
 int
 main(int argc, char *argv[])
@@ -26 +32 @@
 
     psLibInit(NULL);
-    psTimerStart("STREAKSREMOVE");
+    psTimerStart("TOTAL_TIME");
 
     pmConfig *config = parseArguments(argc, argv);
@@ -33 +39 @@
         streaksExit("", PS_EXIT_CONFIG_ERROR);
     }
-
-    // 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");
@@ -50 +52 @@
     streakFiles *sfiles = openFiles(config, true, argv[0]);
     setupAstrometry(sfiles);
+    sfiles->stats = psMetadataAlloc();
 
     // Optionally we can set pixels that are masked to NAN since they couldn't have been
@@ -61 +64 @@
 
     bool exciseAll = false;
-    // --keepnonwarped is a test and debug mode
-    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 ) {
+    // --keepnondiffed is a test and debug mode
+    bool keepNonDiffedPixels = psMetadataLookupBool(&status, config->arguments, "KEEP_NON_DIFFED");
+
+    // we need to check for non diffed pixels unless we've been asked not to or the stage is diff
+    // (By definition pixels in diff images were included in the difference images)
+    bool checkNonDiffedPixels = ! (keepNonDiffedPixels || (sfiles->stage == IPP_STAGE_DIFF) );
+
+    if (checkNonDiffedPixels ) {
         // From magic ICD:
         // In the raw and detrended images, the pixels which were not
@@ -76 +79 @@
         // if no skycells are provided sfiles->exciseAll is set to true
 
-        psTimerStart("COMPUTE_WARPED_PIXELS");
-        if (! computeWarpedPixels(sfiles) ) {
+        psTimerStart("COMPUTE_DIFFED_PIXELS");
+        if (! computeDiffedPixels(sfiles) ) {
             // we have no choice to excise all pixels
             exciseAll = true;
         }
-        psF64 cwp_t = psTimerClear("COMPUTE_WARPED_PIXELS");
-        psLogMsg("streaksremove", PS_LOG_INFO, "time to compute warped pixels: %f\n", cwp_t);
+        psF64 cdp_t = psTimerClear("COMPUTE_DIFFED_PIXELS");
+        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "COMPUTE_UNDIFFED_PIXELS", PS_META_REPLACE, "time to compute non-diffedpixels", cdp_t);
+        psLogMsg("streaksremove", PS_LOG_INFO, "time to compute diffed pixels: %f\n", cdp_t);
     }
 
@@ -100 +104 @@
     long totalPixels = 0;
     long totalStreakPixels = 0;
-
+    long nonDiffedPixels = 0;
+
+    // accumulators for the various timers
+    psF64 gsp_t = 0;
+    psF64 enw_t = 0;
+    psF64 rms_t = 0;
+    psF64 cs_t = 0;
+    psF64 wi_t = 0;
+    psF64 ua_t = 0;
     // Iterate through each component of the input (except for raw images there is only one)
     do {
@@ -118 +130 @@
 
         // now that we've read the input files, lookup the mask values
-        if (maskStreak == 0) {
-            strkGetMaskValues(sfiles, &maskStreak, &maskMask);
-        }
+        strkGetMaskValues(sfiles);
 
         totalPixels += sfiles->inImage->numRows * sfiles->inImage->numCols;
@@ -132 +142 @@
             StreakPixels *pixels = streak_on_component(streaks, sfiles->astrom, sfiles->inImage->numCols,
                                                         sfiles->inImage->numRows);
+            gsp_t +=  psTimerClear("GET_STREAK_PIXELS");
+            psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "GET_STREAK_PIXELS", PS_META_REPLACE, "", gsp_t);
             psLogMsg("streaksremove", PS_LOG_INFO, "time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS"));
 
@@ -137 +149 @@
             // otherwise it contained an image cube (video cell) which is handled in the if block
             if (sfiles->inImage->image) {
-                if (checkNonWarpedPixels) {
-                    psTimerStart("EXCISE_NON_WARPED");
-
-                    // set non-warped pixels and variance to NAN, mask to maskStreak (since the pixel
+                if (checkNonDiffedPixels) {
+                    psTimerStart("EXCISE_NON_DIFFED");
+
+                    // set non-diffed pixels and variance to NAN, mask to maskStreak (since the pixel
                     // is excised as part of the destreaking process)
-                    exciseNonWarpedPixels(sfiles, maskStreak);
-
-                    psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED"));
+                    nonDiffedPixels += exciseNonDiffedPixels(sfiles, sfiles->maskStreak);
+
+                    enw_t +=  psTimerClear("EXCISE_NON_DIFFED");
+                    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "EXCISE_NON_DIFFED", PS_META_REPLACE, "", enw_t);
+                    psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non diffed pixels: %f\n", enw_t);
                 }
 
                 psTimerStart("REMOVE_STREAKS");
 
-                totalStreakPixels += censorPixels(sfiles, pixels, checkNonWarpedPixels, maskStreak);
-
-                psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS"));
+                totalStreakPixels += censorPixels(sfiles, pixels, checkNonDiffedPixels, sfiles->maskStreak);
+
+                rms_t += psTimerClear("REMOVE_STREAKS");
+                psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "REMOVE_STREAKS", PS_META_REPLACE, "", enw_t);
+                psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", rms_t);
 
                 if (nanForRelease) {
                     // set any pixels that were masked, to NAN (unless they are already NAN)
-                    setMaskedToNAN(sfiles, maskMask, true);
+                    setMaskedToNAN(sfiles, sfiles->maskMask, true);
                 }
 
@@ -170 +186 @@
             // chip processed files with the data calcuated by psastro at the camera stage
             // (actually we use a linear approximation)
+            psTimerStart("UPDATE_ASTROMETRY");
             updateAstrometry(sfiles);
-        }
-
-        censorSources(sfiles, maskStreak);
+            ua_t += psTimerClear("UPDATE_ASTROMETRY");
+            psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "UPDATE_ASTROMETRY", PS_META_REPLACE, "", ua_t);
+            psLogMsg("streaksremove", PS_LOG_INFO, "time to update astrometry: %f\n", ua_t);
+        }
+
+        psTimerStart("CENSOR_SOURCES");
+        censorSources(sfiles, sfiles->maskStreak);
+        cs_t += psTimerClear("CENSOR_SOURCES");
+        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "CENSOR_SOURCES", PS_META_REPLACE, "", cs_t);
 
         // write the destreaked "temporary" images and the recovery images
+        psTimerStart("WRITE_IMAGES");
         writeImages(sfiles, exciseImageCube);
-
+        wi_t += psTimerClear("WRITE_IMAGES");
+        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "WRITE_IMAGES", PS_META_REPLACE, "", wi_t);
 
         psLogMsg("streaksremove", PS_LOG_INFO, "time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT"));
+
     } while (streakFilesNextExtension(sfiles));
 
@@ -185 +211 @@
     psFree(streaks);
 
-    psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, 100. * totalStreakPixels / totalPixels);
+    if (exciseAll) {
+        totalStreakPixels = totalPixels;
+    }
+
+    psF64 streakFraction = (double) totalStreakPixels / totalPixels;
+    psLogMsg("streaksremove", PS_LOG_INFO, "   total pixels:  %8ld\n", totalPixels);
+    psLogMsg("streaksremove", PS_LOG_INFO, "  streak pixels:  %8ld %4.2f%%\n", totalStreakPixels, streakFraction * 100);
+    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "STREAK_FRACTION", PS_META_REPLACE, "", streakFraction);
+
+    psF64 nonDiffedFraction = (double) nonDiffedPixels / totalPixels;
+    psLogMsg("streaksremove", PS_LOG_INFO, "nondiffed pixels:  %8ld %4.2f%%\n", nonDiffedPixels, nonDiffedFraction * 100);
+    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "NONDIFFED_FRACTION", PS_META_REPLACE, "", nonDiffedFraction);
 
     // check the weight and mask files for extra extensions that might be in files
@@ -194 +231 @@
 
     psTimerStart("CLOSE_IMAGES");
-
     closeImages(sfiles);
-
-    psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", psTimerClear("CLOSE_IMAGES"));
-
-
+    psF64 ci_t = psTimerClear("CLOSE_IMAGES");
+    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "CLOSE_IMAGES", PS_META_REPLACE, "", ci_t);
+
+    psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", ci_t);
+
+#ifdef DO_REPLICATE
+    psTimerStart("REPLICATE_OUTPUTS");
     if (!replicateOutputs(sfiles)) {
         psErrorStackPrint(stderr, "failed to replicate output files");
@@ -205 +244 @@
         exit(PS_EXIT_UNKNOWN_ERROR);
     }
+    psF64 ro_t = psTimerClear("REPLICATE_OUTPUTS");
+    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "REPLICATE_OUTPUTS", PS_META_REPLACE, "", ro_t);
+#endif
 
     // NOTE: from here on we can't just quit if something goes wrong.
@@ -214 +256 @@
         //     swap the instances for the input and output
         //     Note this is a nebulous database operation. No file I/O is performed
+        psTimerStart("SWAP_INSTANCES");
         if (!swapOutputsToInputs(sfiles)) {
-            // 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 and that the de-streaked
-            // versions are deleted
+            // any original non-destreaked input files that have been swapped are restored and that the de-streaked
+            // versions are deleted.
 
             psErrorStackPrint(stderr, "failed to swap files");
@@ -224 +266 @@
             // XXX: pick a specific error code for this failure
             exit(PS_EXIT_UNKNOWN_ERROR);
+        }
+        psF64 si_t = psTimerClear("SWAP_INSTANCES");
+        psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "SWAP_INSTANCES", PS_META_REPLACE, "", si_t);
+    }
+
+    psF64 total_time = psTimerClear("TOTAL_TIME");
+    psMetadataAddF32(sfiles->stats, PS_LIST_TAIL, "TOTAL_TIME", PS_META_REPLACE, "", total_time);
+    psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", total_time);
+
+    if (sfiles->statsFile) {
+        const char *statsMDC = psMetadataConfigFormat(sfiles->stats);
+        if (!statsMDC || strlen(statsMDC) == 0) {
+            psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n");
+        } else {
+            fprintf(sfiles->statsFile, "%s", statsMDC);
+            psFree(statsMDC);
+            fclose(sfiles->statsFile);
+            sfiles->statsFile = NULL;
+            psFree(sfiles->stats);
+            sfiles->stats = NULL;
         }
     }
@@ -234 +296 @@
     streaksNebulousCleanup();
     pmConfigDone();
-    psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE"));
     psLibFinalize();
 
@@ -243 +304 @@
 
 static long
-censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonWarpedPixels, psU16 maskStreak)
+censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonDiffedPixels, psU16 maskStreak)
 {
     long streakPixels = 0;
@@ -250 +311 @@
         for (int x = 0; x < sfiles->inImage->numCols; x++) {
             if (psImageGet(pixels, x, y)) {
-                ++streakPixels;
-                if (!checkNonWarpedPixels || warpedPixel(sfiles, x, y)) {
+                if (!checkNonDiffedPixels || diffedPixel(sfiles, x, y)) {
+                    ++streakPixels;
 
                     excisePixel(sfiles, x, y, true, maskStreak);
@@ -257 +318 @@
                 } else {
                     // This pixel was not included in any warp and has thus already excised
-                    // by exciseNonWarpedPixels
+                    // by exciseNonDiffedPixels
                 }
             }
@@ -274 +335 @@
     fprintf(stderr, "\t-weight WEIGHT.fits: weight file to de-streak\n");
     fprintf(stderr, "\t-replace: replace the input images with the output\n");
-    fprintf(stderr, "\t-keepnonwarped: do not exise pixels that were not part of difference processing\n");
+    fprintf(stderr, "\t-keepnondiffed: do not exise pixels that were not part of difference processing\n");
     fprintf(stderr, "\t-transparent val: instead of setting excicsed pixel to NAN add val\n");
     fprintf(stderr, "\t-chip_mask MASK.fits: name of mask for chip stage (camera stage mask is passed tih -mask)\n");
@@ -365 +426 @@
     }
 
-    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 ((argnum = psArgumentGet(argc, argv, "-keepnondiffed"))) {
+        psArgumentRemove(argnum, &argc, argv);
+        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "KEEP_NON_DIFFED", 0,
+            "skip excising of non diffed pixels", true);
     }
 
@@ -378 +439 @@
         psArgumentRemove(argnum, &argc, argv);
         double transparentStreaks = atof(argv[argnum]);
-        psMetadataAddF64(config->arguments, PS_LIST_TAIL, "TRANSPARENT_STREAKS", 0,
+        psMetadataAddF32(config->arguments, PS_LIST_TAIL, "TRANSPARENT_STREAKS", 0,
             "value to adjust excised pixels", transparentStreaks);
         psArgumentRemove(argnum, &argc, argv);
     }
 
-    // if skycells are not provided then we have to execise all pixels  unless -keepnonwarped
+    // if skycells are not provided then we have to execise all pixels  unless -keepnondiffed
     pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist");
 
@@ -449 +510 @@
     }
 
+    if ((argnum = psArgumentGet(argc, argv, "-stats"))) {
+        psArgumentRemove(argnum, &argc, argv);
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "STATS", 0,
+                "name of input stats file", argv[argnum]);
+        psArgumentRemove(argnum, &argc, argv);
+    }
+
     if ((argnum = psArgumentGet(argc, argv, "-tmproot"))) {
         psArgumentRemove(argnum, &argc, argv);
@@ -491 +559 @@
 updateAstrometry(streakFiles *sf)
 {
-    // XXX: why do I check this here? Shouldn't it be just around the call to linearizeTransforms?
     if (sf->bilevelAstrometry) {
-
         if (!linearizeTransforms(sf->astrom)) {
-            // fit failed, leave the astrometry unchanged
+            // fit failed, leave the transform in the file unchanged
             return;
         }
-
-        if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
-            psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
-            streaksExit("", PS_EXIT_UNKNOWN_ERROR);
-        }
-        if (sf->outMask) {
-            pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
-        }
-        if (sf->outWeight) {
-            pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
+    }
+    if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
+        psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
+        streaksExit("", PS_EXIT_UNKNOWN_ERROR);
+    }
+    if (sf->outMask) {
+        pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
+    }
+    if (sf->outWeight) {
+        pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
+    }
+}
+
+static void 
+setStreakBits(psImage *maskImage, psU32 maskStreak)
+{
+    for (int y=0 ; y < maskImage->numRows; y++) {
+        for (int x=0 ; x < maskImage->numCols; x++) {
+            maskImage->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskStreak;
         }
     }
@@ -516 +591 @@
 readAndCopyToOutput(streakFiles *sf, bool exciseAll)
 {
-    bool    updateAstrometry = false;
     if (sf->inImage->pmfile) {
         // image data from pmFPAfile (diff or warp file)
@@ -534 +608 @@
         streaksExit("", PS_EXIT_UNKNOWN_ERROR);
     }
-    // For the chip level files, copy the WCS from the astrometry file to the header
-    // XXX: do we want to do this for raw images as well?
-    if (sf->stage == IPP_STAGE_CHIP) {
-        if (!sf->bilevelAstrometry) {
-            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);
-                streaksExit("", PS_EXIT_UNKNOWN_ERROR);
-            }
-        }
-    }
     sf->outImage->header =  psMemIncrRefCounter(sf->inImage->header);
     if (sf->recImage) {
@@ -577 +640 @@
             }
             addDestreakKeyword(sf->outMask->header);
-            if (updateAstrometry) {
-                pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
-            }
-            setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll);
+            // Note: we don't excise the mask pixels even if exciseAll is true.
+            setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, false);
+            if (exciseAll) {
+                strkGetMaskValues(sf);
+                
+                // add the STREAK bit to the mask image pixels
+                setStreakBits(sf->inMask->image, sf->maskStreak);
+            }
             if (sf->outChMask) {
                 sf->outChMask->header = (psMetadata *) psMemIncrRefCounter(sf->outMask->header);
@@ -605 +672 @@
         }
         addDestreakKeyword(sf->outWeight->header);
-        if (updateAstrometry) {
-            pmAstromWriteWCS(sf->inWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
-        }
         setupImageRefs(sf->outWeight, sf->recWeight, sf->inWeight, sf->extnum, exciseAll);
 
@@ -630 +694 @@
     } else {
         // we have an image cube
-        double initValue;
         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 init value
-        psImage *image = psArrayGet (sf->inImage->imagecube, 0);
-        psMemIncrRefCounter(image);
-        psImageInit(image, initValue);
+        }
+
+        // Now deal with the other output image
         if (exciseImageCube) {
-            sf->outImage->image = image;
-            writeImage(sf->outImage, extname, sf->extnum);
+            // Set the values in the imagecube images to NAN and write them to the output image
+            for (int i = 0; i < psArrayLength(sf->inImage->imagecube); i++) {
+                psImage *image = psArrayGet (sf->inImage->imagecube, i);
+                // XXX: NAN isn't right. It should be the integer equivalent. Should use exciseValue
+                // but it isn't set with this code path. Fix that.
+                psImageInit(image, 65535);
+            }
+            writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
         } else {
-            // write zero valued image to reccovery
             if (sf->recImage) {
-                sf->recImage->image = image;
-                writeImage(sf->recImage, extname, sf->extnum);
+                // Set the values in the imagecube images to zero
+                for (int i = 0; i < psArrayLength(sf->inImage->imagecube); i++) {
+                    psImage *image = psArrayGet (sf->inImage->imagecube, i);
+                    psImageInit(image, 0);
+                }
+                // copy the entire zeroed image to the recovery image
+                writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum);
             }
         }
@@ -683 +751 @@
     }
 
-    float imageValue  = sfiles->inImage->image->data.F32[y][x];
-    if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
-        sfiles->recImage->image->data.F32[y][x] = imageValue;
-    }
-
-    if (sfiles->transparentStreaks == 0) {
-        sfiles->outImage->image->data.F32[y][x] = exciseValue;
+    if (sfiles->inImage->image->type.type == PS_TYPE_U16) {
+        psU16 imageValue  = sfiles->inImage->image->data.U16[y][x];
+        if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
+            sfiles->recImage->image->data.U16[y][x] = imageValue;
+        }
+
+        if (sfiles->transparentStreaks == 0) {
+            sfiles->outImage->image->data.U16[y][x] = exciseValue;
+        } else {
+            if (streak) {
+                // as a visualization aid don't mask the pixel, just change the intensity
+                sfiles->outImage->image->data.U16[y][x] = imageValue + sfiles->transparentStreaks;
+            } else {
+                sfiles->outImage->image->data.U16[y][x] = exciseValue;
+            }
+        }
     } else {
-        if (streak) {
-            // as a visualization aid don't mask the pixel, just change the intensity
-            sfiles->outImage->image->data.F32[y][x] = imageValue + sfiles->transparentStreaks;
+        float imageValue  = sfiles->inImage->image->data.F32[y][x];
+        if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
+            sfiles->recImage->image->data.F32[y][x] = imageValue;
+        }
+
+        if (sfiles->transparentStreaks == 0) {
+            sfiles->outImage->image->data.F32[y][x] = exciseValue;
         } else {
-            sfiles->outImage->image->data.F32[y][x] = exciseValue;
+            if (streak) {
+                // as a visualization aid don't mask the pixel, just change the intensity
+                sfiles->outImage->image->data.F32[y][x] = imageValue + sfiles->transparentStreaks;
+            } else {
+                sfiles->outImage->image->data.F32[y][x] = exciseValue;
+            }
         }
     }
@@ -711 +797 @@
                 sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
         }
-        sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] =
-            sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] | newMaskValue;
-    }
-}
-
-static void
-exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue)
+        sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= newMaskValue;
+    }
+}
+
+static long
+exciseNonDiffedPixels(streakFiles *sfiles, psImageMaskType newMaskValue)
 {
     int cell_x0 = sfiles->astrom->cell_x0;
@@ -725 +810 @@
     int numCols = sfiles->inImage->numCols; // for raw images this was calculated from the width of datasec
     int numRows = sfiles->inImage->numRows; // for raw images this was calculated from the height of datasec
+
+    long excisedPixels = 0;
 
 //    printf("%2d x0: %4d y0: %4d xpar: %d ypar: %d\n", sfiles->extnum, cell_x0, cell_y0, xParity, yParity);
@@ -739 +826 @@
         }
 
-        psU8 *pixels = sfiles->warpedPixels->data.U8[yChip];
+        psU8 *pixels = sfiles->diffedPixels->data.U8[yChip];
 
         if (xParity == 1) {
@@ -746 +833 @@
                 if (! *pixels ) {
                     excisePixel(sfiles, xCell, yCell, false, newMaskValue);
+                    excisedPixels++;
                 }
             }
@@ -755 +843 @@
                 if (!*pixels) {
                     excisePixel(sfiles, xCell, yCell, false, newMaskValue);
+                    excisedPixels++;
                 }
             }
         }
     }
+    return excisedPixels;
 }
 
 static bool
-warpedPixel(streakFiles *sfiles, int x, int y)
+diffedPixel(streakFiles *sfiles, int x, int y)
 {
     PixelPos chipCoord;
 
     if (!CHIP_LEVEL_INPUT(sfiles->stage)) {
-        // if we're here on a skycell image by definition this pixel was warped
+        // if we're here on a skycell image by definition this pixel was diffed
         return true;
     }
@@ -780 +870 @@
     cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, x, y);
 
-    if (chipCoord.x < 0 || chipCoord.x >= sfiles->warpedPixels->numCols) {
+    if (chipCoord.x < 0 || chipCoord.x >= sfiles->diffedPixels->numCols) {
         return false;
     }
-    if (chipCoord.y < 0 || chipCoord.y >= sfiles->warpedPixels->numRows) {
+    if (chipCoord.y < 0 || chipCoord.y >= sfiles->diffedPixels->numRows) {
         return false;
     }
 
-    return psImageGet(sfiles->warpedPixels, chipCoord.x, chipCoord.y) ? true : false;
+    return psImageGet(sfiles->diffedPixels, chipCoord.x, chipCoord.y) ? true : false;
 }
 

@@ -1 +1 @@
 magic
 ssa-core-cpp
-Makefile
 Makefile.bak

@@ -4 +4 @@
 streakscompare
 streaksrelease
-makefile
+Makefile
+Makefile.in
+config.h
+.deps
+streaksVersionDefinitions.h
+config.h.in
+stamp-h1