Changeset 25027 for branches/pap/magic/remove/src

Timestamp:

Aug 7, 2009, 4:08:25 PM (17 years ago)

Author:

Paul Price

Message:

Merging trunk (r25026) to get up-to-date on old branch.

Location:

branches/pap

Files:

• . (modified) (1 prop)
• magic (modified) (1 prop)
• magic/remove/src (modified) (1 prop)
• magic/remove/src/Line.c (modified) (5 diffs)
• magic/remove/src/Line.h (modified) (2 diffs)
• magic/remove/src/Makefile.simple (modified) (3 diffs)
• magic/remove/src/isdestreaked.c (copied) (copied from trunk/magic/remove/src/isdestreaked.c )
• magic/remove/src/streaksastrom.c (modified) (2 diffs)
• magic/remove/src/streaksastrom.h (modified) (1 diff)
• magic/remove/src/streaksextern.c (modified) (2 diffs)
• magic/remove/src/streaksextern.h (modified) (2 diffs)
• magic/remove/src/streaksio.c (modified) (31 diffs)
• magic/remove/src/streaksio.h (modified) (2 diffs)
• magic/remove/src/streaksrelease.c (modified) (7 diffs)
• magic/remove/src/streaksremove.c (modified) (24 diffs)
• magic/remove/src/streaksremove.h (modified) (4 diffs)
• magic/remove/src/streaksreplace.c (modified) (5 diffs)
• magic/remove/src/streaksutil.c (modified) (1 diff)

branches/pap/magic/remove/src

@@ +1 @@
+isdestreaked
 streaksremove
 streaksreplace
 streakscompare
 streaksrelease
+makefile

@@ +1 @@
+isdestreaked
 streaksremove
 streaksreplace
 streakscompare
 streaksrelease
+makefile

branches/pap/magic/remove/src/Line.c

@@ -16 +16 @@
 {
     double temp = *first;
+    *first = *second;
+    *second = temp;
+}
+
+/** Internal routine to swap integer values */
+
+void SwapInt (int* first, int* second)
+{
+    int temp = *first;
     *first = *second;
     *second = temp;
@@ -255 +264 @@
 }
 
-/** Map a line to an image for its specified width and store as a list
-    of pixel positions
+/** Clip the line between (minX,minY) and (maxX,maxY)
+
+    @param[in,out] line line to be clipped within the bounds
+    @param[in] minX minimum X (columns) for the line
+    @param[in] minY minimum Y (rows) for the line
+    @param[in] maxX maximum X (columns) for the line
+    @param[in] maxY maximum Y (rows) for the line
+    @return true if line overlaps the clip boundaries           */
+
+bool LineClipFull (Line *line, int minX, int minY, int maxX, int maxY)
+{
+    unsigned int i, found = 0;
+    Line boundLine, clipLine;
+    strkPt tuple1, tuple2, vertices[4];
+    vertices[0].x = minX; vertices[0].y = minY;
+    vertices[1].x = maxX; vertices[1].y = minY;
+    vertices[2].x = maxX; vertices[2].y = maxY;
+    vertices[3].x = minX; vertices[3].y = maxY;
+
+    for (i = 0; i < 4 && found < 2; ++i)
+    {
+        boundLine.begin = vertices[i];
+        boundLine.end   = vertices[(i + 1) % 4];
+        if (LineIntercept (line, &boundLine, &tuple1, &tuple2, false, true))
+        {
+            if (found == 0)
+            {
+                clipLine.begin = tuple1;
+                ++found;
+            }
+            else if (tuple1.x != clipLine.begin.x || 
+                     tuple1.y != clipLine.begin.y)
+            {
+                clipLine.end = tuple1;
+                ++found;
+            }
+        }
+    }
+    
+    // If two endpoints are found, clip the line
+    
+    if (found > 1)
+    {
+        if (clipLine.begin.x <= clipLine.end.x)
+        {
+            line->begin = clipLine.begin;
+            line->end   = clipLine.end;
+        }
+        else
+        {
+            line->begin = clipLine.end;
+            line->end   = clipLine.begin;
+        }
+    }
+    return found > 1;
+}
+
+/** Move a line by the specified X and Y offsets
+    @param[in,out] line Line to move by X and Y offsets
+    @param[in] xOffset X shift applied to both endpoints
+    @param[in] yOffset Y shift applied to both endpoints        */
+
+void LineMove (Line *line, double xOffset, double yOffset)
+{
+    line->begin.x += xOffset;
+    line->begin.y += yOffset;
+    line->end.x   += xOffset;
+    line->end.y   += yOffset;
+}
+
+/** Return the maximum bounds between the line endpoints and
+    current bounds
+    @param[in] line Line endpoints to compare
+    @param[in,out] xMin minimum X value to update
+    @param[in,out] xMax maximum X value to update
+    @param[in,out] yMin minimum Y value to update
+    @param[in,out] yMax maximum Y value to update               */
+
+void MaxBounds (Line *line, int *xMin, int *xMax,  int *yMin, int *yMax)
+{
+    if (line->begin.x < *xMin) *xMin = (int) floor (line->begin.x);
+    if (line->end.x   < *xMin) *xMin = (int) floor (line->end.x);
+    if (line->begin.y < *yMin) *yMin = (int) floor (line->begin.y);
+    if (line->end.y   < *yMin) *yMin = (int) floor (line->end.y);
+
+    if (line->begin.x > *xMax) *xMax = (int) ceil (line->begin.x);
+    if (line->end.x   > *xMax) *xMax = (int) ceil (line->end.x);
+    if (line->begin.y > *yMax) *yMax = (int) ceil (line->begin.y);
+    if (line->end.y   > *yMax) *yMax = (int) ceil (line->end.y);
+}
+
+/** Map a line to an image for its specified width and store as
+    a list of pixel positions
 
     @param[out] pixels list of PixelPos pointers corresponding
                        based on the line settings
-    @param[in] line Line to map to pixels                       */
-
-void PixelsFromLine (StreakPixels* pixels, Line *line)
-{
+    @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            */
+
+void PixelsFromLine (StreakPixels* pixels, Line *line, int numCols, int numRows)
+{
+    Line offsetLine;
     PixelPos *pixel;
-    double slope, xOffset, yOffset, xMid, yMid, xBegin, yBegin, xEnd, yEnd, x, y;
+    double slope, xOffset, yOffset, xMid, yMid;
+    int x, y, xBegin = numCols, yBegin = numRows, xEnd = 0, yEnd = 0;
 
     // Extract the endpoints
@@ -280 +384 @@
     double dr = sqrt (dx * dx + dy * dy);
     double halfWidth  = line->width / 2.0;
-    double halfWidth2 = halfWidth * halfWidth;
     if (!dr) return;
-    
+
+    // Compute the intercepts of line width bounds and determine maximum
+    // bounds in each axis
+
+    xOffset = -halfWidth * dy / dr;
+    yOffset =  halfWidth * dx / dr;
+
+    offsetLine = *line;
+    LineMove (&offsetLine, xOffset, yOffset);
+    if (LineClip (&offsetLine, numCols, numRows))
+        MaxBounds (&offsetLine, &xBegin, &xEnd, &yBegin, &yEnd);
+
+    offsetLine = *line;
+    LineMove (&offsetLine, -xOffset, -yOffset);
+    if (LineClip (&offsetLine, numCols, numRows))
+        MaxBounds (&offsetLine, &xBegin, &xEnd, &yBegin, &yEnd);
+
     // Step point by point based on the dominate axis
     
@@ -307 +426 @@
         // Compute the x and y offsets for the line width extent
 
-        xOffset = halfWidth * dy / dr;
-        yOffset = halfWidth * dr / dx;
-        yMid   = y1 + slope * (floor (x1 - xOffset) - x1);
-        xBegin = floor (x1 - xOffset);
-        xEnd   = ceil  (x2 + xOffset) + 1.0;
-
-        for (x = xBegin; x < xEnd; ++x)
-        {
-            yBegin = floor (yMid - yOffset);
-            yEnd   = ceil  (yMid + yOffset) + 1.0;
-            for (y = yBegin; y < yEnd; ++y)
-            {
-                if (DistanceSquared (line, x, y) <= halfWidth2)
+        if (xBegin > xEnd) 
+            SwapInt (&xBegin, &xEnd);
+        else
+            ++xEnd;
+        if (xBegin < 0) xBegin = 0;
+        if (xEnd > numCols) xEnd = numCols;
+
+        yMid = y1 + slope * (xBegin - x1);
+        yOffset = fabs (halfWidth * dr / dx);
+
+        for (x = xBegin; x != xEnd; ++x)
+        {
+            yBegin = (int) floor (yMid - yOffset);
+            yEnd   = (int) ceil  (yMid + yOffset) + 1;
+            for (y = yBegin; y != yEnd; ++y)
+            {
+                if (y >= 0 && y < numRows)
                 {
-                    if (x >=0 && y >= 0) {
-                        pixel = psAlloc (sizeof(PixelPos));
-                        pixel->x = (unsigned int) x;
-                        pixel->y = (unsigned int) y;
-                        psArrayAdd (pixels, 1024, pixel);
-                        psFree (pixel);
-                    }
+                    psImageSet(pixels, x, y, 1);
                 }
             }
@@ -354 +471 @@
         
         // Compute the x and y offsets for the line width extent
-        
-        xOffset = halfWidth * dr / dy;
-        yOffset = halfWidth * dx / dr;
-
-        xMid   = x1 + slope * (floor (y1 - yOffset) - y1);
-        yBegin = floor (y1 - yOffset);
-        yEnd   = ceil  (y2 + yOffset) + 1.0;
-        
-        for (y = yBegin; y < yEnd; ++y)
-        {
-            xBegin = floor (xMid - xOffset);
-            xEnd   = ceil  (xMid + xOffset) + 1.0;
-            for (x = xBegin; x < xEnd; ++x)
-            {
-                if (DistanceSquared (line, x, y) <= halfWidth2)
+
+        if (yBegin > yEnd)
+            SwapInt (&yBegin, &yEnd);
+        else
+            ++yEnd;
+        if (yBegin < 0) yBegin = 0;
+        if (yEnd > numRows) yEnd = numRows;
+
+        xMid = x1 + slope * (yBegin - y1);
+        xOffset = fabs (halfWidth * dr / dy);
+        
+        for (y = yBegin; y != yEnd; ++y)
+        {
+            xBegin = (int) floor (xMid - xOffset);
+            xEnd   = (int) ceil  (xMid + xOffset) + 1;
+            for (x = xBegin; x != xEnd; ++x)
+            {
+                if (x >=0 && x < numCols)
                 {
-                    if (x >=0 && y >= 0) {
-                        pixel = psAlloc (sizeof(PixelPos));
-                        pixel->x = (unsigned int) x;
-                        pixel->y = (unsigned int) y;
-                        psArrayAdd (pixels, 1024, pixel);
-                        psFree (pixel);
-                    }
+                    psImageSet(pixels, x, y, 1);
                 }
             }

branches/pap/magic/remove/src/Line.h

@@ -21 +21 @@
 extern bool LineClip (Line *line, int numCols, int numRows);
 
+/** Clip the line between (minX,minY) and (maxX,maxY)
+
+    @param[in,out] line line to be clipped within the bounds
+    @param[in] minX minimum X (columns) for the line
+    @param[in] minY minimum Y (rows) for the line
+    @param[in] maxX maximum X (columns) for the line
+    @param[in] maxY maximum Y (rows) for the line
+    @return true if line overlaps the clip boundaries           */
+
+extern bool LineClipFull (Line *line, int minX, int minY, int maxX, int maxY);
+
 /** Map a line to an image for its specified width and append as
     a list of pixel positions
@@ -26 +37 @@
     @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            */
 
-extern void PixelsFromLine (StreakPixels* pixels, Line *line);
+extern void PixelsFromLine (StreakPixels* pixels, Line *line,
+                            int numCols, int numRows);
 
 #endif /* STREAK_LINE_H */

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

@@ -28 +28 @@
     streaksrelease.o
 
-# STREAKSFLAGS=-DSTREAKS_COMPRESS_OUTPUT=1
+ISDESTREAKED_OBJECTS=      \
+    ${COMMON_OBJECTS} \
+    isdestreaked.o
+
+
+HEADERS= Line.h  streaksastrom.h  streaksextern.h  streaksio.h  streaksremove.h
+
+STREAKSFLAGS=-DSTREAKS_COMPRESS_OUTPUT=1
 OPTFLAGS= -g -O2
-OPTFLAGS= -g
+#OPTFLAGS= -g
 CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS} ${STREAKSFLAGS}
 #CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS}
 LDFLAGS=`psmodules-config --libs`
 
-PROGRAMS= streaksremove streaksreplace streakscompare streaksrelease
+PROGRAMS= streaksremove streaksreplace streakscompare streaksrelease isdestreaked
+HEADERS=Line.h streaksastrom.h streaksextern.h streaksio.h streaksremove.h
 
 all:    ${PROGRAMS}
 
+${REMOVE_OBJECTS}:      ${HEADERS}
 streaksremove:  ${REMOVE_OBJECTS}
 
@@ -47 +56 @@
 streaksrelease:  ${RELEASE_OBJECTS}
 
+isdestreaked:   ${ISDESTREAKED_OBJECTS}
+
+
 install:        ${PROGRAMS}
         install -t  $(PSCONFDIR)/$(PSCONFIG)/bin streaksremove
@@ -52 +64 @@
         install -t  $(PSCONFDIR)/$(PSCONFIG)/bin streakscompare
         install -t  $(PSCONFDIR)/$(PSCONFIG)/bin streaksrelease
+        install -t  $(PSCONFDIR)/$(PSCONFIG)/bin isdestreaked
 
 clean:

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

@@ -150 +150 @@
  
 bool
+SkyToLocal(strkPt *outPt, strkAstrom *astrom, double ra, double dec)
+{
+    // generate a local project using the RA, DEC of the 0,0 pixel of the chip as the
+    // projection center with the same plate scale as the nominal TP->Sky astrometry.
+
+    pmFPA *fpa   = (pmFPA *) astrom->fpa;
+    pmChip *chip = (pmChip *) astrom->chip;
+
+    // find the RA,DEC coords of the 0,0 pixel for this chip:
+
+    psPlane ptTP;
+    psSphere ptSky;
+
+    ptSky.r = ra;
+    ptSky.d = dec;
+    ptSky.rErr = 0.0;
+    ptSky.dErr = 0.0;
+
+    psProject(&ptTP, &ptSky, fpa->toSky);
+
+    outPt->x = ptTP.x;
+    outPt->y = ptTP.y;
+
+    return true;
+}
+
+bool
+LocalToSky(strkPt *outPt, strkAstrom *astrom, strkPt *inPt)
+{
+    // generate a local project using the RA, DEC of the 0,0 pixel of the chip as the
+    // projection center with the same plate scale as the nominal TP->Sky astrometry.
+
+    pmFPA *fpa   = (pmFPA *) astrom->fpa;
+    pmChip *chip = (pmChip *) astrom->chip;
+
+    // find the RA,DEC coords of the 0,0 pixel for this chip:
+
+    psPlane ptTP;
+    psSphere ptSky;
+
+    ptTP.x = inPt->x;
+    ptTP.y = inPt->y;
+    ptTP.xErr = 0.0;
+    ptTP.yErr = 0.0;
+
+    psDeproject(&ptSky, &ptTP, fpa->toSky);
+
+    outPt->x = ptSky.r;
+    outPt->y = ptSky.d;
+
+    return true;
+}
+
+bool
+componentBounds(int *minX, int *minY, int *maxX, int *maxY, strkAstrom *astrom, int numCols, int numRows)
+{
+    // find the bounds of the (padded) chip region in tangent-plane coordinates
+
+    pmFPA *fpa   = (pmFPA *) astrom->fpa;
+    pmChip *chip = (pmChip *) astrom->chip;
+
+    psPlane ptCH, ptFP, TPo, TPx, TPy;
+
+    // coordinate of the chip center:
+    ptCH.x = 0.5*numCols;
+    ptCH.y = 0.5*numRows;
+    psPlaneTransformApply(&ptFP, chip->toFPA, &ptCH);
+    psPlaneTransformApply(&TPo, fpa->toTPA,  &ptFP);
+
+    // coordinate of the chip center + dX/2:
+    ptCH.x = numCols;
+    ptCH.y = 0.5*numRows;
+    psPlaneTransformApply(&ptFP, chip->toFPA, &ptCH);
+    psPlaneTransformApply(&TPx, fpa->toTPA,  &ptFP);
+
+    // coordinate of the chip center + dY/2:
+    ptCH.x = 0.5*numCols;
+    ptCH.y = numRows;
+    psPlaneTransformApply(&ptFP, chip->toFPA, &ptCH);
+    psPlaneTransformApply(&TPy, fpa->toTPA,  &ptFP);
+
+    // half-lengths of the two sides in tangent-plane coords:
+    double xSize = hypot (TPx.x - TPo.x, TPx.y - TPo.y);
+    double ySize = hypot (TPy.x - TPo.x, TPy.y - TPo.y);
+    double radius = hypot (xSize, ySize);
+
+    // define the region encompassed by the radius with some padding:
+    *minX = TPo.x - 1.1*radius;
+    *minY = TPo.y - 1.1*radius;
+    *maxX = TPo.x + 1.1*radius;
+    *maxY = TPo.y + 1.1*radius;
+
+    return true;
+}
+
+bool
 skyToCell(strkPt *outPt, strkAstrom *astrom, double ra, double dec)
 {
@@ -281 +377 @@
 }
 
-void
+bool
 linearizeTransforms(strkAstrom *astrom)
 {
     if (!pmAstromLinearizeTransforms((pmFPA *) astrom->fpa, (pmChip *) astrom->chip)) {
-        streaksExit("linear fit to astrometry failed\n", PS_EXIT_UNKNOWN_ERROR);
-    }
-}
+        // streaksExit("linear fit to astrometry failed\n", PS_EXIT_UNKNOWN_ERROR);
+        psError(PS_ERR_UNKNOWN, false, "linear fit to astrometry failed ignoring");
+        return false;
+    }
+    return true;
+}

branches/pap/magic/remove/src/streaksastrom.h

@@ -30 +30 @@
 extern bool skyToCell(strkPt *, strkAstrom *astrom, double ra, double dec);
 extern bool cellToSky(strkPt *, strkAstrom *astrom, double x, double y);
-extern void linearizeTransforms(strkAstrom *astrom);
+extern bool linearizeTransforms(strkAstrom *astrom);
+
+extern bool SkyToLocal(strkPt *outPt, strkAstrom *astrom, double ra, double dec);
+extern bool LocalToSky(strkPt *outPt, strkAstrom *astrom, strkPt *inPt);
+extern bool componentBounds(int *minX, int *minY, int *maxX, int *maxY, strkAstrom *astrom, int numCols, int numRows);
 
 #endif // STREAKS_ASTROM_H

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

@@ -34 +34 @@
     int i;
     Line line;
-    StreakPixels *pixels = psArrayAllocEmpty (1024);
+    StreakPixels *pixels = psImageAlloc(numCols, numRows, PS_TYPE_U8);
+    psImageInit(pixels, 0);
     int streaksOnComponent = 0;
+
+    int minX, minY, maxX, maxY;
+
+    // find the chip dimensions in the tangent-plane coordinates (length of hypotenuse)
+    componentBounds (&minX, &minY, &maxX, &maxY, astrom, numCols, numRows);
+
     for (i = 0; i != streaks->size; ++i)
     {
@@ -41 +48 @@
 
         line.width = streaks->list[i].width;
-        if (skyToCell (&line.begin, astrom,
-                       streaks->list[i].ra1, streaks->list[i].dec1) &&
-            skyToCell (&line.end, astrom,
-                       streaks->list[i].ra2, streaks->list[i].dec2) &&
+
+        /* Use tangent plane coordinates to narrow down the ra,dec range of the line closer to
+         * the chip boundaries.  Use these new ra,dec positions to generate the line on the
+         * chip using the full non-linear astrometry */
+           
+        // project the ends of the line using a linear projection centered on the chip center:
+        Line full;
+        SkyToLocal (&full.begin, astrom, streaks->list[i].ra1, streaks->list[i].dec1);
+        SkyToLocal (&full.end,   astrom, streaks->list[i].ra2, streaks->list[i].dec2);
+
+        // clip the line to a square box with diameter = hypotenuse of the chip image centerd
+        // on the chip center in tangent-plane coordinates.  skip the rest of this streak if
+        // the line does not intersect this region
+        if (!LineClipFull (&full, minX, minY, maxX, maxY)) {
+            continue;
+        }
+
+        // convert the end points back into ra, dec pairs:
+        strkPt sky1, sky2;
+        LocalToSky (&sky1, astrom, &full.begin);
+        LocalToSky (&sky2, astrom, &full.end);
+
+        if (skyToCell (&line.begin, astrom, sky1.x, sky1.y) &&
+            skyToCell (&line.end,   astrom, sky2.x, sky2.y) &&
             LineClip (&line, numCols, numRows))
         {
-            PixelsFromLine (pixels, &line);
+            PixelsFromLine (pixels, &line, numCols, numRows);
             streaksOnComponent++;
         }

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

@@ -4 +4 @@
 /** psLib includes */
 
-#include "psArray.h"
+#include "psImage.h"
 
 /** streakremove includes */
@@ -36 +36 @@
 /** For now, use psArray to define a list of PixelPos pointers */
 
-typedef psArray StreakPixels;
+typedef psImage StreakPixels;
 
 /** Create a list of pixel positions from a Streaks pointer list

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

@@ -10 +10 @@
 static nebServer *ourNebServer = NULL;
 
+// Assumptions about the file structure of non-raw files
+// The 'image' for each file (image, mask weight) is contained in the first
+// extension. 
+
+
 // open the files required for streaks procesing.
 // if remove is true the calling program is streaksremove and the recovery files are outputs
 // if false the recovery files are inputs
-streakFiles *openFiles(pmConfig *config, bool remove)
+streakFiles *openFiles(pmConfig *config, bool remove, char *program_name)
 {
     bool status;
@@ -21 +26 @@
 
     sf->config = config;
+    sf->program_name = basename(program_name);
+
+    if (remove) {
+        // remember pointer so that streaksExit can delete temps
+        setStreakFiles(sf);
+    }
+
 
     // error checking is done by sFileOpen. If a file can't be opened we just exit
@@ -59 +71 @@
     // if we are passed a mask it is camera stage mask and we also
     // need to destreak the chip level mask file as well
+    // If it doesn't exist, we didn't have a camera mask
     if (remove && sf->inMask && (stage == IPP_STAGE_CHIP)) {
-        sf->inChMask = sFileOpen(config, stage, "INPUT.CHMASK", NULL, true);
-        inputBasename = basename(sf->inChMask->name);
-        sf->outChMask = sFileOpen(config, stage, "OUTPUT", inputBasename, true);
-        sf->recChMask = sFileOpen(config, stage, "RECOVERY", inputBasename, false);
+        sf->inChMask = sFileOpen(config, stage, "INPUT.CHMASK", NULL, false);
+        if (sf->inChMask) {
+            inputBasename = basename(sf->inChMask->name);
+            sf->outChMask = sFileOpen(config, stage, "OUTPUT", inputBasename, true);
+            sf->recChMask = sFileOpen(config, stage, "RECOVERY", inputBasename, false);
+        }
     }
 
@@ -74 +89 @@
         } else {
             sf->recWeight = sFileOpen(config, stage, "RECOVERY.WEIGHT", NULL, true);
+        }
+    }
+    if (remove) {
+        sf->inSources = sFileOpen(config, stage, "SOURCES", NULL, false);
+        if (sf->inSources) {
+            inputBasename = basename(sf->inSources->name);
+            sf->outSources = sFileOpen(config, stage, "OUTPUT", inputBasename, true);
         }
     }
@@ -159 +181 @@
 }
 
+// figure out if a nebulous instance is a non-destreaked file
+static bool
+nebFileIsDestreaked(sFile *sfile)
+{
+    if (!sfile->resolved_name) {
+        psError(PS_ERR_PROGRAMMING, true, "resolved name is null");
+        return false;
+    }
+    psFits *fits = psFitsOpen(sfile->resolved_name, "r");
+    if (!fits) {
+        psError(PS_ERR_IO, true, "failed open %s", sfile->name);
+        // can't tell if it is a destreaked file
+        return false;
+    }
+    psMetadata *header = psFitsReadHeader(NULL, fits);
+    if (!header) {
+        psError(PS_ERR_IO, true, "failed to read header for: %s", sfile->name);
+        return false;
+    }
+    bool mdok;
+    bool isDestreaked = psMetadataLookupBool(&mdok, header, "PSDESTRK");
+    if (mdok && isDestreaked) {
+        return true;
+    } else {
+        psError(PS_ERR_IO, false, "output file already exists and may not be de-streaked: %s", sfile->name);
+        return false;
+    }
+}
+
 static psString
 resolveFilename(pmConfig *config, sFile *sfile, bool create)
@@ -170 +221 @@
             // delete the existing file, since there may be more than one
             // instance. It will get created below in pmConfigConvertFilename
-            if (nebFind(server, sfile->name)) {
+            if ((sfile->resolved_name = nebFind(server, sfile->name)) != NULL) {
+                if (!nebFileIsDestreaked(sfile)) {
+                    psError(PS_ERR_IO, false, "attempting to delete file that has not been destreaked %s", sfile->name);
+                    return NULL;
+                }
+                nebFree(sfile->resolved_name);
+                sfile->resolved_name = NULL;
                 nebDelete(server, sfile->name);
             }
@@ -192 +249 @@
     // all of the keywords in the raw image files written to the output destreaked files
 
-    if (!CHIP_LEVEL_INPUT(stage) && !strcmp(fileSelect, "INPUT")) {
+    if (!outputFilename && !CHIP_LEVEL_INPUT(stage) && !strcmp(fileSelect, "INPUT")) {
         // stage is warp or diff AND fileSelect eq "INPUT"
         // get data from pmFPAfile.
@@ -234 +291 @@
         sfile->header = (psMetadata*) psMemIncrRefCounter(sfile->pmfile->fpa->hdu->header);
 
+        sfile->nHDU = psFitsGetSize(sfile->pmfile->fits);
+
         return sfile;
     }
@@ -249 +308 @@
     }
 
-    // if outputFilename is not null name it contains the "directory"
+    // if outputFilename is not null name it contains the "directory" (perhaps with a prefix like SR_)
     // and outputFilename is the basename name of the file (or nebulous key)
     // and the file is to be opened for writing
@@ -333 +392 @@
 
 void
+addDestreakKeyword(psMetadata *header)
+{
+    psMetadataAddBool(header, PS_LIST_TAIL, "PSDESTRK", PS_META_REPLACE,
+        "Have streaks been removed from image?", true);
+}
+
+void
+addRecoveryKeyword(psMetadata *header)
+{
+    psMetadataAddBool(header, PS_LIST_TAIL, "PSRECOVR", PS_META_REPLACE,
+        "Does this image contain excised streak pixels?", true);
+}
+
+void
 copyPHU(streakFiles *sfiles, bool remove)
 {
@@ -343 +416 @@
         streaksExit("", PS_EXIT_DATA_ERROR);
     }
-
-    // TODO: add keyword indicating that streaks have been removed
+    psMetadata *recHeader = NULL;
+    if (remove && sfiles->recImage) {
+       recHeader = psMetadataCopy(NULL, imageHeader);
+       addRecoveryKeyword(recHeader);
+    }
+
+    // add keyword indicating that streaks have been removed
+    addDestreakKeyword(imageHeader);
+
     if (!psFitsWriteBlank(sfiles->outImage->fits, imageHeader, NULL)) {
         psError(PS_ERR_IO, false, "failed to write primary header to %s",
@@ -350 +430 @@
         streaksExit("", PS_EXIT_DATA_ERROR);
     }
-    // TODO: add keyword indicating that this is the recovery image
-    if (remove && sfiles->recImage && !psFitsWriteBlank(sfiles->recImage->fits, imageHeader, NULL)) {
+    if (recHeader && !psFitsWriteBlank(sfiles->recImage->fits, recHeader, NULL)) {
         psError(PS_ERR_IO, false, "failed to write primary header to %s",
             sfiles->recImage->resolved_name);
         streaksExit("", PS_EXIT_DATA_ERROR);
     }
+    psFree(recHeader);
+    recHeader = NULL;
     psFree(imageHeader);
 
@@ -366 +447 @@
             streaksExit("", 1);
         }
-        // TODO: add keyword indicating that streaks have been removed
+        if (remove && sfiles->recMask) {
+            recHeader = psMetadataCopy(NULL, maskHeader);
+            // add keyword indicating that this is the recovery image
+           addRecoveryKeyword(recHeader);
+        }
+        // add keyword indicating that streaks have been removed
+        addDestreakKeyword(maskHeader);
         if (!psFitsWriteBlank(sfiles->outMask->fits, maskHeader, NULL)) {
             psError(PS_ERR_IO, false, "failed to write primary header to %s",
@@ -372 +459 @@
             streaksExit("", PS_EXIT_DATA_ERROR);
         }
-        // TODO: add keyword indicating that this is the recovery image
-        if (remove && sfiles->recMask && !psFitsWriteBlank(sfiles->recMask->fits, maskHeader, NULL)) {
+        if (recHeader && !psFitsWriteBlank(sfiles->recMask->fits, recHeader, NULL)) {
             psError(PS_ERR_IO, false, "failed to write primary header to %s",
                 sfiles->recMask->resolved_name);
             streaksExit("", PS_EXIT_DATA_ERROR);
         }
+        psFree(recHeader);
+        recHeader = NULL;
         psFree(maskHeader);
     }
@@ -388 +476 @@
             streaksExit("", 1);
         }
-        // TODO: add keyword indicating that streaks have been removed
+        if (remove && sfiles->recWeight) {
+            recHeader = psMetadataCopy(NULL, weightHeader);
+            // add keyword indicating that this is a recovery image
+           addRecoveryKeyword(recHeader);
+        }
+
+        // add keyword indicating that streaks have been removed
+        addDestreakKeyword(weightHeader);
         if (!psFitsWriteBlank(sfiles->outWeight->fits, weightHeader, NULL)) {
             psError(PS_ERR_IO, false, "failed to write primary header to %s",
@@ -394 +489 @@
             streaksExit("", PS_EXIT_DATA_ERROR);
         }
-        // TODO: add keyword indicating that this is a recovery image
-        if (remove && sfiles->recWeight && !psFitsWriteBlank(sfiles->recWeight->fits, weightHeader, NULL)) {
+        if (recHeader && !psFitsWriteBlank(sfiles->recWeight->fits, recHeader, NULL)) {
             psError(PS_ERR_IO, false, "failed to write primary header to %s",
                 sfiles->recWeight->resolved_name);
@@ -401 +495 @@
         }
         psFree(weightHeader);
+        psFree(recHeader);
     }
 }
@@ -530 +625 @@
         streaksExit("", PS_EXIT_DATA_ERROR);
     }
-
+ 
     bool status;
     in->extname = psMetadataLookupStr(&status, in->header, "EXTNAME");
@@ -565 +660 @@
 
 static void
-setFitsOptions(sFile *sfile, int bitpix, float bscale, float bzero)
+setFitsOptions(sFile *sfile, int bitpix, float bscale, float bzero, psFitsCompressionType compType,
+    psVector *tiles)
 {
     if (!sfile) {
@@ -578 +674 @@
     sfile->fits->options->bscale = bscale;
     sfile->fits->options->bzero = bzero;
-}
-
-void
-copyFitsOptions(sFile *out, sFile *rec, sFile *in)
-{
+
+    psFitsSetCompression(sfile->fits, compType, tiles, 8, 0, 0);
+}
+
+void
+copyFitsOptions(sFile *out, sFile *rec, sFile *in, psVector *tiles)
+{
+    bool mdok;
+    psString compTypeStr = psMetadataLookupStr(&mdok, in->header, "ZCMPTYPE");
+    psFitsCompressionType compType = psFitsCompressionTypeFromString(compTypeStr);
+    if (compType == PS_FITS_COMPRESS_NONE) {
+        return;
+    }
     // Get current BITPIX, BSCALE, BZERO, EXTNAME
     // Probably not necessary to look the numerical values up in this
@@ -609 +713 @@
 
 #ifdef STREAKS_COMPRESS_OUTPUT
-    // Paul says that I should be able to leave this blank
-    bitpix = 0;
-    setFitsOptions(out, bitpix, bscale, bzero);
-    setFitsOptions(rec, bitpix, bscale, bzero);
+    // printf("%d %f %f\n", bitpix, bscale, bzero);
+    setFitsOptions(out, bitpix, bscale, bzero, compType, tiles);
+    setFitsOptions(rec, bitpix, bscale, bzero, compType, tiles);
 #endif
 }
@@ -789 +892 @@
 
 bool
-replicate(sFile *sfile, void *xattr)
-{
-    if (!sfile->inNebulous) {
+replicate(sFile *outFile, sFile *inFile)
+{
+    if (!outFile->inNebulous) {
         return true;
     }
     nebServer *server = getNebServer(NULL);
 
-    // for now just set "user.copies" to 2
-    if (!nebSetXattr(server, sfile->name, "user.copies", "2",  NEB_REPLACE)) {
-        psError(PM_ERR_UNKNOWN, true, "nebSetXattr failed for %s\n%s", sfile->name, nebErr(server));
+    char *user_copies = nebGetXattr(server, inFile->name, "user.copies");
+    bool free_user_copies = true;
+    if (user_copies == NULL) {
+        user_copies = "2";
+        free_user_copies = false;
+    }
+    if (!nebSetXattr(server, outFile->name, "user.copies", user_copies,  NEB_REPLACE)) {
+        psError(PM_ERR_UNKNOWN, true, "nebSetXattr failed for %s\n%s", outFile->name, nebErr(server));
         return false;
     }
-    if (!nebReplicate(server, sfile->name, NULL, NULL)) {
-        psError(PM_ERR_UNKNOWN, true, "nebSetXattr failed for %s\n%s", sfile->name, nebErr(server));
+    if (!nebReplicate(server, outFile->name, "any", NULL)) {
+        psError(PM_ERR_UNKNOWN, true, "nebReplicate failed for %s\n%s", outFile->name, nebErr(server));
         return false;
+    }
+    if (free_user_copies) {
+        nebFree(user_copies);
     }
     return true;
@@ -816 +927 @@
     bool status = false;
 
-    // XXX: TODO: need a nebGetXatrr function, but there isn't one
-    // another option would be to take the number of copies to be
-    // created as an option. That way the system could decide
-    // whether to replicate anything other than raw Image files
-    void *xattr = NULL;
-
-    if (!replicate(sfiles->outImage, xattr)) {
+    if (!replicate(sfiles->outImage, sfiles->inImage)) {
         psError(PM_ERR_SYS, false, "failed to replicate outImage.");
         return false;
     }
 
-#ifdef notyet
-    // XXX: don't replicate mask and weight images until we can look up
-    // the input's xattr. There may be a perl program that can getXattr
     if (sfiles->outMask) {
-        // get xattr from input to see if we need to replicate
-        if (!replicate(sfiles->outMask, xattr)) {
+        if (!replicate(sfiles->outMask, sfiles->inMask)) {
             psError(PM_ERR_SYS, false, "failed to replicate outImage.");
             return false;
         }
     }
-    if (sfiles->outWeight) {
-        // get xattr from input to see if we need to replicate
-        if (!replicate(sfiles->outWeight, xattr)) {
+    if (sfiles->outChMask) {
+        if (!replicate(sfiles->outChMask, sfiles->inChMask)) {
             psError(PM_ERR_SYS, false, "failed to replicate outImage.");
             return false;
         }
     }
-#endif
-
-//      replicate the recovery images (if in nebulous)
+    if (sfiles->outWeight) {
+        if (!replicate(sfiles->outWeight, sfiles->inWeight)) {
+            psError(PM_ERR_SYS, false, "failed to replicate outImage.");
+            return false;
+        }
+    }
+
+    if (sfiles->outSources) {
+        if (!replicate(sfiles->outSources, sfiles->inSources)) {
+            psError(PM_ERR_SYS, false, "failed to replicate outSources.");
+            return false;
+        }
+    }
+
+//      XXX: replicate the recovery images (if in nebulous)
 //      perhaps whether we do that or not should be configurable.
 //      Sounds like we need a recipe
@@ -903 +1015 @@
     }
 
+    if (sfiles->outChMask) {
+        if (!swapOutputToInput(sfiles->inChMask, sfiles->outChMask)) {
+            psError(PM_ERR_SYS, false, "failed to swap instances for chip mask.");
+            return false;
+        }
+    }
+
+    if (sfiles->outSources) {
+        if (!swapOutputToInput(sfiles->inSources, sfiles->outSources)) {
+            psError(PM_ERR_SYS, false, "failed to swap instances for sources.");
+            return false;
+        }
+    }
+
     if (!swapOutputToInput(sfiles->inImage, sfiles->outImage)) {
         psError(PM_ERR_SYS, false, "failed to swap instances for Image.");
@@ -908 +1034 @@
     }
 
+
     return true;
 }
@@ -914 +1041 @@
 deleteFile(sFile *sfile)
 {
-#if 0
-    // XXX API for nebDelete has changed; need to fix this later
     if (sfile->inNebulous) {
         nebServer *server = getNebServer(NULL);
@@ -930 +1055 @@
         }
     }
-#endif
     return true;
 }
@@ -938 +1062 @@
 {
     if (sfiles->outMask) {
-        if (!deleteFile(sfiles->outMask)) {
-            psError(PM_ERR_SYS, false, "failed to delete Mask.");
-            return false;
-        }
+        deleteFile(sfiles->outMask);
+    }
+
+    if (sfiles->outChMask) {
+        deleteFile(sfiles->outChMask);
     }
 
     if (sfiles->outWeight) {
-        if (!deleteFile(sfiles->outWeight)) {
-            psError(PM_ERR_SYS, false, "failed to delete Weight.");
-            return false;
-        }
-    }
-
-    if (!deleteFile(sfiles->outImage)) {
-        psError(PM_ERR_SYS, false, "failed to delete Image.");
-        return false;
+        deleteFile(sfiles->outWeight);
+    }
+
+    if (sfiles->outImage) {
+        deleteFile(sfiles->outImage);
+    }
+
+    if (sfiles->outSources) {
+        deleteFile(sfiles->outSources);
     }
 
@@ -1018 +1143 @@
         }
         if (printCounts) {
-            printf("time to NAN mask pixels: %f\n", psTimerClear("NAN_MASKED"));
+            psLogMsg(sfiles->program_name, PS_LOG_INFO, "time to NAN mask pixels: %f\n", psTimerClear("NAN_MASKED"));
             int totalPixels = image->numRows * image->numCols;
-            printf("pixels:        %10ld\n", totalPixels);
-            printf("masked pixels: %10ld %4.2f%%\n", maskedPixels, 100. * maskedPixels / totalPixels);
-            printf("nand pixels:   %10ld %4.2f%%\n", nandPixels, 100. * nandPixels / totalPixels);
+            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);
+            psLogMsg(sfiles->program_name, PS_LOG_INFO, "nand pixels:   %10ld %4.2f%%\n", nandPixels, 100. * nandPixels / totalPixels);
             if (weight) {
-                printf("nand weights:  %10ld %4.2f%%\n", nandWeights, 100. * nandWeights / totalPixels);
+                psLogMsg(sfiles->program_name, PS_LOG_INFO, "nand weights:  %10ld %4.2f%%\n", nandWeights, 100. * nandWeights / totalPixels);
             }
         }
@@ -1050 +1175 @@
 strkGetMaskValues(streakFiles *sfiles, psU32 *maskStreak, psU32 *maskMask)
 {
-    if (sfiles->inMask->header) {
+    if (sfiles->inMask && sfiles->inMask->header) {
         if (!pmConfigMaskReadHeader(sfiles->config, sfiles->inMask->header)) {
             streaksExit("failed to read mask values from file", PS_EXIT_CONFIG_ERROR);
@@ -1080 +1205 @@
     *maskMask = ~convPoor;
 }
+
+static void
+doCopyExtraExtensions(streakFiles *sf, sFile *in, sFile *out)
+{
+    for (int extnum = sf->extnum; extnum < in->nHDU; extnum++) {
+        moveExt(in, extnum);
+        readImage(in, extnum, sf->stage, false);
+
+        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;
+            }
+            if (!psFitsWriteImage(out->fits, in->header, in->image, 0, in->extname)) {
+                psError(PM_ERR_UNKNOWN, false, "failed to write image for extension %d\n", extnum);
+                streaksExit("", PS_EXIT_UNKNOWN_ERROR);
+            }
+            if (out->fits->options) {
+                out->fits->options->bitpix = bitpix;
+            }
+            if (!psFitsCompressionApply(out->fits, compress)) {
+                psError(PM_ERR_UNKNOWN, false, "failed to rest compression image for extension %d\n", extnum);
+                streaksExit("", PS_EXIT_UNKNOWN_ERROR);
+            }
+#endif
+        } else {
+            copyTable(out, in, extnum);
+        }
+    }
+}
+
+void
+copyExtraExtensions(streakFiles *sf)
+{
+    // XXX: Bad assumption, the begining of the mask and weight files have the same structure as the image
+    // So we are currently at the image portion
+    sf->extnum = sf->inImage->nHDU;
+
+    if (sf->inWeight && sf->outWeight) {
+        doCopyExtraExtensions(sf, sf->inWeight, sf->outWeight);
+    }
+    if (sf->inMask && sf->outMask) {
+        doCopyExtraExtensions(sf, sf->inMask, sf->outMask);
+    }
+    if (sf->inChMask && sf->outChMask) {
+        doCopyExtraExtensions(sf, sf->inChMask, sf->outChMask);
+    }
+}

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

@@ -2 +2 @@
 #define STREAKS_IO_H 1
 
-streakFiles *openFiles(pmConfig *config, bool remove);
+streakFiles *openFiles(pmConfig *config, bool remove, char *program_name);
 
 sFile *sFileOpen(pmConfig *config, ippStage stage, psString fileSelect,
@@ -14 +14 @@
 void copyPHU(streakFiles *sfiles, bool remove);
 void copyTable(sFile *out, sFile *in, int extnum);
-void copyFitsOptions(sFile *out, sFile *rec, sFile *in);
+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 copyExtraExtensions(streakFiles *sfiles);
 
 void writeImage(sFile *sfile, psString extname, int extnum);
 void writeImageCube(sFile *sfile, psArray *imagecube, psString extname, int extnum);
-bool replicate(sFile *sfile, void *xattr);
+bool replicate(sFile *outFile, sFile *inFile);
 void readImageFrom_pmFile(streakFiles *sf);
+
+void addDestreakKeyword(psMetadata *);
+void addRecoveryKeyword(psMetadata *);
 
 bool streakFilesNextExtension(streakFiles *sf);

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

@@ -23 +23 @@
     }
 
-    psMetadata *masks = psMetadataLookupMetadata(&status, config->recipes, "MASKS");
-    if (!status) {
-        psError(PM_ERR_CONFIG, false, "failed to lookup MASKS in recipes\n");
-        return PS_EXIT_CONFIG_ERROR;
-    }
-    psU8 poorWarp = (double) psMetadataLookupU8(&status, masks, "POOR.WARP");
-    if (!status) {
-        psError(PM_ERR_CONFIG, false, "failed to lookup mask value for POOR.WARP in recipes\n");
-        return PS_EXIT_CONFIG_ERROR;
-    }
-    // we're setting pixels with any mask bits execpt POOR.WARP to NAN
-    psU8 maskMask = ~poorWarp;
+    // 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?
-    streakFiles *sfiles = openFiles(config, true);
+    streakFiles *sfiles = openFiles(config, true, argv[0]);
 
     if (sfiles->stage == IPP_STAGE_RAW) {
@@ -64 +55 @@
         }
 
+        // 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);
 
@@ -71 +67 @@
         printf("time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT"));
     } while (streakFilesNextExtension(sfiles));
+
+    // check the weight and mask files for extra extensions that might be in files
+    // (covariance matrix for example)
+    copyExtraExtensions(sfiles);
+
+
 
     psTimerStart("CLOSE_IMAGES");
@@ -77 +79 @@
     printf("time to close images: %f\n", psTimerClear("CLOSE_IMAGES"));
 
-#ifdef NOTYET
-    if (!replicateOutputs(sfiles)) {
-        psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
-        psErrorStackPrint(stderr, "");
-        exit(PS_EXIT_UNKNOWN_ERROR);
-    }
-
-    if (psMetadataLookupBool(&status, config->arguments, "REPLACE")) {
-        //     swap the instances for the input and output
-        //     Note this is a database operation. No file I/O is performed
-        if (!swapOutputsToInputs(sfiles)) {
-            psError(PS_ERR_UNKNOWN, false, "failed to swap files");
-
-            // XXX: Now what? I guess swapOutputsToInputs will need to undo anything that
-            // it has done and give a detailed report of what happened
-
-            psErrorStackPrint(stderr, "");
-            exit(PS_EXIT_UNKNOWN_ERROR);
-        }
-
-        if (psMetadataLookupBool(&status, config->arguments, "REMOVE")) {
-            //      delete the temporary storage objects (which now points to the original image(s)
-            if (!deleteTemps(sfiles)) {
-                psError(PS_ERR_UNKNOWN, false, "failed to delete temporary files");
-                // XXX: Now what? At this point the output files have been swapped, so we can't
-                // repeat the operation.
-
-                // Returning error status here is problematic. The inputs have been streak removed
-                // but they're still lying around
-                // Maybe just print an error message and
-                // let other system tools clean up
-                psErrorStackPrint(stderr, "");
-                exit(PS_EXIT_UNKNOWN_ERROR);
-            }
-        }
-    }
-#endif  // REPLACE, REMOVE
     printf("time to run streaksrelease: %f\n", psTimerClear("STREAKSREMOVE"));
 
@@ -275 +240 @@
 
     // set up the compression parameters
-#ifdef STREAKS_COMPRESS_OUTPUT
-    // compression of the image pixels is disabled for now. Some consortium members
-    // have problems reading them
-    copyFitsOptions(sf->outImage, sf->recImage, sf->inImage);
-
-    // XXX: TODO: can we derive these values from the input header?
-    // psFitsCompressionGet(sf->inImage->image) gives compression none
-    // perhaps we should just use the definition of COMP_IMG in the configuration 
-    psFitsSetCompression(sf->outImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-    if (sf->recImage) {
-        psFitsSetCompression(sf->recImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-    }
-#endif
+    copyFitsOptions(sf->outImage, sf->recImage, sf->inImage, sf->tiles);
 
     if (sf->inMask) {
@@ -306 +259 @@
             }
 
-#ifdef STREAKS_COMPRESS_OUTPUT
-            // XXX: see note above
-            copyFitsOptions(sf->outMask, sf->recMask, sf->inMask);
-            psFitsSetCompression(sf->outMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-            if (sf->recMask) {
-                psFitsSetCompression(sf->recMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-            }
+            copyFitsOptions(sf->outMask, sf->recMask, sf->inMask, sf->tiles);
             if (sf->outChMask) {
-                copyFitsOptions(sf->outChMask, sf->recChMask, sf->inMask);
-                psFitsSetCompression(sf->outChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-                if (sf->recChMask) {
-                    psFitsSetCompression(sf->recChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-                }
-            }
-#endif
+                copyFitsOptions(sf->outChMask, sf->recChMask, sf->inMask, sf->tiles);
+            }
         }
     }
@@ -332 +274 @@
         setupImageRefs(sf->outWeight, sf->recWeight, sf->inWeight, sf->extnum, exciseAll);
 
-#ifdef STREAKS_COMPRESS_OUTPUT
-        copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight);
-        // XXX: see note above
-        psFitsSetCompression(sf->outWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-        if (sf->recWeight) {
-            psFitsSetCompression(sf->recWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-        }
-#endif
+        copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight, sf->tiles);
     }
     // and for raw images, create sub images that represent the actual image

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

@@ +1 @@
+/*
+ * streaksremove
+ *
+ * Convert satellite streak detctions into masks and remove the covered pixels from the
+ * input images.
+ * Optionally swap the inputs with the outputs so that subsequent references to the original
+ * images access the destreaked versions.
+ */
+
 #include "streaksremove.h"
 
@@ -4 +13 @@
 static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll);
 static void exciseNonWarpedPixels(streakFiles *sfiles, double newMaskValue);
-static bool warpedPixel(streakFiles *sfiles, PixelPos *cellCoord);
+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 void writeImages(streakFiles *sf, bool exciseImageCube);
 static void updateAstrometry(streakFiles *sfiles);
+static void censorSources(streakFiles *sfiles, psU32 maskStreak);
 
 int
@@ -23 +33 @@
     }
 
-    psU32 maskStreak = 0;
-    psU32 maskMask = 0;
+    // Values to set for masked pixels
+    psU32 maskStreak = 0;           // for the image and weight (usually NAN, MAXINT for integer images)
+    psU32 maskMask = 0;             // value looked up for MASK.STREAK 
 
     psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS");
 
+    // call Paul Sydney's code to parse the streaks file that DetectStreaks produced
     Streaks *streaks = readStreaksFile(streaksFileName);
     if (!streaks) {
-        psErrorStackPrint(stderr, "failed to read streaks file: %s", streaksFileName);
+        psError(PS_ERR_UNKNOWN, "failed to read streaks file: %s", streaksFileName);
         streaksExit("", PS_EXIT_PROG_ERROR);
     }
 
-    streakFiles *sfiles = openFiles(config, true);
+    // open all of the input and output files, save their descriptions in the streakFiles struct
+    streakFiles *sfiles = openFiles(config, true, argv[0]);
     setupAstrometry(sfiles);
 
+    // Optionally we can set pixels that are masked to NAN since they couldn't have been
+    // examined for streaks. Usually this is done by the distribution system just prior
+    // to release
     bool nanForRelease = psMetadataLookupBool(&status, config->arguments, "NAN_FOR_RELEASE");
     if (nanForRelease && (sfiles->inMask == NULL)) {
@@ -52 +68 @@
 
     if (checkNonWarpedPixels ) {
-        // From ICD:
+        // From magic ICD:
         // In the raw and detrended images, the pixels which were not
         // included in any of the streak-processed warps must also be masked.
@@ -65 +81 @@
         }
         psF64 cwp_t = psTimerClear("COMPUTE_WARPED_PIXELS");
-        printf("time to compute warped pixels: %f\n", cwp_t);
+        psLogMsg("streaksremove", PS_LOG_INFO, "time to compute warped pixels: %f\n", cwp_t);
     }
     
     if (sfiles->stage == IPP_STAGE_RAW) {
-        // copy PHU to output files
+        // Except for raw stage, all of our (GPC1) files have one image extension.
+        // Raw files have a phu and multiple extensions, one per chip
+        // Since this is a raw file, copy it's PHU to output files
         copyPHU(sfiles, true);
 
@@ -82 +100 @@
     int totalStreakPixels = 0;
 
-    // Iterate through each component of the input (there is only one except for raw images)
+    // Iterate through each component of the input (except for raw images there is only one)
     do {
         bool exciseImageCube = false;
@@ -110 +128 @@
             psTimerStart("GET_STREAK_PIXELS");
 
-            StreakPixels *pixels = streak_on_component (streaks, sfiles->astrom,
-                                      sfiles->inImage->numCols, sfiles->inImage->numRows);
-
-            printf("time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS"));
-
+            // call Paul Sydney's code to compute the set of pixels that are covered by the detected streaks
+            StreakPixels *pixels = streak_on_component(streaks, sfiles->astrom, sfiles->inImage->numCols,
+                                                        sfiles->inImage->numRows);
+            psLogMsg("streaksremove", PS_LOG_INFO, "time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS"));
             
+            // if this extension contained an image, excise the streaked pixels.
+            // otherwise it contained an image cube (video cell) which is handled in the if block
             if (sfiles->inImage->image) {
                 if (checkNonWarpedPixels) {
@@ -124 +143 @@
                     exciseNonWarpedPixels(sfiles, maskStreak);
 
-                    printf("time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED"));
+                    psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED"));
                 }
-                totalStreakPixels +=  psArrayLength(pixels);
+
+
                 psTimerStart("REMOVE_STREAKS");
-                for (int i = 0; i < psArrayLength (pixels); ++i) {
-                    PixelPos *pixelPos = psArrayGet (pixels, i);
-
-                    if (!checkNonWarpedPixels || warpedPixel(sfiles, pixelPos)) {
-
-                        excisePixel(sfiles, pixelPos->x, pixelPos->y, true, maskStreak);
-
-                    } else {
-                        // This pixel was not included in any warp and has thus already excised
-                        // by exciseNonWarpedPixels
+
+                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
+                            }
+                        }
                     }
                 }
-                printf("time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS"));
+
+                psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS"));
 
                 if (nanForRelease) {
+                    // set any pixels that were masked, to NAN (unless they are already NAN)
                     setMaskedToNAN(sfiles, maskMask, true);
                 }
 
             } else { 
-                // this component contains an image cube, excise it completely
+                // this component contains an image cube
+                // For now excise it completely
                 exciseImageCube = true;
             }
-            psArrayElementsFree (pixels);
             psFree(pixels);
         }
 
         if (sfiles->stage == IPP_STAGE_CHIP) {
+            // as a convience to the user of the output, replace the bogus WCS transform in the
+            // chip processed files with the data calcuated by psastro at the camera stage
+            // (actually we use a linear approximation)
             updateAstrometry(sfiles);
         }
 
-        // write out the destreaked temporary images and the recovery images
+        censorSources(sfiles, maskStreak);
+
+        // write the destreaked "temporary" images and the recovery images
         writeImages(sfiles, exciseImageCube);
 
-        printf("time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT"));
+
+        psLogMsg("streaksremove", PS_LOG_INFO, "time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT"));
     } while (streakFilesNextExtension(sfiles));
 
+
     psFree(streaks);
 
-    printf("pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, 100. * totalStreakPixels / totalPixels);
+    psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, 100. * totalStreakPixels / totalPixels);
+
+    // check the weight and mask files for extra extensions that might be in files
+    // (covariance matrix for example)
+    copyExtraExtensions(sfiles);
+
+    // all done close the files. This is where the files are written so it can take a long time.
 
     psTimerStart("CLOSE_IMAGES");
-    // close all files
+
     closeImages(sfiles);
-    printf("time to close images: %f\n", psTimerClear("CLOSE_IMAGES"));
+
+    psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", psTimerClear("CLOSE_IMAGES"));
+
+
+    if (!replicateOutputs(sfiles)) {
+        psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
+        deleteTemps(sfiles);
+        psErrorStackPrint(stderr, "");
+        exit(PS_EXIT_UNKNOWN_ERROR);
+    }
 
     // NOTE: from here on we can't just quit if something goes wrong.
     // especially if we're working at the raw stage
-
-    if (!replicateOutputs(sfiles)) {
-        psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
-        psErrorStackPrint(stderr, "");
-        exit(PS_EXIT_UNKNOWN_ERROR);
-    }
+    // turn off automatic deletion of output files by streaksExit
+    setStreakFiles(NULL);
 
     if (psMetadataLookupBool(&status, config->arguments, "REPLACE")) {
         //     swap the instances for the input and output
-        //     Note this is a database operation. No file I/O is performed
+        //     Note this is a nebulous database operation. No file I/O is performed
         if (!swapOutputsToInputs(sfiles)) {
-            psError(PS_ERR_UNKNOWN, false, "failed to swap files");
-
-            // XXX: Now what? I guess swapOutputsToInputs will need to undo anything that
-            // it has done and give a detailed report of what happened
-
-            psErrorStackPrint(stderr, "");
+            // XXX: Now what?
+            // It is up to the program that reverts failed destreak runs to insure that
+            // any input files that have been swapped are restored and that the de-streaked
+            // versions are deleted
+
+            psErrorStackPrint(stderr, "failed to swap files");
+
+            // XXX: pick a specific error code for this failure
             exit(PS_EXIT_UNKNOWN_ERROR);
         }
-
-        if (psMetadataLookupBool(&status, config->arguments, "REMOVE")) {
-            //      delete the temporary storage objects (which now points to the original image(s)
-            if (!deleteTemps(sfiles)) {
-                psError(PS_ERR_UNKNOWN, false, "failed to delete temporary files");
-                // XXX: Now what? At this point the output files have been swapped, so we can't
-                // repeat the operation.
-
-                // Returning error status here is problematic. The inputs have been streak removed
-                // but they're still lying around
-                // Maybe just print an error message and
-                // let other system tools clean up
-                psErrorStackPrint(stderr, "");
-                exit(PS_EXIT_UNKNOWN_ERROR);
-            }
-        }
-    }
+    }
+    // all done. Clean up to look for memory leaks.
 
     psFree(sfiles);
@@ -218 +249 @@
     streaksNebulousCleanup(); 
     pmConfigDone();
-    printf("time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE"));
+    psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE"));
     psLibFinalize();
 
@@ -235 +266 @@
     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-remove: remove the original image after processing (requires -replace)\n");
     fprintf(stderr, "\t-keepnonwarped: 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");
@@ -311 +341 @@
     bool gotReplace = false;
     if ((argnum = psArgumentGet(argc, argv, "-replace"))) {
+        if (!nebulousImage) {
+            psError(PS_ERR_UNKNOWN, true, "replace is only supported for nebulous files");
+            usage();
+        }
         gotReplace = true;
         psArgumentRemove(argnum, &argc, argv);
@@ -394 +428 @@
         psArgumentRemove(argnum, &argc, argv);
     }
+    if ((argnum = psArgumentGet(argc, argv, "-sources"))) {
+        psArgumentRemove(argnum, &argc, argv);
+        bool nebulousSources = IN_NEBULOUS(argv[argnum]);
+        if (nebulousSources != nebulousImage) {
+            psError(PS_ERR_UNKNOWN, true, "sources file must have %snebulous path with %s image path\n",
+                nebulousImage ? "" : "non ", nebulousImage ? "nebulous" : "regular");
+            usage();
+        }
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "SOURCES", 0,
+                "name of input sources file", argv[argnum]);
+        psArgumentRemove(argnum, &argc, argv);
+    }
 
     if ((argnum = psArgumentGet(argc, argv, "-tmproot"))) {
@@ -403 +449 @@
             usage();
         }
-        psString dir = pathToDirectory(argv[argnum]);
-        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "directory for temporary files",
-            dir);
-        psFree(dir);
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "path for (temporary if replae) output files",
+            argv[argnum]);
         psArgumentRemove(argnum, &argc, argv);
     } else {
@@ -415 +459 @@
     if ((argnum = psArgumentGet(argc, argv, "-recovery"))) {
         psArgumentRemove(argnum, &argc, argv);
-        psString dir = pathToDirectory(argv[argnum]);
-        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "RECOVERY", 0, "directory for recovery files",
-            dir);
-        psFree(dir);
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "RECOVERY", 0, "path for recovery files",
+            argv[argnum]);
         psArgumentRemove(argnum, &argc, argv);
     } else if ((stage == IPP_STAGE_RAW) && gotReplace) {
         psError(PS_ERR_UNKNOWN, true, "-recovery is required for -stage raw with -replace\n");
         usage();
-    }
-
-    if ((argnum = psArgumentGet(argc, argv, "-remove"))) {
-        if (!gotReplace) {
-            psError(PS_ERR_UNKNOWN, true, "-replace is required with -remove\n");
-            usage();
-        }
-        psArgumentRemove(argnum, &argc, argv);
-        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "REMOVE", 0, "remove original files",
-            true);
     }
 
@@ -451 +483 @@
 updateAstrometry(streakFiles *sf)
 {
+    // XXX: why do I check this here? Shouldn't it be just around the call to linearizeTransforms?
     if (sf->bilevelAstrometry) {
 
-        linearizeTransforms(sf->astrom);
+        if (!linearizeTransforms(sf->astrom)) {
+            // fit failed, leave the astrometry unchanged
+            return;
+        }
 
         if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
@@ -501 +537 @@
         }
     }
-    sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header);
+    sf->outImage->header =  psMemIncrRefCounter(sf->inImage->header);
     if (sf->recImage) {
-        sf->recImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header);
-    }
+        sf->recImage->header = psMetadataCopy(NULL, sf->inImage->header);
+        addRecoveryKeyword(sf->recImage->header);
+    }
+    addDestreakKeyword(sf->outImage->header);
 
     if (!SFILE_IS_IMAGE(sf->inImage)) {
@@ -520 +558 @@
 
     // set up the compression parameters
-#ifdef STREAKS_COMPRESS_OUTPUT
-    // compression of the image pixels is disabled for now. Some consortium members
-    // have problems reading them
-    copyFitsOptions(sf->outImage, sf->recImage, sf->inImage);
-
-    // XXX: TODO: can we derive these values from the input header?
-    // psFitsCompressionGet(sf->inImage->image) gives compression none
-    // perhaps we should just use the definition of COMP_IMG in the configuration 
-    psFitsSetCompression(sf->outImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-    if (sf->recImage) {
-        psFitsSetCompression(sf->recImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-    }
-#endif
+    copyFitsOptions(sf->outImage, sf->recImage, sf->inImage, sf->tiles);
 
     if (sf->inMask) {
@@ -539 +565 @@
             sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
             if (sf->recMask) {
-                sf->recMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
-            }
+                sf->recMask->header = psMetadataCopy(NULL, sf->outMask->header);
+                addRecoveryKeyword(sf->recMask->header);
+            }
+            addDestreakKeyword(sf->outMask->header);
             if (updateAstrometry) {
                 pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
@@ -554 +582 @@
             }
 
-#ifdef STREAKS_COMPRESS_OUTPUT
-            // XXX: see note above
-            copyFitsOptions(sf->outMask, sf->recMask, sf->inMask);
-            psFitsSetCompression(sf->outMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-            if (sf->recMask) {
-                psFitsSetCompression(sf->recMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-            }
+            copyFitsOptions(sf->outMask, sf->recMask, sf->inMask, sf->tiles);
             if (sf->outChMask) {
-                copyFitsOptions(sf->outChMask, sf->recChMask, sf->inMask);
-                psFitsSetCompression(sf->outChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-                if (sf->recChMask) {
-                    psFitsSetCompression(sf->recChMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
-                }
-            }
-#endif
+                copyFitsOptions(sf->outChMask, sf->recChMask, sf->inMask, sf->tiles);
+            }
         }
     }
@@ -576 +593 @@
         sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);
         if (sf->recWeight) {
-            sf->recWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);
-        }
+            sf->recWeight->header = psMetadataCopy(NULL, sf->outWeight->header);
+            addRecoveryKeyword(sf->recWeight->header);
+        }
+        addDestreakKeyword(sf->outWeight->header);
         if (updateAstrometry) {
             pmAstromWriteWCS(sf->inWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
@@ -583 +602 @@
         setupImageRefs(sf->outWeight, sf->recWeight, sf->inWeight, sf->extnum, exciseAll);
 
-#ifdef STREAKS_COMPRESS_OUTPUT
-        copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight);
-        // XXX: see note above
-        psFitsSetCompression(sf->outWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-        if (sf->recWeight) {
-            psFitsSetCompression(sf->recWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-        }
-#endif
-    }
-    // and for raw images, create sub images that represent the actual image
-    // area (no overscan)
+        copyFitsOptions(sf->outWeight, sf->recWeight, sf->inWeight, sf->tiles);
+    }
 
     return true;
 }
-
-
 
 static void
@@ -746 +754 @@
 
 static bool
-warpedPixel(streakFiles *sfiles, PixelPos *cellCoord)
+warpedPixel(streakFiles *sfiles, int x, int y)
 {
     PixelPos chipCoord;
@@ -757 +765 @@
     // we clip so that the streak calculation code doesn't have to
     // clipping here insures that we don't touch the overscan regions
-    if ((cellCoord->x < 0) || (cellCoord->x >= sfiles->inImage->numCols) ||
-        (cellCoord->y < 0) || (cellCoord->y >= sfiles->inImage->numRows)) {
+    if ((x < 0) || (x >= sfiles->inImage->numCols) ||
+        (y < 0) || (y >= sfiles->inImage->numRows)) {
         return false;
     }
 
-    cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, cellCoord->x, cellCoord->y);
+    cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, x, y);
 
     if (chipCoord.x < 0 || chipCoord.x >= sfiles->warpedPixels->numCols) {
@@ -774 +782 @@
 }
 
+// read a sources file (.cmf) and remove any rows whose coordinate is convered by a
+// streak mask
+static void
+censorSources(streakFiles *sfiles, psU32 maskStreak)
+{
+    if ((!sfiles->inSources) || (!sfiles->outMask)) {
+        return;
+    }
+    psImage *maskImage = sfiles->outMask->image;
+    if (!maskImage) {
+        psError(PS_ERR_IO, false, "maskImage is null!");
+        streaksExit("", PS_EXIT_PROG_ERROR);
+    }
+
+    sFile *in = sfiles->inSources;
+    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);
+        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/pap/magic/remove/src/streaksremove.h

@@ -62 +62 @@
     sFile *outChMask;
     sFile *recChMask;
+    sFile *inSources;
+    sFile *outSources;
     psString class_id;
     pmFPAfile  *inAstrom;
@@ -72 +74 @@
     psVector    *tiles;
     double  transparentStreaks;
+    psString    program_name;
 } streakFiles;
 
@@ -87 +90 @@
 extern ippStage parseStage(psString);
 extern psString pathToDirectory(char *path);
+extern void setStreakFiles( streakFiles *);
 
 #define CHIP_LEVEL_INPUT(_stage) ((_stage == IPP_STAGE_RAW) || (_stage == IPP_STAGE_CHIP))
@@ -94 +98 @@
 #define IN_NEBULOUS(_filename) (!strncasecmp(_filename, "neb://", strlen("neb://")))
 
+
 #endif // STREAKS_H

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

@@ -33 +33 @@
     }
 
-    streakFiles *sfiles = openFiles(config, false);
+    streakFiles *sfiles = openFiles(config, false, argv[0]);
 
     if (sfiles->stage == IPP_STAGE_RAW) {
@@ -97 +97 @@
     }
 
-#ifdef NOTYET
-    if (psMetadataLookupBool(&status, config->arguments, "REPLACE")) {
-        //     swap the instances for the input and output
-        //     Note this is a database operation. No file I/O is performed
-        if (!swapOutputsToInputs(sfiles)) {
-            psError(PS_ERR_UNKNOWN, false, "failed to swap files");
-
-            // XXX: Now what? I guess swapOutputsToInputs will need to undo anything that
-            // it has done and give a detailed report of what happened
-
-            psErrorStackPrint(stderr, "");
-            exit(PS_EXIT_UNKNOWN_ERROR);
-        }
-
-        if (psMetadataLookupBool(&status, config->arguments, "REMOVE")) {
-            //      delete the temporary storage objects (which now points to the original image(s)
-            if (!deleteTemps(sfiles)) {
-                psError(PS_ERR_UNKNOWN, false, "failed to delete temporary files");
-                // XXX: Now what? At this point the output files have been swapped, so we can't
-                // repeat the operation.
-
-                // Returning error status here is problematic. The inputs have been streak removed
-                // but they're still lying around
-                // Maybe just print an error message and
-                // let other system tools clean up
-                psErrorStackPrint(stderr, "");
-                exit(PS_EXIT_UNKNOWN_ERROR);
-            }
-        }
-    }
-#endif  // REPLACE, REMOVE
 //    nebServerFree(ourNebServer);
     psFree(config);
@@ -344 +313 @@
 
     // set up the compression parameters
-#ifdef STREAKS_COMPRESS_OUTPUT
-    // compression of the image pixels is disabled for now. Some consortium members
-    // have problems reading compressed images.
-    copyFitsOptions(sf->outImage, sf->recImage, sf->inImage);
-
-    // XXX: TODO: can we derive these values from the input header?
-    // psFitsCompressionGet(sf->inImage->image) gives compression none
-    // perhaps we should just use the definition of COMP_IMG in the configuration 
-    psFitsSetCompression(sf->outImage->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
-#endif
+    copyFitsOptions(sf->outImage, sf->recImage, sf->inImage, sf->tiles);
+
     if (sf->inMask) {
         readImage(sf->inMask, sf->extnum, sf->stage, true);
@@ -362 +323 @@
         setupImageRefs(sf->outMask, NULL, sf->inMask, sf->extnum, false);
 
-        // XXX: see note above
-        copyFitsOptions(sf->outMask, NULL, sf->inMask);
-        psFitsSetCompression(sf->outMask->fits, PS_FITS_COMPRESS_PLIO, sf->tiles, 8, 0, 0);
+        copyFitsOptions(sf->outMask, NULL, sf->inMask, sf->tiles);
     }
 
@@ -374 +333 @@
         setupImageRefs(sf->outMask, NULL, sf->inMask, sf->extnum, false);
 
-        copyFitsOptions(sf->outWeight, NULL, sf->inWeight);
-        // XXX: see note above
-        psFitsSetCompression(sf->outWeight->fits, PS_FITS_COMPRESS_RICE, sf->tiles, 8, 0, 0);
+        copyFitsOptions(sf->outWeight, NULL, sf->inWeight, sf->tiles);
     }
 

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

@@ -33 +33 @@
 }
 
+streakFiles *ourStreakFiles = NULL;
+
+void
+setStreakFiles(streakFiles *sfiles)
+{
+    ourStreakFiles = sfiles;
+}
+
 // to enhance clarity in these programs we don't propagate errors up the stack
 // we just bail out
 void streaksExit(psString str, int exitCode) {
     psErrorStackPrint(stderr, str);
+    if (ourStreakFiles) {
+        deleteTemps(ourStreakFiles);
+    }
     exit(exitCode);
 }

@@ +1 @@
+isdestreaked
 streaksremove
 streaksreplace
 streakscompare
 streaksrelease
+makefile