Changeset 24557 for branches/eam_branches/20090522/magic/remove/src

Timestamp:

Jun 25, 2009, 2:00:56 PM (17 years ago)

Author:

eugene

Message:

merging changes from head

Location:

branches/eam_branches/20090522

Files:

• . (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) (4 diffs)
• magic/remove/src/isdestreaked.c (copied) (copied from trunk/magic/remove/src/isdestreaked.c )
• magic/remove/src/streaksastrom.c (modified) (1 diff)
• magic/remove/src/streaksastrom.h (modified) (1 diff)
• magic/remove/src/streaksextern.c (modified) (2 diffs)
• magic/remove/src/streaksio.c (modified) (20 diffs)
• magic/remove/src/streaksio.h (modified) (2 diffs)
• magic/remove/src/streaksrelease.c (modified) (6 diffs)
• magic/remove/src/streaksremove.c (modified) (19 diffs)
• magic/remove/src/streaksremove.h (modified) (2 diffs)
• magic/remove/src/streaksreplace.c (modified) (4 diffs)
• magic/remove/src/streaksutil.c (modified) (1 diff)

branches/eam_branches/20090522/magic/remove/src

@@ -3 +3 @@
 streakscompare
 streaksrelease
+makefile

@@ -3 +3 @@
 streakscompare
 streaksrelease
+makefile

branches/eam_branches/20090522/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 (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 (DistanceSquared (line, x, y) <= halfWidth2)
+                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);
-                    }
+                    pixel = psAlloc (sizeof(PixelPos));
+                    pixel->x = (unsigned int) x;
+                    pixel->y = (unsigned int) y;
+                    psArrayAdd (pixels, 1024, pixel);
+                    psFree (pixel);
                 }
             }
@@ -354 +475 @@
         
         // 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 (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 (DistanceSquared (line, x, y) <= halfWidth2)
+                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);
-                    }
+                    pixel = psAlloc (sizeof(PixelPos));
+                    pixel->x = (unsigned int) x;
+                    pixel->y = (unsigned int) y;
+                    psArrayAdd (pixels, 1024, pixel);
+                    psFree (pixel);
                 }
             }

branches/eam_branches/20090522/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/eam_branches/20090522/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
@@ -35 +42 @@
 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/eam_branches/20090522/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)
 {

branches/eam_branches/20090522/magic/remove/src/streaksastrom.h

@@ -32 +32 @@
 extern void 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/eam_branches/20090522/magic/remove/src/streaksextern.c

@@ -36 +36 @@
     StreakPixels *pixels = psArrayAllocEmpty (1024);
     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 +47 @@
 
         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/eam_branches/20090522/magic/remove/src/streaksio.c

@@ -19 +19 @@
     psMemSetDeallocator(sf, (psFreeFunc) streakFilesFree);
     memset(sf, 0, sizeof(*sf));
+
+    if (remove) {
+        // remember pointer so that streaksExit can delete temps
+        setStreakFiles(sf);
+    }
 
     sf->config = config;
@@ -160 +165 @@
 }
 
+// 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)
@@ -171 +205 @@
             // 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);
             }
@@ -193 +233 @@
     // 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.
@@ -250 +290 @@
     }
 
-    // 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
@@ -334 +374 @@
 
 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)
 {
@@ -344 +398 @@
         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",
@@ -351 +412 @@
         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);
 
@@ -367 +429 @@
             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",
@@ -373 +441 @@
             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);
     }
@@ -389 +458 @@
             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",
@@ -395 +471 @@
             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);
@@ -402 +477 @@
         }
         psFree(weightHeader);
+        psFree(recHeader);
     }
 }
@@ -566 +642 @@
 
 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) {
@@ -579 +656 @@
     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
@@ -610 +695 @@
 
 #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
 }
@@ -790 +874 @@
 
 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;
@@ -817 +909 @@
     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;
+        }
+    }
+
+//      XXX: replicate the recovery images (if in nebulous)
 //      perhaps whether we do that or not should be configurable.
 //      Sounds like we need a recipe
@@ -904 +990 @@
     }
 
+    if (sfiles->outChMask) {
+        if (!swapOutputToInput(sfiles->inChMask, sfiles->outChMask)) {
+            psError(PM_ERR_SYS, false, "failed to swap instances for chip mask.");
+            return false;
+        }
+    }
+
     if (!swapOutputToInput(sfiles->inImage, sfiles->outImage)) {
         psError(PM_ERR_SYS, false, "failed to swap instances for Image.");
@@ -936 +1029 @@
 {
     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);
     }
 

branches/eam_branches/20090522/magic/remove/src/streaksio.h

@@ -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);
@@ -20 +20 @@
 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/eam_branches/20090522/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?
@@ -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);
 
@@ -77 +73 @@
     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 +234 @@
 
     // 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 +253 @@
             }
 
-#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 +268 @@
         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/eam_branches/20090522/magic/remove/src/streaksremove.c

@@ +1 @@
+/*
+ * streaksremove
+ *
+ * Convert satellite streak detctions into masks and remove the covered pixels from the
+ * input images.
+ * Optionally swap the inputs with the outputs so that subsequent references to the original
+ * images access the destreaked versions.
+ */
+
 #include "streaksremove.h"
 
@@ -23 +32 @@
     }
 
-    psU32 maskStreak = 0;
-    psU32 maskMask = 0;
+    // Values to set for masked pixels
+    psU32 maskStreak = 0;           // for the image and weight (usually NAN, MAXINT for integer images)
+    psU32 maskMask = 0;             // value looked up for MASK.STREAK 
 
     psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS");
 
+    // call Paul Sydney's code to parse the streaks file that DetectStreaks produced
     Streaks *streaks = readStreaksFile(streaksFileName);
     if (!streaks) {
-        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);
     }
 
+    // open all of the input and output files, save their descriptions in the streakFiles struct
     streakFiles *sfiles = openFiles(config, true, argv[0]);
     setupAstrometry(sfiles);
 
+    // Optionally we can set pixels that are masked to NAN since they couldn't have been
+    // examined for streaks. Usually this is done by the distribution system just prior
+    // to release
     bool nanForRelease = psMetadataLookupBool(&status, config->arguments, "NAN_FOR_RELEASE");
     if (nanForRelease && (sfiles->inMask == NULL)) {
@@ -52 +67 @@
 
     if (checkNonWarpedPixels ) {
-        // From ICD:
+        // From magic ICD:
         // In the raw and detrended images, the pixels which were not
         // included in any of the streak-processed warps must also be masked.
@@ -69 +84 @@
     
     if (sfiles->stage == IPP_STAGE_RAW) {
-        // copy PHU to output files
+        // Except for raw stage, all of our (GPC1) files have one image extension.
+        // Raw files have a phu and multiple extensions, one per chip
+        // Since this is a raw file, copy it's PHU to output files
         copyPHU(sfiles, true);
 
@@ -82 +99 @@
     int totalStreakPixels = 0;
 
-    // Iterate through each component of the input (there is only one except for raw images)
+    // Iterate through each component of the input (except for raw images there is only one)
     do {
         bool exciseImageCube = false;
@@ -110 +127 @@
             psTimerStart("GET_STREAK_PIXELS");
 
-            StreakPixels *pixels = streak_on_component (streaks, sfiles->astrom,
-                                      sfiles->inImage->numCols, sfiles->inImage->numRows);
-
+            // call Paul Sydney's code to compute the set of pixels that are covered by the detected streaks
+            StreakPixels *pixels = streak_on_component(streaks, sfiles->astrom, sfiles->inImage->numCols,
+                                                        sfiles->inImage->numRows);
             psLogMsg("streaksremove", PS_LOG_INFO, "time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS"));
-
             
+            // if this extension contained an image, excise the streaked pixels.
+            // otherwise it contained an image cube (video cell) which is handled in the if block
             if (sfiles->inImage->image) {
                 if (checkNonWarpedPixels) {
@@ -126 +144 @@
                     psLogMsg("streaksremove", PS_LOG_INFO, "time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED"));
                 }
+
                 totalStreakPixels +=  psArrayLength(pixels);
+
                 psTimerStart("REMOVE_STREAKS");
+
+                // for each pixel covered by the streak
                 for (int i = 0; i < psArrayLength (pixels); ++i) {
                     PixelPos *pixelPos = psArrayGet (pixels, i);
 
+                    // if this pixel was not part of the warp, skip because it has already been
+                    // excised (unless we weren't asked to check)
                     if (!checkNonWarpedPixels || warpedPixel(sfiles, pixelPos)) {
 
@@ -140 +164 @@
                     }
                 }
+
                 psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS"));
 
                 if (nanForRelease) {
+                    // set any pixels that were masked, to NAN (unless they are already NAN)
                     setMaskedToNAN(sfiles, maskMask, true);
                 }
 
             } else { 
-                // this component contains an image cube, excise it completely
+                // this component contains an image cube
+                // For now excise it completely
                 exciseImageCube = true;
             }
@@ -155 +182 @@
 
         if (sfiles->stage == IPP_STAGE_CHIP) {
+            // as a convience to the user of the output, replace the bogus WCS transform in the
+            // chip processed files with the data calcuated by psastro at the camera stage
+            // (actually we use a linear approximation)
             updateAstrometry(sfiles);
         }
 
-        // write out the destreaked temporary images and the recovery images
+        // write the destreaked "temporary" images and the recovery images
         writeImages(sfiles, exciseImageCube);
 
@@ -168 +198 @@
     psLogMsg("streaksremove", PS_LOG_INFO, "pixels: %ld streak pixels: %ld %4.2f%%\n", totalPixels, totalStreakPixels, 100. * totalStreakPixels / totalPixels);
 
+    // all done close the files. This is where the files are written so it can take a long time.
+
     psTimerStart("CLOSE_IMAGES");
-    // close all files
+
     closeImages(sfiles);
+
     psLogMsg("streaksremove", PS_LOG_INFO, "time to close images: %f\n", psTimerClear("CLOSE_IMAGES"));
+
+
+    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
+
+            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);
@@ -403 +426 @@
             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 +436 @@
     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 +460 @@
 updateAstrometry(streakFiles *sf)
 {
+    // XXX: why do I check this here? Shouldn't it be just around the call to linearizeTransforms?
     if (sf->bilevelAstrometry) {
 
@@ -501 +511 @@
         }
     }
-    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 +532 @@
 
     // 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 +539 @@
             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 +556 @@
             }
 
-#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 +567 @@
         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 +576 @@
         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

branches/eam_branches/20090522/magic/remove/src/streaksremove.h

@@ -88 +88 @@
 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))
@@ -95 +96 @@
 #define IN_NEBULOUS(_filename) (!strncasecmp(_filename, "neb://", strlen("neb://")))
 
+
 #endif // STREAKS_H

branches/eam_branches/20090522/magic/remove/src/streaksreplace.c

@@ -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/eam_branches/20090522/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);
 }

@@ -3 +3 @@
 streakscompare
 streaksrelease
+makefile