Changeset 27840 for branches/simtest_nebulous_branches/magic/remove/src/streaksremove.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_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/simtest_nebulous_branches/magic

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

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

branches/simtest_nebulous_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/simtest_nebulous_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, double newMaskValue);
-static bool warpedPixel(streakFiles *sfiles, int x, int y);
-static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, double newMaskValue);
+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, psU32 maskStreak);
-
+static void censorSources(streakFiles *sfiles, psImageMaskType 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[])
@@ -25 +32 @@
 
     psLibInit(NULL);
-    psTimerStart("STREAKSREMOVE");
+    psTimerStart("TOTAL_TIME");
 
     pmConfig *config = parseArguments(argc, argv);
@@ -33 +40 @@
     }
 
-    // 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");
 
@@ -42 +45 @@
     Streaks *streaks = readStreaksFile(streaksFileName);
     if (!streaks) {
-        psError(PS_ERR_UNKNOWN, "failed to read streaks file: %s", streaksFileName);
+        psError(PS_ERR_UNKNOWN, false, "failed to read streaks file: %s", streaksFileName);
         streaksExit("", PS_EXIT_PROG_ERROR);
     }
@@ -49 +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
@@ -60 +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
@@ -75 +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);
+    }
+
     if (sfiles->stage == IPP_STAGE_RAW) {
         // Except for raw stage, all of our (GPC1) files have one image extension.
@@ -97 +102 @@
     }
 
-    int totalPixels = 0;
-    int totalStreakPixels = 0;
-
+    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 {
@@ -117 +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;
@@ -131 +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"));
-            
+
             // 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) {
-                    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");
 
-                for (int y=0 ; y < sfiles->inImage->numRows; y++) {
-                    for (int x = 0; x < sfiles->inImage->numCols; x++) {
-                        if (psImageGet(pixels, x, y)) {
-                            ++totalStreakPixels;
-                            if (!checkNonWarpedPixels || warpedPixel(sfiles, x, y)) {
-
-                                excisePixel(sfiles, x, y, true, maskStreak);
-
-                            } else {
-                                // This pixel was not included in any warp and has thus already excised
-                                // by exciseNonWarpedPixels
-                            }
-                        }
-                    }
-                }
-
-                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);
                 }
 
-            } else { 
+            } else {
                 // this component contains an image cube
                 // For now excise it completely
@@ -184 +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));
 
@@ -199 +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
@@ -208 +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)) {
-        psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
+        psErrorStackPrint(stderr, "failed to replicate output files");
         deleteTemps(sfiles);
-        psErrorStackPrint(stderr, "");
         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.
@@ -229 +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");
@@ -239 +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;
         }
     }
@@ -247 +294 @@
     pmConceptsDone();
     pmModelClassCleanup();
-    streaksNebulousCleanup(); 
+    streaksNebulousCleanup();
     pmConfigDone();
-    psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE"));
     psLibFinalize();
 
@@ -255 +301 @@
 
     return 0;
+}
+
+static long
+censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonDiffedPixels, psU16 maskStreak)
+{
+    long streakPixels = 0;
+
+    for (int y=0 ; y < sfiles->inImage->numRows; y++) {
+        for (int x = 0; x < sfiles->inImage->numCols; x++) {
+            if (psImageGet(pixels, x, y)) {
+                if (!checkNonDiffedPixels || diffedPixel(sfiles, x, y)) {
+                    ++streakPixels;
+
+                    excisePixel(sfiles, x, y, true, maskStreak);
+
+                } else {
+                    // This pixel was not included in any warp and has thus already excised
+                    // by exciseNonDiffedPixels
+                }
+            }
+        }
+    }
+    return streakPixels;
 }
 
@@ -266 +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");
@@ -356 +425 @@
             true);
     }
-    
-    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);
     }
 
@@ -370 +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");
 
@@ -441 +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);
@@ -483 +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;
         }
     }
@@ -508 +591 @@
 readAndCopyToOutput(streakFiles *sf, bool exciseAll)
 {
-    bool    updateAstrometry = false;
     if (sf->inImage->pmfile) {
         // image data from pmFPAfile (diff or warp file)
@@ -526 +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) {
@@ -569 +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);
@@ -597 +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);
 
@@ -622 +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 
+            // 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);
             }
         }
@@ -666 +742 @@
 
 static void
-excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, double newMaskValue)
+excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, psImageMaskType newMaskValue)
 {
     double exciseValue = sfiles->inImage->exciseValue;
@@ -675 +751 @@
     }
 
-    double imageValue  = psImageGet (sfiles->inImage->image,  x, y);
-    if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
-        psImageSet (sfiles->recImage->image,  x, y, imageValue);
-    }
-
-    if (sfiles->transparentStreaks == 0) {
-        psImageSet (sfiles->outImage->image,  x, y, 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
-            psImageSet (sfiles->outImage->image,  x, y, 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 {
-            psImageSet (sfiles->outImage->image,  x, y, 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;
+            }
         }
     }
@@ -693 +787 @@
     if (sfiles->outWeight) {
         if (sfiles->recWeight) {
-            double weightValue = psImageGet (sfiles->inWeight->image, x, y);
-            psImageSet (sfiles->recWeight->image, x, y, weightValue);
+            sfiles->recWeight->image->data.F32[y][x] = sfiles->inWeight->image->data.F32[y][x];
         }
         // Assume that weight images are always a floating point type
-        psImageSet (sfiles->outWeight->image, x, y, NAN);
+        sfiles->outWeight->image->data.F32[y][x] = NAN;
     }
     if (sfiles->outMask) {
         if (sfiles->recMask) {
-            double maskValue   = psImageGet (sfiles->inMask->image,   x, y);
-            psImageSet (sfiles->recMask->image,   x, y, maskValue);
-        }
-        psImageSet (sfiles->outMask->image,   x, y, newMaskValue);
-    }
-}
-
-static void
-exciseNonWarpedPixels(streakFiles *sfiles, double newMaskValue)
+            sfiles->recMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] =
+                sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
+        }
+        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;
@@ -717 +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);
@@ -731 +826 @@
         }
 
-        psU8 *pixels = sfiles->warpedPixels->data.U8[yChip];
+        psU8 *pixels = sfiles->diffedPixels->data.U8[yChip];
 
         if (xParity == 1) {
@@ -738 +833 @@
                 if (! *pixels ) {
                     excisePixel(sfiles, xCell, yCell, false, newMaskValue);
+                    excisedPixels++;
                 }
             }
@@ -747 +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;
     }
@@ -772 +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;
 }
 
@@ -785 +883 @@
 // streak mask
 static void
-censorSources(streakFiles *sfiles, psU32 maskStreak)
+censorSources(streakFiles *sfiles, psImageMaskType maskStreak)
 {
     if ((!sfiles->inSources) || (!sfiles->outMask)) {
@@ -799 +897 @@
     sFile *out = sfiles->outSources;
 
-    in->header = psFitsReadHeader(NULL, in->fits);
-    if (!in->header) {
-        psError(PS_ERR_IO, false, "failed to read header from %s", in->resolved_name);
+
+    // Primary header, should be "something.hdr"
+    {
+        psMetadata *header = psFitsReadHeader(NULL, in->fits);
+        if (!header) {
+            psError(PS_ERR_IO, false, "failed to read header from %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+
+        bool status;
+        psString extname = psMetadataLookupStr(&status, header, "EXTNAME");
+        if (!extname) {
+            psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        addDestreakKeyword(header);
+
+        if (!psFitsWriteBlank(out->fits, header, extname)) {
+            psError(PS_ERR_IO, false, "failed to write blank in header of %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        psFree(header);
+    }
+
+    // Extension with PSF fits, should be "something.psf"
+    {
+        if (!psFitsMoveExtNum(in->fits, 1, true)) {
+            psErrorStackPrint(stderr, "failed to read header from %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+
+        psMetadata *header = psFitsReadHeader(NULL, in->fits);
+        if (!header) {
+            psErrorStackPrint(stderr, "failed to read header from %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        psString extname = psMetadataLookupStr(NULL, header, "EXTNAME");
+        if (!extname) {
+            psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+
+        psArray *inTable = psFitsReadTable(in->fits);
+        if (!inTable->n) {
+            psErrorStackPrint(stderr, "table in %s is empty", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+
+        psArray *outTable = psArrayAllocEmpty(inTable->n);
+        int j = 0;
+        int numCensored = 0;
+        for (int i = 0 ; i < inTable->n; i++) {
+            psMetadata *row = inTable->data[i];
+
+            psF32 x = psMetadataLookupF32(NULL, row, "X_PSF");
+            psF32 y = psMetadataLookupF32(NULL, row, "Y_PSF");
+
+            psImageMaskType mask;
+            if ((x >= maskImage->numCols) || (y >= maskImage->numRows) || (x <  0) || (y < 0)) {
+                mask = maskStreak;
+            } else {
+                mask = maskImage->data.PS_TYPE_IMAGE_MASK_DATA[(int)y][(int)x];
+            }
+
+            // Key the source if the center pixel is not masked with maskStreak
+            if (!(mask & maskStreak) ) {
+                psArraySet(outTable, j++, row);
+            } else {
+                numCensored++;
+            }
+        }
+
+        // get rid of unused elements (don't know if this is necessary)
+        psArrayRealloc(outTable, j);
+
+        addDestreakKeyword(header);
+        if (psArrayLength(outTable) > 0) {
+            printf("Censored %d sources\n", numCensored);
+            if (! psFitsWriteTable(out->fits, header, outTable, extname)) {
+                psErrorStackPrint(stderr, "failed to write table to %s", out->resolved_name);
+                streaksExit("", PS_EXIT_DATA_ERROR);
+            }
+        } else {
+            printf("Censored ALL %d sources\n", numCensored);
+            if (! psFitsWriteTableEmpty(out->fits, header, inTable->data[0], extname)) {
+                psErrorStackPrint(stderr, "failed to write empty table to %s", out->resolved_name);
+                streaksExit("", PS_EXIT_DATA_ERROR);
+            }
+        }
+        psFree(header);
+        psFree(outTable);
+        psFree(inTable);
+    }
+
+    // XXX Will need to update to handle extension with extended sources, etc.
+
+    if (!psFitsClose(out->fits)) {
+        psErrorStackPrint(stderr, "failed to close table %s", out->resolved_name);
         streaksExit("", PS_EXIT_DATA_ERROR);
     }
-
-    bool status;
-    psString extname = psMetadataLookupStr(&status, in->header, "EXTNAME");
-    if (!extname) {
-        psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-
-    psArray *inTable = psFitsReadTable(in->fits);
-    if (!inTable->n) {
-        psError(PS_ERR_IO, false, "table in %s is empty", in->resolved_name);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-
-    psArray *outTable = psArrayAllocEmpty(inTable->n);
-    int j = 0;
-    int numCensored = 0;
-    for (int i = 0 ; i < inTable->n; i++) {
-        psMetadata *row = inTable->data[i];
-
-        psF32 x = psMetadataLookupF32 (&status, row, "X_PSF");
-        psF32 y = psMetadataLookupF32 (&status, row, "Y_PSF");
-        
-        psU32 mask = psImageGet(maskImage, x, y);
-
-        // Key the source if the center pixel is not masked with maskStreak
-        if (! (mask & maskStreak) ) {
-            psArraySet(outTable, j++, row);
-        } else {
-            numCensored++;
-        }
-    }
-
-    // get rid of unused elements (don't know if this is necessary)
-    psArrayRealloc(outTable, j);
-
-    addDestreakKeyword(in->header);
-    if (psArrayLength(outTable) > 0) {
-        printf("Censored %d sources\n", numCensored);
-        if (! psFitsWriteTable(out->fits, in->header, outTable, extname)) {
-            psError(PS_ERR_IO, false, "failed to write table to %s", out->resolved_name);
-            streaksExit("", PS_EXIT_DATA_ERROR);
-        }
-    } else {
-        printf("Censored ALL %d sources\n", numCensored);
-        if (! psFitsWriteTableEmpty(out->fits, in->header, inTable->data[0], extname)) {
-            psError(PS_ERR_IO, false, "failed to write empty table to %s", out->resolved_name);
-            streaksExit("", PS_EXIT_DATA_ERROR);
-        }
-    }
-
-    if (!psFitsClose(out->fits)) {
-        psError(PS_ERR_IO, false, "failed to close table %s", out->resolved_name);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-}
+}

@@ -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