Changeset 27840 for branches/simtest_nebulous_branches/magic/remove/src

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/Line.c (modified) (23 diffs)
• magic/remove/src/Makefile.am (copied) (copied from trunk/magic/remove/src/Makefile.am )
• magic/remove/src/Makefile.simple (modified) (2 diffs)
• magic/remove/src/diffedpixels.c (copied) (copied from trunk/magic/remove/src/diffedpixels.c )
• magic/remove/src/isdestreaked.c (modified) (3 diffs)
• magic/remove/src/streaksVersion.c (copied) (copied from trunk/magic/remove/src/streaksVersion.c )
• magic/remove/src/streaksVersionDefinitions.h.in (copied) (copied from trunk/magic/remove/src/streaksVersionDefinitions.h.in )
• magic/remove/src/streaksastrom.c (modified) (10 diffs)
• magic/remove/src/streakscompare.c (modified) (5 diffs)
• magic/remove/src/streaksextern.c (modified) (1 diff)
• magic/remove/src/streaksextern.h (modified) (1 diff)
• magic/remove/src/streaksio.c (modified) (23 diffs)
• magic/remove/src/streaksio.h (modified) (1 diff)
• magic/remove/src/streaksrelease.c (modified) (7 diffs)
• magic/remove/src/streaksremove.c (modified) (37 diffs)
• magic/remove/src/streaksremove.h (modified) (4 diffs)
• magic/remove/src/streaksreplace.c (modified) (4 diffs)
• magic/remove/src/streaksutil.c (modified) (1 diff)
• magic/remove/src/warpedpixels.c (deleted)

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/Line.c

@@ -36 +36 @@
     // If the current line is not vertical, check to see if the point
     // is inside the two endpoints independent of direction
-    
+
     if (line->begin.x != line->end.x)
     {
@@ -59 +59 @@
     ux = line->end.x - line->begin.x;
     uy = line->end.y - line->begin.y;
-    
+
     vx = x - line->begin.x;
     vy = y - line->begin.y;
-    
+
     double u_u = ux * ux + uy * uy;                         // (u . u)
     double u_v = ux * vx + uy * vy;                         // (u . v)
-    
+
     if (u_v <= 0) return vx * vx + vy * vy;                 // (v . v)
     if (u_u <= u_v)
@@ -73 +73 @@
         return wx * wx + wy * wy;                           // (w . w)
     }
-    
+
     // Compute P(b) is the base of the perpendicular dropped from tuple to
     // the line
-    
+
     b = u_v / u_u;
     px = vx - b * ux;
@@ -109 +109 @@
     double wx  = first->begin.x - second->begin.x;
     double wy  = first->begin.y - second->begin.y;
-    
+
     // Calculate the perpendicular product, dt
-    
+
     double dt = dx1 * dy2 - dy1 * dx2;
-    
+
     // Are the lines parallel?
-    
+
     if (fabs (dt) <= 1e-8)
     {
         // Check to see if they are overlap
-        
+
         if (dx1 * wy - dy1 * wx != 0 || dx2 * wy - dy2 * wx != 0) return 0;
-        
+
         // The line are coplanar, so check to see if they are degenerate
-        
+
         len1 = dx1 * dx1 + dy1 * dy1;
         len2 = dx2 * dx2 + dy2 * dy2;
-        
+
         // First, check to see if both lines are points and if they are
         // distinct tuples
-        
+
         if (len1 == 0 && len2 == 0)
         {
@@ -137 +137 @@
             return 1;
         }
-        
+
         // Check to see if the first line is a point and is inside the
         // second line
-        
+
         if (len1 == 0)
         {
@@ -148 +148 @@
                 return 0;
         }
-        
+
         // Check to see if the second line is a point and is inside the
         // first line
-        
+
         if (len2 == 0)
         {
@@ -159 +159 @@
                 return 0;
         }
-        
+
         // Since both lines are coplanar and have length, search for
         // overlap tuples
-        
+
         w2x = first->end.x - second->begin.x;
         w2y = first->end.y - second->begin.y;
-        
+
         t0 = (dx2 != 0) ? wx  / dx2 : wy  / dy2;
         t1 = (dx2 != 0) ? w2x / dx2 : w2y / dy2;
-        
+
         if (t0 > t1) SwapDouble (&t0, &t1);
         if ((inFirstSegment || inSecondSegment) && (t0 > 1 || t1 < 0))
@@ -176 +176 @@
         t0 = (t0 < 0) ? 0 : t0;             // Clip to min 0
         t1 = (t1 > 1) ? 1 : t1;             // Clip to max 1
-        
+
         // Set the first tuple of intersection
-        
+
         tuple1->x = second->begin.x + t0 * dx2;
         tuple1->y = second->begin.y + t0 * dy2;
-        
+
         // Check to see if the intersection is a single tuple
-        
+
         if (t0 == t1) return 1;
-        
+
         // Otherwise, set the second tuple of intersection as well
-        
+
         tuple2->x = second->begin.x + t1 * dx2;
         tuple2->y = second->begin.y + t1 * dy2;
         return 2;
     }
-    
+
     // The segments are skew and may intersect in a tuple.
     // Get the intersect parameter for the first line
-    
+
     double i1 = (dx2 * wy - dy2 * wx) / dt;
     if (inFirstSegment && (i1 < 0 || i1 > 1)) return 0;
-    
+
     // Get the intersect parameter for the second line
 
@@ -226 +226 @@
     vertices[3].x = 0;       vertices[3].y = numRows;
 
+    clipLine.begin.x = 0;
+    clipLine.begin.y = 0;
+    clipLine.end.x = 0;
+    clipLine.end.y = 0;
+
     for (i = 0; i < 4 && found < 2; ++i)
     {
@@ -237 +242 @@
                 ++found;
             }
-            else if (tuple1.x != clipLine.begin.x || 
+            else if (tuple1.x != clipLine.begin.x ||
                      tuple1.y != clipLine.begin.y)
             {
@@ -245 +250 @@
         }
     }
-    
+
     // If two endpoints are found, clip the line
-    
+
     if (found > 1)
     {
@@ -283 +288 @@
     vertices[3].x = minX; vertices[3].y = maxY;
 
+    clipLine.begin.x = 0;
+    clipLine.begin.y = 0;
+    clipLine.end.x = 0;
+    clipLine.end.y = 0;
+
     for (i = 0; i < 4 && found < 2; ++i)
     {
@@ -294 +304 @@
                 ++found;
             }
-            else if (tuple1.x != clipLine.begin.x || 
+            else if (tuple1.x != clipLine.begin.x ||
                      tuple1.y != clipLine.begin.y)
             {
@@ -302 +312 @@
         }
     }
-    
+
     // If two endpoints are found, clip the line
-    
+
     if (found > 1)
     {
@@ -360 +370 @@
     @param[out] pixels list of PixelPos pointers corresponding
                        based on the line settings
-    @param[in] line Line to map to pixels   
+    @param[in] line Line to map to pixels
     @param[in] numCols maximum X (columns) for the line
     @param[in] numRows maximum Y (rows) for the line            */
@@ -367 +377 @@
 {
     Line offsetLine;
-    PixelPos *pixel;
     double slope, xOffset, yOffset, xMid, yMid;
     int x, y, xBegin = numCols, yBegin = numRows, xEnd = 0, yEnd = 0;
 
     // Extract the endpoints
-    
+
     double x1 = line->begin.x;
     double y1 = line->begin.y;
     double x2 = line->end.x;
     double y2 = line->end.y;
-    
+
     // Determine the width and height of the line
-    
+
     double dx = x2 - x1;
     double dy = y2 - y1;
@@ -403 +412 @@
 
     // Step point by point based on the dominate axis
-    
+
     if (fabs (dx) > fabs (dy))
     {
@@ -411 +420 @@
         // If line is back to front, reorder endpoints and recompute
         // the line width and height
-        
+
         if (x1 > x2)
         {
@@ -421 +430 @@
 
         // Compute the slope of the line
-        
+
         slope = dy / dx;
-        
+
         // Compute the x and y offsets for the line width extent
 
-        if (xBegin > xEnd) 
+        if (xBegin > xEnd)
             SwapInt (&xBegin, &xEnd);
         else
@@ -451 +460 @@
     }
     else
-    {       
+    {
         // -------------------
         // vertical(ish) lines
@@ -457 +466 @@
         // If line is back to front, reorder endpoints and recompute
         // the line width and height
-        
+
         if (y1 > y2)
         {
@@ -465 +474 @@
             dy = y2 - y1;
         }
-        
+
         // Compute the slope of the line
-        
+
         slope = dx / dy;
-        
+
         // Compute the x and y offsets for the line width extent
 
@@ -481 +490 @@
         xMid = x1 + slope * (yBegin - y1);
         xOffset = fabs (halfWidth * dr / dy);
-        
+
         for (y = yBegin; y != yEnd; ++y)
         {

branches/simtest_nebulous_branches/magic/remove/src/Makefile.simple

@@ -10 +10 @@
 
 REMOVE_OBJECTS=    \
-    ${COMMON_OBJECTS} \
     streaksremove.o \
     streaksextern.o \
-    warpedpixels.o \
+    diffedpixels.o \
     Line.o
 
 REPLACE_OBJECTS=      \
-    ${COMMON_OBJECTS} \
     streaksreplace.o
 
 COMPARE_OBJECTS=      \
-    ${COMMON_OBJECTS} \
     streakscompare.o
 
 RELEASE_OBJECTS=      \
-    ${COMMON_OBJECTS} \
     streaksrelease.o
 
 ISDESTREAKED_OBJECTS=      \
-    ${COMMON_OBJECTS} \
     isdestreaked.o
 
@@ -36 +31 @@
 
 STREAKSFLAGS=-DSTREAKS_COMPRESS_OUTPUT=1
-OPTFLAGS= -g -O2
+OPTFLAGS= -g -O2 -Wall -Werror
 #OPTFLAGS= -g
-CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS} ${STREAKSFLAGS}
-#CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS}
-LDFLAGS=`psmodules-config --libs`
+CFLAGS=`psmodules-config --cflags` `pslib-config --cflags` -std=gnu99 ${OPTFLAGS} ${STREAKSFLAGS}
+#CFLAGS=`psmodules-config --cflags` `pslib-config --cflags` -std=gnu99 ${OPTFLAGS}
+LDFLAGS=`psmodules-config --libs` `pslib-config --libs`
 
 PROGRAMS= streaksremove streaksreplace streakscompare streaksrelease isdestreaked
 HEADERS=Line.h streaksastrom.h streaksextern.h streaksio.h streaksremove.h
 
-all:    ${PROGRAMS}
+all: ${PROGRAMS}
 
-${REMOVE_OBJECTS}:      ${HEADERS}
-streaksremove:  ${REMOVE_OBJECTS}
+# programs all depend on the common objects and the common headers
+PROGRAM_OBJECTS = \
+ ${REMOVE_OBJECTS} \
+ ${REPLACE_OBJECTS} \
+ ${COMPARE_OBJECTS} \
+ ${RELEASE_OBJECTS} \
+ ${ISDESTREAKED_OBJECTS}
 
-streaksreplace:  ${REPLACE_OBJECTS}
+${COMMON_OBJECTS}: ${HEADERS}
 
-streakscompare:  ${COMPARE_OBJECTS}
+${PROGRAM_OBJECTS}: ${HEADERS} ${COMMON_OBJECTS}
 
-streaksrelease:  ${RELEASE_OBJECTS}
+streaksremove:  ${REMOVE_OBJECTS} ${COMMON_OBJECTS}
 
-isdestreaked:   ${ISDESTREAKED_OBJECTS}
+streaksreplace:  ${REPLACE_OBJECTS} ${COMMON_OBJECTS}
+
+streakscompare:  ${COMPARE_OBJECTS} ${COMMON_OBJECTS}
+
+streaksrelease:  ${RELEASE_OBJECTS} ${COMMON_OBJECTS}
+
+isdestreaked:   ${ISDESTREAKED_OBJECTS} ${COMMON_OBJECTS}
 
 

branches/simtest_nebulous_branches/magic/remove/src/isdestreaked.c

@@ -1 +1 @@
 #include "streaksremove.h"
+
+char *streaksProgram = "isdestreaked";
 
 int
@@ -18 +20 @@
     psFits *fits = psFitsOpen(fileName, "r");
     if (!fits) {
-        psError(PS_ERR_UNKNOWN, false, "failed to open fits file: %s\n", fileName);
-        psErrorStackPrint(stderr, "");
+        psErrorStackPrint(stderr, "failed to open fits file: %s\n", fileName);
         exit(PS_EXIT_DATA_ERROR);
     }
@@ -25 +26 @@
     psMetadata *header = psFitsReadHeader(NULL, fits);
     if (!header) {
-        psError(PS_ERR_IO, false, "failed to read header for: %s", fileName);
-        psErrorStackPrint(stderr, "");
+        psErrorStackPrint(stderr, "failed to read header for: %s", fileName);
         exit(PS_EXIT_DATA_ERROR);
     }

branches/simtest_nebulous_branches/magic/remove/src/streaksastrom.c

@@ -93 +93 @@
     // since the transform already accounts for the chip parity we just use the cell parity
     astrom->xParity = (double) xParityCell;
-    astrom->yParity = (double) yParityCell; 
+    astrom->yParity = (double) yParityCell;
 
     astrom->numCols = numCols;
@@ -135 +135 @@
 
 void
-cellToChipInt(int *xChip, int *yChip, strkAstrom *astrom, int xCell, int yCell)
+cellToChipInt(unsigned int *xChip, unsigned int *yChip, strkAstrom *astrom, int xCell, int yCell)
 {
     if (astrom->xParity > 0) {
@@ -148 +148 @@
     }
 }
- 
+
 bool
 SkyToLocal(strkPt *outPt, strkAstrom *astrom, double ra, double dec)
@@ -156 +156 @@
 
     pmFPA *fpa   = (pmFPA *) astrom->fpa;
-    pmChip *chip = (pmChip *) astrom->chip;
 
     // find the RA,DEC coords of the 0,0 pixel for this chip:
@@ -183 +182 @@
 
     pmFPA *fpa   = (pmFPA *) astrom->fpa;
-    pmChip *chip = (pmChip *) astrom->chip;
 
     // find the RA,DEC coords of the 0,0 pixel for this chip:
@@ -251 +249 @@
     pmChip *chip = (pmChip *) astrom->chip;
     pmAstromObj *pt = (pmAstromObj *) astrom->pt;
-    
+
     pt->sky->r = ra;
     pt->sky->d = dec;
@@ -276 +274 @@
     pmChip *chip = (pmChip *) astrom->chip;
     pmAstromObj *pt = (pmAstromObj *) astrom->pt;
-    
+
     cellToChip(&pt->chip->x, &pt->chip->y, astrom, x, y);
 
@@ -294 +292 @@
 }
 
- 
+
 static bool
 readAstrometry(streakFiles *sf)
@@ -362 +360 @@
             break;
         }
-    } 
+    }
     if (!sf->chip) {
         psError(PS_ERR_UNKNOWN, true, "Failed to find chip with data.");
@@ -380 +378 @@
 linearizeTransforms(strkAstrom *astrom)
 {
-    if (!pmAstromLinearizeTransforms((pmFPA *) astrom->fpa, (pmChip *) astrom->chip)) {
-        // streaksExit("linear fit to astrometry failed\n", PS_EXIT_UNKNOWN_ERROR);
-        psError(PS_ERR_UNKNOWN, false, "linear fit to astrometry failed ignoring");
+    if (!pmAstromLinearizeTransforms((pmFPA *) astrom->fpa, (pmChip *) astrom->chip, NULL, NULL, NULL, 0, 0)) {
+        psErrorStackPrint(stderr, "linear fit to astrometry failed. ignoring\n");
         return false;
     }

branches/simtest_nebulous_branches/magic/remove/src/streakscompare.c

@@ -1 +1 @@
 #include "streaksremove.h"
+
+char *streaksProgram = "streakscompare";
 
 static pmConfig * parseArguments(int, char **);
@@ -5 +7 @@
 int main(int argc, char *argv[])
 {
-    long i;
     bool status;
 
@@ -35 +36 @@
     int numErrors = 0;
     for (int component = 0; component < ncomponents; component++) {
-        if (component && !(psFitsMoveExtNum(file1->fits, 1, true) 
+        if (component && !(psFitsMoveExtNum(file1->fits, 1, true)
                            && psFitsMoveExtNum(file2->fits, 1, true))) {
             streaksExit("failed to advance to next extesion\n", PS_EXIT_DATA_ERROR);
@@ -45 +46 @@
         psImage *image2 = file2->image;
 
-        // TODO: do more sanity checking. For example check that extname's  (if any) match 
+        // TODO: do more sanity checking. For example check that extname's  (if any) match
         // check for matching image cubes
         if (image1 && image2) {
@@ -68 +69 @@
                         error = ! isnan(value1) && isnan(value2);
                     }
-                    
+
                     if (error) {
                         numErrors++;

branches/simtest_nebulous_branches/magic/remove/src/streaksextern.c

@@ -78 +78 @@
         }
     }
-    printf("%d streaks on this component\n", streaksOnComponent);
+    fprintf(stderr, "%d streaks on this component\n", streaksOnComponent);
     return pixels;
 }

branches/simtest_nebulous_branches/magic/remove/src/streaksextern.h

@@ -54 +54 @@
             pixels
 */
-    
+
 extern StreakPixels *streak_on_component(Streaks *streaks, strkAstrom *astrom,
                                          int numCols, int numRows);

branches/simtest_nebulous_branches/magic/remove/src/streaksio.c

@@ -12 +12 @@
 // Assumptions about the file structure of non-raw files
 // The 'image' for each file (image, mask weight) is contained in the first
-// extension. 
+// extension.
 
 
@@ -121 +121 @@
     sf->transparentStreaks = psMetadataLookupF64(&status, config->arguments, "TRANSPARENT_STREAKS");
 
+    sf->stats = psMetadataAlloc();
+    psString statsFileName= psMetadataLookupStr(&status, config->arguments, "STATS");
+    if (statsFileName) {
+        sf->statsFile = fopen(statsFileName, "w");
+        if (!sf->statsFile) {
+            psError(PS_ERR_IO, true, "failed to open stats file %s", statsFileName);
+            streaksExit("", PS_EXIT_CONFIG_ERROR);
+        }
+    }
+
     return sf;
 }
@@ -128 +138 @@
 {
     freeImages(sf);
-    psFree(sf->warpedPixels);
+    psFree(sf->diffedPixels);
     psFree(sf->tiles);
     psFree(sf->view);
@@ -288 +298 @@
         sfile->resolved_name = psStringCopy(sfile->pmfile->filename);
 
+        sfile->nHDU = psFitsGetSize(sfile->pmfile->fits);
+
         // copy header from fpu
         sfile->header = (psMetadata*) psMemIncrRefCounter(sfile->pmfile->fpa->hdu->header);
-
-        sfile->nHDU = psFitsGetSize(sfile->pmfile->fits);
 
         return sfile;
@@ -625 +635 @@
         streaksExit("", PS_EXIT_DATA_ERROR);
     }
- 
+
     bool status;
     in->extname = psMetadataLookupStr(&status, in->header, "EXTNAME");
@@ -654 +664 @@
             streaksExit("", PS_EXIT_DATA_ERROR);
         }
-        psImage *image = (psImage *) (in->imagecube->data[0]);
+
+        // Ensure input is of the expected type
+        psDataType expected = isMask ? PS_TYPE_IMAGE_MASK : PS_TYPE_F32; // Expected type for image
+        for (int i = 0; i < in->imagecube->n; i++) {
+            psImage *image = in->imagecube->data[i]; // Image of interest
+            if (image->type.type != expected) {
+                psImage *temp = psImageCopy(NULL, image, expected);
+                psFree(image);
+                in->imagecube->data[i] = temp;
+            }
+        }
     }
     setDataExtent(stage, in, (stage == IPP_STAGE_RAW) && !isMask);
@@ -671 +691 @@
     sfile->fits->options = psFitsOptionsAlloc();
     sfile->fits->options->scaling = PS_FITS_SCALE_MANUAL;
+    sfile->fits->options->fuzz = false;
     sfile->fits->options->bitpix = bitpix;
     sfile->fits->options->bscale = bscale;
@@ -789 +810 @@
         return;
     }
+    streaksVersionHeaderFull(sfile->header);
+
     if (!psFitsWriteImage(sfile->fits, sfile->header, sfile->image, 0, extname)) {
         psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d",
@@ -806 +829 @@
         return;
     }
+    streaksVersionHeaderFull(sfile->header);
     if (!psFitsWriteImageCube(sfile->fits, sfile->header, imagecube, extname)) {
         psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d",
@@ -840 +864 @@
     sfile->header = NULL;
     psFree(sfile->image);
-    sfile->header = NULL;
+    sfile->image = NULL;
     psFree(sfile->imagecube);
-    sfile->header = NULL;
+    sfile->imagecube = NULL;
 }
 
@@ -848 +872 @@
 closeImages(streakFiles *sf)
 {
+    if (sf->stage == IPP_STAGE_RAW) {
+        if (sf->view) {
+            sf->view->readout = -1;
+            pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_AFTER);
+            sf->view->cell = -1;
+            pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_AFTER);
+            sf->view->chip = -1;
+            pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_AFTER);
+        }
+    }
     closeImage(sf->inImage);
     closeImage(sf->outImage);
@@ -902 +936 @@
     bool free_user_copies = true;
     if (user_copies == NULL) {
-        user_copies = "2";
-        free_user_copies = false;
+        // input does not have replication requested
+        return true;
     }
     if (!nebSetXattr(server, outFile->name, "user.copies", user_copies,  NEB_REPLACE)) {
@@ -909 +943 @@
         return false;
     }
-    if (!nebReplicate(server, outFile->name, "any", NULL)) {
+    if (!nebReplicate(server, outFile->name, NULL, NULL)) {
         psError(PM_ERR_UNKNOWN, true, "nebReplicate failed for %s\n%s", outFile->name, nebErr(server));
         return false;
@@ -925 +959 @@
 replicateOutputs(streakFiles *sfiles)
 {
-    bool status = false;
-
     if (!replicate(sfiles->outImage, sfiles->inImage)) {
         psError(PM_ERR_SYS, false, "failed to replicate outImage.");
@@ -973 +1005 @@
 
         if (!nebSwap(server, in->name, out->name)) {
-            psError(PM_ERR_SYS, true, "failed to swap files for: %s.",
+            psError(PM_ERR_SYS, true, "failed to swap files for %s: %s.",
                 in->name, nebErr(server));
             return false;
@@ -1097 +1129 @@
 setMaskedToNAN(streakFiles *sfiles, psU32 maskMask, bool printCounts)
 {
-        int maskedPixels = 0;
-        int nandPixels = 0;
-        int nandWeights = 0;
+        long maskedPixels = 0;
+        long nandPixels = 0;
+        long nandWeights = 0;
 
         psImage *image = sfiles->outImage->image;
+        int imageType = image->type.type;
+        double exciseValue = sfiles->inImage->exciseValue;
+
         psImage *mask = sfiles->inMask->image;
         psImage *weight = NULL;
+        int weightType = 0;
+        double weightExciseValue = NAN;
         if (sfiles->outWeight) {
             weight = sfiles->outWeight->image;
-        }
-        double exciseValue = sfiles->inImage->exciseValue;
+            weightType = weight->type.type;
+            weightExciseValue = sfiles->inWeight->exciseValue;
+        }
 
         if (printCounts) {
@@ -1117 +1155 @@
                 // these gets are not necessary, we could just set the pixels to nan
                 // but I want to get the counts
-                double imageVal  = psImageGet(image, x, y);
-                psU32 maskVal;
+                double imageVal = 0; // avoid compiler warning
+                if (imageType == PS_TYPE_F32) {
+                    imageVal  = image->data.F32[y][x];
+                } else if (imageType == PS_TYPE_S16) {
+                    imageVal = image->data.S16[y][x];
+                } else if (imageType == PS_TYPE_U16) {
+                    imageVal = image->data.U16[y][x];
+                } else {
+                    psError(PS_ERR_PROGRAMMING, true, "unexpected image type found: %d\n",
+                                imageType);
+                    streaksExit("", PS_EXIT_PROG_ERROR);
+                }
+                psU32 maskVal = 0;
                 if (sfiles->stage == IPP_STAGE_RAW) {
-                    int xChip, yChip;
+                    unsigned int xChip, yChip;
                     cellToChipInt(&xChip, &yChip, sfiles->astrom, x, y);
-                    maskVal = psImageGet(mask, xChip, yChip);
+                    maskVal = mask->data.PS_TYPE_IMAGE_MASK_DATA[yChip][xChip];
                 } else {
-                    maskVal = psImageGet(mask, x, y);
+                    maskVal = mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
                 }
                 if (maskVal & maskMask) {
                     ++maskedPixels;
-                    if (!isExciseValue(imageVal, sfiles->inImage->exciseValue)) {
+                    if (!isExciseValue(imageVal, exciseValue)) {
                         ++nandPixels;
-                        psImageSet(image, x, y, exciseValue);
+                        if (imageType == PS_TYPE_F32) {
+                            image->data.F32[y][x] = exciseValue;
+                        } else if (imageType == PS_TYPE_S16) {
+                            image->data.S16[y][x] = exciseValue;
+                        } else if (imageType == PS_TYPE_U16) {
+                            image->data.U16[y][x] = exciseValue;
+                        } else {
+                            psError(PS_ERR_PROGRAMMING, true, "unexpected image type found: %d\n",
+                                imageType);
+                            streaksExit("", PS_EXIT_PROG_ERROR);
+                        }
                     }
                     if (weight) {
-                        double weightVal = weight ? psImageGet(weight, x, y) : 0;
-                        if (!isnan(weightVal)) {
-                            ++nandWeights;
-                            psImageSet(weight, x, y, NAN);
+                        if (weightType == PS_TYPE_F32) {
+                            double weightVal = weight->data.F32[y][x];
+                            if (!isnan(weightVal)) {
+                                ++nandWeights;
+                                weight->data.F32[y][x] = NAN;
+                            }
+                        } else if(weightType == PS_TYPE_S16) {
+                            double weightVal = weight->data.S16[y][x];
+                            if (!isExciseValue(weightVal, weightExciseValue)) {
+                                ++nandWeights;
+                                image->data.S16[y][x] = weightExciseValue;
+                            }
+                        } else if (weightType == PS_TYPE_U16) {
+                            double weightVal = weight->data.U16[y][x];
+                            if (!isExciseValue(weightVal, weightExciseValue)) {
+                                ++nandWeights;
+                                image->data.U16[y][x] = weightExciseValue;
+                            }
+                        } else {
+                            psError(PS_ERR_PROGRAMMING, true, "unexpected image type found: %d\n",
+                                    weightType);
+                                streaksExit("", PS_EXIT_PROG_ERROR);
                         }
                     }
@@ -1144 +1221 @@
         if (printCounts) {
             psLogMsg(sfiles->program_name, PS_LOG_INFO, "time to NAN mask pixels: %f\n", psTimerClear("NAN_MASKED"));
-            int totalPixels = image->numRows * image->numCols;
+            long totalPixels = image->numRows * image->numCols;
             psLogMsg(sfiles->program_name, PS_LOG_INFO, "pixels:        %10ld\n", totalPixels);
             psLogMsg(sfiles->program_name, PS_LOG_INFO, "masked pixels: %10ld %4.2f%%\n", maskedPixels, 100. * maskedPixels / totalPixels);
@@ -1167 +1244 @@
         psMetadataAddU16(in->header, PS_LIST_TAIL, "BLANK", 0, "", 65535);
         psMetadataAddU16(in->header, PS_LIST_TAIL, "ZBLANK", 0, "", 65535);
-    } else {
+    } else if (in->image->type.type == PS_TYPE_S16) {
+        in->exciseValue = 32767;
+        psMetadataAddU16(in->header, PS_LIST_TAIL, "BLANK", 0, "", 32767);
+        psMetadataAddU16(in->header, PS_LIST_TAIL, "ZBLANK", 0, "", 32767);
+    } else if (in->image->type.type == PS_TYPE_F32) {
         in->exciseValue = NAN;
-    }
-}
-
-void
-strkGetMaskValues(streakFiles *sfiles, psU32 *maskStreak, psU32 *maskMask)
-{
+    } else {
+        psError(PS_ERR_PROGRAMMING, true, "unexpected image type found: %d\n", in->image->type.type);
+        streaksExit("", PS_EXIT_PROG_ERROR);
+    }
+}
+
+// Get the mask values that we need.
+// If an input mask file is provided get them from there.
+// Otherwise get them from the recipe MASKS
+void
+strkGetMaskValues(streakFiles *sfiles)
+{
+    if (sfiles->maskStreak) {
+        // already done
+        return;
+    }
     if (sfiles->inMask && sfiles->inMask->header) {
         if (!pmConfigMaskReadHeader(sfiles->config, sfiles->inMask->header)) {
@@ -1187 +1278 @@
         streaksExit("", PS_EXIT_CONFIG_ERROR);
     }
-    *maskStreak = psMetadataLookupU32(&status, masks, "STREAK");
+    sfiles->maskStreak = psMetadataLookupU32(&status, masks, "STREAK");
     if (!status) {
         psError(PM_ERR_CONFIG, false, "failed to lookup mask value for STREAK in recipes\n");
@@ -1203 +1294 @@
         }
     }
-    *maskMask = ~convPoor;
+    sfiles->maskMask = ~convPoor;
 }
 
@@ -1215 +1306 @@
         if (SFILE_IS_IMAGE(in)) {
             // Turn off compression (code adapted from pmFPAWrite)
-#ifdef SAVE_AND_RESTORE_COMPRESSION
             int bitpix = out->fits->options ? out->fits->options->bitpix : 0; // Desired bits per pixel
             psFitsCompression *compress = psFitsCompressionGet(out->fits); // Current compression options
-#endif
             if (!psFitsSetCompression(out->fits, PS_FITS_COMPRESS_NONE, NULL, 0, 0, 0)) {
                 psError(PM_ERR_UNKNOWN, false, "failed to turn off compression for extension %d\n", extnum);
                 streaksExit("", PS_EXIT_UNKNOWN_ERROR);
             }
-#ifdef SAVE_AND_RESTORE_COMPRESSION
             if (out->fits->options) {
                 out->fits->options->bitpix = 0;
@@ -1238 +1326 @@
                 streaksExit("", PS_EXIT_UNKNOWN_ERROR);
             }
-#endif
         } else {
             copyTable(out, in, extnum);

branches/simtest_nebulous_branches/magic/remove/src/streaksio.h

@@ -16 +16 @@
 void copyFitsOptions(sFile *out, sFile *rec, sFile *in, psVector *tiles);
 void setupImageRefs(sFile *out, sFile *recoveryOut, sFile *in, int extnum, bool exciseAll);
-void strkGetMaskValues(streakFiles *sfiles, psU32 *maskStreak, psU32 *maskMask);
+void strkGetMaskValues(streakFiles *sfiles);
 void copyExtraExtensions(streakFiles *sfiles);
 

branches/simtest_nebulous_branches/magic/remove/src/streaksrelease.c

@@ -3 +3 @@
 static pmConfig *parseArguments(int argc, char **argv);
 static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll);
-static StreakPixels * getStreakPixels(streakFiles *sfiles, Streaks *streaks);
-static void exciseNonWarpedPixels(streakFiles *sfiles, double newMaskValue);
-static bool warpedPixel(streakFiles *sfiles, PixelPos *cellCoord);
-static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, double newMaskValue);
 static void writeImages(streakFiles *sf, bool exciseImageCube);
+
+char *streaksProgram = "streaksrelease";
 
 int
 main(int argc, char *argv[])
 {
-    bool status;
-
     psLibInit(NULL);
     psTimerStart("STREAKSREMOVE");
@@ -22 +18 @@
         return 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 
 
     // Does true work here?
@@ -56 +48 @@
 
         // now that we've read the input files, lookup the mask values that we read
-        if (maskStreak == 0) {
-            strkGetMaskValues(sfiles, &maskStreak, &maskMask);
-        }
-
-        setMaskedToNAN(sfiles, maskMask, true);
+        strkGetMaskValues(sfiles);
+
+        if (!sfiles->inImage->imagecube) {
+            setMaskedToNAN(sfiles, sfiles->maskMask, true);
+        }
 
         // write out the destreaked temporary images and the recovery images
@@ -216 +208 @@
         // image data directly from psFits
         readImage(sf->inImage, sf->extnum, sf->stage, false);
-        
+
         // astrom struct is only needed for raw stage, and only the chip to cell parameters are used
         sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, NULL, NULL, false,
@@ -232 +224 @@
     if (sf->inImage->image) {
         setupImageRefs(sf->outImage, sf->recImage, sf->inImage, sf->extnum, exciseAll);
-    } else if (sf->inImage->imagecube) {
-        // Image cubes should have been excised in the destreaking process
-        streaksExit("unexpected imagecube found", PS_EXIT_CONFIG_ERROR);
     }  else {
-        return false;
+        // image cubes are handled specially
     }
 
@@ -284 +273 @@
 
 
+// XXXX: why isn't this in streaksio.c ?? See also streakremove.c
 static void
 writeImages(streakFiles *sf, bool exciseImageCube)
@@ -304 +294 @@
             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);
-        if (exciseImageCube) {
-            sf->outImage->image = image;
-            writeImage(sf->outImage, extname, sf->extnum);
-        } else {
-            // 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;

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);
-    }
-}
+}

branches/simtest_nebulous_branches/magic/remove/src/streaksremove.h

@@ -33 +33 @@
     int         yParity;
     double      exciseValue;
+    psString    program;
 } sFile;
 
-
-typedef enum {
-    IPP_STAGE_NONE = 0,
-    IPP_STAGE_RAW,
-    IPP_STAGE_CHIP,
-    IPP_STAGE_WARP,
-    IPP_STAGE_DIFF
-} ippStage;
-
+// used for error messages
+extern char * streaksProgram;
 
 typedef struct {
@@ -64 +58 @@
     sFile *inSources;
     sFile *outSources;
+    psMetadata *stats;
+    FILE *statsFile;
     psString class_id;
     pmFPAfile  *inAstrom;
@@ -71 +67 @@
     pmChip  *chip;  // current chip
     pmCell  *cell;  // current cell
-    psImage *warpedPixels;
+    psImage *diffedPixels;
     psVector    *tiles;
     double  transparentStreaks;
     psString    program_name;
+    psU32   maskStreak;
+    psU32   maskMask;
 } streakFiles;
 
@@ -84 +82 @@
 // can't declare this in streaksastrom due to header file ordering
 extern void cellToChip(double *xChip, double *yChip, strkAstrom *astrom, double xCell, double yCell);
-extern void cellToChipInt(int *xChip, int *yChip, strkAstrom *astrom, int xCell, int yCell);
+extern void cellToChipInt(unsigned int *xChip, unsigned int *yChip, strkAstrom *astrom, int xCell, int yCell);
 
-extern bool computeWarpedPixels(streakFiles *sf);
+extern bool computeDiffedPixels(streakFiles *sf);
 extern void streaksExit(psString, int);
 extern ippStage parseStage(psString);
 extern psString pathToDirectory(char *path);
 extern void setStreakFiles( streakFiles *);
+
+extern bool streaksVersionHeaderFull(psMetadata *header);
+extern psString streaksVersionLong(void);
 
 #define CHIP_LEVEL_INPUT(_stage) ((_stage == IPP_STAGE_RAW) || (_stage == IPP_STAGE_CHIP))

branches/simtest_nebulous_branches/magic/remove/src/streaksreplace.c

@@ -6 +6 @@
 static bool readAndCopyToOutput(streakFiles *sf, bool restoreImageCube);
 
+char *streaksProgram = "streaksreplace";
+
 int main(int argc, char *argv[])
 {
-    long i;
     bool status;
-    StreakPixels *pixels;
-    PixelPos *pixelPos;
 
     psLibInit(NULL);
@@ -92 +91 @@
 
     if (!replicateOutputs(sfiles)) {
-        psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
-        psErrorStackPrint(stderr, "");
+        psErrorStackPrint(stderr, "failed to replicate output files");
         exit(PS_EXIT_UNKNOWN_ERROR);
     }
@@ -172 +170 @@
             true);
     }
-    
+
     bool gotMask = false;
     if ((argnum = psArgumentGet(argc, argv, "-mask"))) {
@@ -353 +351 @@
         // we have an image cube
         double initValue;
-        // otherwise write it to the output 
+        // otherwise write it to the output
         writeImageCube(sf->outImage, sf->recImage->imagecube, extname, sf->extnum);
 

branches/simtest_nebulous_branches/magic/remove/src/streaksutil.c

@@ -44 +44 @@
 // we just bail out
 void streaksExit(psString str, int exitCode) {
-    psErrorStackPrint(stderr, str);
+    psErrorStackPrint(stderr, "%s", str);
     if (ourStreakFiles) {
         deleteTemps(ourStreakFiles);

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