Changeset 21156 for trunk/magic/remove/src

Timestamp:

Jan 23, 2009, 3:05:43 PM (17 years ago)

Author:

bills

Message:

Fixes to destreaking of raw images

Location:

trunk/magic/remove/src

Files:

• . (modified) (1 prop)
• .cvsignore (modified) (1 diff)
• Line.c (modified) (2 diffs)
• Makefile.simple (modified) (2 diffs)
• streaksastrom.c (modified) (10 diffs)
• streakscompare.c (modified) (1 diff)
• streaksio.c (modified) (12 diffs)
• streaksio.h (modified) (2 diffs)
• streaksrelease.c (modified) (6 diffs)
• streaksremove.c (modified) (18 diffs)
• streaksremove.h (modified) (2 diffs)
• streaksreplace.c (modified) (3 diffs)
• warpedpixels.c (modified) (2 diffs)

trunk/magic/remove/src

@@ -2 +2 @@
 streaksreplace
 streakscompare
+streaksrelease

@@ -2 +2 @@
 streaksreplace
 streakscompare
+streaksrelease

trunk/magic/remove/src/.cvsignore

@@ -2 +2 @@
 streaksreplace
 streakscompare
+streaksrelease

trunk/magic/remove/src/Line.c

@@ -321 +321 @@
                 if (DistanceSquared (line, x, y) <= halfWidth2)
                 {
-                    pixel = psAlloc (sizeof(PixelPos));
-                    pixel->x = (int) x;
-                    pixel->y = (int) y;
-                    psArrayAdd (pixels, 1024, pixel);
-                    psFree (pixel);
+                    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);
+                    }
                 }
             }
@@ -368 +370 @@
                 if (DistanceSquared (line, x, y) <= halfWidth2)
                 {
-                    pixel = psAlloc (sizeof(PixelPos));
-                    pixel->x = (int) x;
-                    pixel->y = (int) y;
-                    psArrayAdd (pixels, 1024, pixel);
-                    psFree (pixel);
+                    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);
+                    }
                 }
             }

trunk/magic/remove/src/Makefile.simple

@@ -6 +6 @@
 COMMON_OBJECTS = \
         streaksio.o \
-        streaksutil.o
+        streaksutil.o \
+        streaksastrom.o
 
 REMOVE_OBJECTS=    \
     ${COMMON_OBJECTS} \
     streaksremove.o \
-    streaksastrom.o \
     streaksextern.o \
     warpedpixels.o \
@@ -30 +30 @@
 STREAKSFLAGS=-DSTREAKS_COMPRESS_OUTPUT=1
 OPTFLAGS= -g
-#CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS} ${STREAKSFLAGS}
-CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS}
+CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS} ${STREAKSFLAGS}
+#CFLAGS=`psmodules-config --cflags` -std=gnu99 ${OPTFLAGS}
 LDFLAGS=`psmodules-config --libs`
 

trunk/magic/remove/src/streaksastrom.c

@@ -17 +17 @@
 
     if (fromCell) {
-        // The metadata is the cell concepts
+        // XXX: I think that this is not needed
+        // We are only asked to get the concepts from the cell when we're reading from a pmfile (skycell)
+        // In that case we (x0,y0) = (0,0) and the parirties are positive
+        // If the x0, y0 are not zero they are probably in the FITS unity based array format
+        // We should probably assert that these are zero
         cell_x0 =  psMetadataLookupS32(&mdok, md, "CELL.X0");
         if (!mdok) {
@@ -28 +32 @@
             return NULL;
         }
-
         xParityCell = psMetadataLookupS32(&mdok, md, "CELL.XPARITY");
         if (!mdok || (xParityCell != 1 && xParityCell != -1)) {
@@ -65 +68 @@
             return false;
         }
-
     } else {
         // no metadata regular cell
+        // used by warpedpixels
         cell_x0 = 0;
         cell_y0 = 0;
@@ -89 +92 @@
 }
 
+void
+chipToCell(double *xCell, double *yCell, strkAstrom *astrom, double xChip, double yChip)
+{
+    // convert from chip to cell
+    if (astrom->xParity > 0) {
+        *xCell = xChip - astrom->cell_x0;
+    } else {
+        *xCell = astrom->cell_x0 - 1 - xChip;
+    }
+    if (astrom->yParity > 0) {
+        *yCell = yChip - astrom->cell_y0;
+    } else {
+        *yCell = astrom->cell_y0 - 1 - yChip;
+    }
+}
+
+void
+cellToChip(double *xChip, double *yChip, strkAstrom *astrom, double xCell, double yCell)
+{
+    if (astrom->xParity > 0) {
+        *xChip = xCell + astrom->cell_x0;
+    } else {
+        *xChip = astrom->cell_x0 - 1 - xCell;
+    }
+    if (astrom->yParity > 0) {
+        *yChip = yCell + astrom->cell_y0;
+    } else {
+        *yChip =  astrom->cell_y0 - 1 - yCell;
+    }
+}
+
+void
+cellToChipInt(int *xChip, int *yChip, strkAstrom *astrom, int xCell, int yCell)
+{
+    if (astrom->xParity > 0) {
+        *xChip = xCell + astrom->cell_x0;
+    } else {
+        *xChip = astrom->cell_x0 - 1 - xCell;
+    }
+    if (astrom->yParity > 0) {
+        *yChip = yCell + astrom->cell_y0;
+    } else {
+        *yChip =  astrom->cell_y0 - 1 - yCell;
+    }
+}
+ 
 bool
 skyToCell(strkPt *outPt, strkAstrom *astrom, double ra, double dec)
@@ -107 +156 @@
     psPlaneTransformApply(pt->chip, chip->fromFPA, pt->FP);
 
-    // convert from chip to cell
-    outPt->x = (pt->chip->x - astrom->cell_x0) * astrom->xParity;
-    outPt->y = (pt->chip->y - astrom->cell_y0) * astrom->yParity;
+    chipToCell(&outPt->x, &outPt->y, astrom, pt->chip->x, pt->chip->y);
 
 //    printf("cell: %f %f chip: %f %f\n", outPt->x, outPt->y, pt->chip->x, pt->chip->y);
@@ -115 +162 @@
     return true;
 }
- 
+
 bool
 cellToSky(strkPt *outPt, strkAstrom *astrom, double x, double y)
@@ -123 +170 @@
     pmAstromObj *pt = (pmAstromObj *) astrom->pt;
     
-    // XXX: TODO: confirm that this cell to chip transformation is correct
-    pt->chip->x = (x * astrom->xParity) + astrom->cell_x0;
-    pt->chip->y = (y * astrom->yParity) + astrom->cell_y0;
+    cellToChip(&pt->chip->x, &pt->chip->y, astrom, x, y);
 
     pt->chip->xErr = 0;  // Is setting these errors to zero the right thing to do?
@@ -142 +187 @@
 }
 
-void
-cellToChip(PixelPos *chip, strkAstrom *astrom, PixelPos *cell)
-{
-    // TODO: do I need to worry about going from int to float and back again here?
-    chip->x = (cell->x * astrom->xParity) + astrom->cell_x0;
-    chip->y = (cell->y * astrom->yParity) + astrom->cell_y0;
-}
  
 static bool
@@ -166 +204 @@
     PS_ASSERT_PTR_NON_NULL(hdu->header, 1)
     if (sf->bilevelAstrometry) {
-        // Do we get here for GPC1 ? I don't necessarily have an fpa for the input image
         if (!pmAstromReadBilevelMosaic(sf->inAstrom->fpa, phu->header)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to read bilevel mosaic astrometry for input FPA.");
@@ -232 +269 @@
     }
 }
+
+void
+linearizeTransforms(strkAstrom *astrom)
+{
+    if (!pmAstromLinearizeTransforms((pmFPA *) astrom->fpa, (pmChip *) astrom->chip)) {
+        streaksExit("linear fit to astrometry failed\n", PS_EXIT_UNKNOWN_ERROR);
+    }
+}

trunk/magic/remove/src/streakscompare.c

@@ -39 +39 @@
             streaksExit("failed to advance to next extesion\n", PS_EXIT_DATA_ERROR);
         }
-        readImage(file1, 0, stage);
-        readImage(file2, 0, stage);
+        readImage(file1, 0, stage, false);
+        readImage(file2, 0, stage, false);
 
         psImage *image1 = file1->image;

trunk/magic/remove/src/streaksio.c

@@ -41 +41 @@
 
     sf->inMask = sFileOpen(config, stage, "INPUT.MASK", NULL, false);
-    if (sf->inMask) {
+    if (sf->inMask && (sf->stage != IPP_STAGE_RAW)) {
         inputBasename = basename(sf->inMask->name);
         sf->outMask = sFileOpen(config, stage, "OUTPUT", inputBasename, true);
@@ -295 +295 @@
 streakFilesNextExtension(streakFiles *sf)
 {
-    // unless stage is raw or chip when we get here we're done
+    // unless stage is raw when we get here we're done
     if (sf->stage != IPP_STAGE_RAW) {
         if (sf->view) {
@@ -315 +315 @@
     if (sf->extnum < sf->nHDU) {
         moveExt(sf->inImage, sf->extnum);
-        if (sf->inMask) {
-            moveExt(sf->inMask, sf->extnum);
-        }
-        if (sf->inWeight) {
-            moveExt(sf->inWeight, sf->extnum);
-        }
+        // No mask and weight images for raw stage
         return true;
     } else {
+        // we're all done
         return false;
     }
@@ -354 +350 @@
 
     sFile *inMask = sfiles->inMask;
-    if (inMask) {
+    if (inMask && sfiles->outMask) {
         psMetadata *maskHeader = psFitsReadHeader(NULL, inMask->fits);
         if (!maskHeader) {
@@ -469 +465 @@
 }
 
-void
-setDataExtent(ippStage stage, sFile *in)
-{
-    if (stage == IPP_STAGE_RAW) {
+static void
+setDataExtent(ippStage stage, sFile *in, bool rawImage)
+{
+    if (rawImage) {
         psString datasec = psMetadataLookupStr(NULL, in->header, "DATASEC");
         if (!datasec) {
@@ -480 +476 @@
         int xmin, xmax, ymin, ymax;
         sscanf(datasec, "[%d:%d,%d:%d]", &xmin, &xmax, &ymin, &ymax);
+
+        // I've inadvertantly introduced an implicit assumption that cell coordinates
+        // are the same as image coordinates. This is true for GPC1
+        // Test for this condition not being true
+        if (xmin != 1) {
+            psError(PS_ERR_IO, true, "can't handle datasec with minimum x != 1: %s", datasec);
+            streaksExit("", PS_EXIT_PROG_ERROR);
+        }
+        if (ymin != 1) {
+            psError(PS_ERR_IO, true, "can't handle datasec with minimum y != 1: %s", datasec);
+            streaksExit("", PS_EXIT_PROG_ERROR);
+        }
+
+        // These values are different than the size of the actual image which includes
+        // the bias area.
         in->numCols = xmax - xmin + 1;
         in->numRows = ymax - ymin + 1;
@@ -489 +500 @@
 
 void
-readImage(sFile *in, int extnum, ippStage stage)
+readImage(sFile *in, int extnum, ippStage stage, bool isMask)
 {
     psRegion region = {0, 0, 0, 0};
@@ -526 +537 @@
             streaksExit("", PS_EXIT_DATA_ERROR);
         }
+        if (in->image->type.type == PS_TYPE_U16) {
+            in->exciseValue = 65535;
+            psMetadataAddU16(in->header, PS_LIST_TAIL, "BLANK", 0, "", 65535);
+            psMetadataAddU16(in->header, PS_LIST_TAIL, "ZBLANK", 0, "", 65535);
+        } else {
+            in->exciseValue = NAN;
+        }
     }  else {
         if (stage != IPP_STAGE_RAW) {
@@ -539 +557 @@
         psImage *image = (psImage *) (in->imagecube->data[0]);
     }
-    setDataExtent(stage, in);
+    setDataExtent(stage, in, (stage == IPP_STAGE_RAW) && !isMask);
 }
 
@@ -676 +694 @@
 writeImageCube(sFile *sfile, psArray *imagecube, psString extname, int extnum)
 {
-    if (!imagecube) {
+    if (!sfile || !imagecube) {
         return;
     }
@@ -687 +705 @@
     sfile->header = NULL;
 }
-
-#ifdef notdef
-void
-writeImages(streakFiles *sf, bool exciseImageCube)
-{
-    psString extname = NULL;
-    if (sf->nHDU > 1) {
-        bool mdok;
-        extname = psMetadataLookupStr(&mdok, sf->inImage->header, "EXTNAME");
-    }
-    if (sf->inImage->numZPlanes == 0)  {
-        // note exciseing complete images is handled in readAndCopyToOutput
-        writeImage(sf->outImage, extname, sf->extnum);
-        writeImage(sf->recImage, extname, sf->extnum);
-    } else {
-        // we have an image cube
-        double initValue;
-        if (exciseImageCube) {
-            // copy the entire input image to the recovery image
-            writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum);
-            initValue = NAN;
-        } else {
-            // otherwise write it to the output 
-            writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
-            initValue = 0;
-        }
-
-        // borrow one of the images from the imagecube and set it to init value
-        psImage *image = psArrayGet (sf->inImage->imagecube, 0);
-        psMemIncrRefCounter(image);
-        psImageInit(image, initValue);
-        if (exciseImageCube) {
-            sf->outImage->image = image;
-            writeImage(sf->outImage, extname, sf->extnum);
-        } else {
-            // write zero valued image to reccovery
-            if (sf->recImage) {
-                sf->recImage->image = image;
-                writeImage(sf->recImage, extname, sf->extnum);
-            }
-        }
-        // Assumption: there are no mask and weight images for video cells
-        return;
-    }
-    if (sf->outMask) {
-        writeImage(sf->outMask, extname, sf->extnum);
-        writeImage(sf->recMask, extname, sf->extnum);
-    }
-    if (sf->outWeight) {
-        writeImage(sf->outWeight, extname, sf->extnum);
-        writeImage(sf->recWeight, extname, sf->extnum);
-    }
-}
-#endif
 
 void
@@ -952 +916 @@
 }
 
+bool
+isExciseValue(double value, double exciseValue)
+{
+    if (isnan(exciseValue)) {
+        return isnan(value);
+    } else {
+        return value == exciseValue;
+    }
+}
+
+void
+setMaskedToNAN(streakFiles *sfiles, psU8 maskMask, bool printCounts)
+{
+        int maskedPixels = 0;
+        int nandPixels = 0;
+        int nandWeights = 0;
+
+        psImage *image = sfiles->outImage->image;
+        psImage *mask = sfiles->inMask->image;
+        psImage *weight = NULL;
+        if (sfiles->outWeight) {
+            weight = sfiles->outWeight->image;
+        }
+        double exciseValue = sfiles->inImage->exciseValue;
+
+        if (printCounts) {
+            psTimerStart("NAN_MASKED");
+        }
+
+        for (int y=0; y < sfiles->inImage->numRows; y++) {
+            for (int x=0; x < sfiles->inImage->numCols; x++) {
+                // these gets are not necessary, we could just set the pixels to nan
+                // but I want to get the counts
+                double imageVal  = psImageGet(image, x, y);
+                psU8 maskVal;
+                if (sfiles->stage == IPP_STAGE_RAW) {
+                    int xChip, yChip;
+                    cellToChipInt(&xChip, &yChip, sfiles->astrom, x, y);
+                    maskVal = psImageGet(mask, xChip, yChip);
+                } else {
+                    maskVal = psImageGet(mask, x, y);
+                }
+                if (maskVal & maskMask) {
+                    ++maskedPixels;
+                    if (!isExciseValue(imageVal, sfiles->inImage->exciseValue)) {
+                        ++nandPixels;
+                        psImageSet(image, x, y, exciseValue);
+                    }
+                    if (weight) {
+                        double weightVal = weight ? psImageGet(weight, x, y) : 0;
+                        if (!isnan(weightVal)) {
+                            ++nandWeights;
+                            psImageSet(weight, x, y, NAN);
+                        }
+                    }
+                }
+            }
+        }
+        if (printCounts) {
+            printf("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);
+            if (weight) {
+                printf("nand weights:  %10ld %4.2f%%\n", nandWeights, 100. * nandWeights / totalPixels);
+            }
+        }
+}

trunk/magic/remove/src/streaksio.h

@@ -12 +12 @@
 void closeImage(sFile *sfile);
 
-void readImage(sFile *sfile, int extnum, ippStage stage);
+void readImage(sFile *sfile, int extnum, ippStage stage, bool isMask);
 void copyPHU(streakFiles *sfiles, bool remove);
 void copyTable(sFile *out, sFile *in, int extnum);
@@ -25 +25 @@
 bool streakFilesNextExtension(streakFiles *sf);
 
+bool isExciseValue(double, double);
+void setMaskedToNAN(streakFiles *sfiles, psU8 maskMask, bool printCounts);
 #endif // STREAKS_IO_H

trunk/magic/remove/src/streaksrelease.c

@@ -51 +51 @@
     }
 
-    int totalPixels = 0;
-    int totalMaskedPixels = 0;
-    int totalNANdPixels = 0;
-    int totalNANdWeights = 0;
-
     // Iterate through each component of the input (there is only one except for raw images)
     do {
@@ -70 +65 @@
         }
 
-        int maskedPixels = 0;
-        int nandPixels = 0;
-        int nandWeights = 0;
-
-        psImage *image = sfiles->outImage->image;
-        psImage *mask = sfiles->outMask->image;
-        psImage *weight= sfiles->outWeight->image;
-
-        for (int y=0; y < image->numRows; y++) {
-            for (int x=0; x < image->numCols; x++) {
-                // these gets are not necessary, we could just set the pixels to nan
-                // but I want to get the counts
-                double imageVal  = psImageGet(image, x, y);
-                double weightVal = psImageGet(weight, x, y);
-                psU8 maskVal   = psImageGet(mask, x, y);
-                if (maskVal & maskMask) {
-                    ++maskedPixels;
-                    if (!isnan(imageVal)) {
-                        ++nandPixels;
-                        psImageSet(image, x, y, NAN);
-                    }
-                    if (!isnan(weightVal)) {
-                        ++nandWeights;
-                        psImageSet(weight, x, y, NAN);
-                    }
-                }
-            }
-        }
-        totalPixels += sfiles->inImage->numRows * sfiles->inImage->numCols;
-        totalMaskedPixels += maskedPixels;
-        totalNANdPixels += nandPixels;
-        totalNANdWeights += nandWeights;
+        setMaskedToNAN(sfiles, maskMask, true);
 
         // write out the destreaked temporary images and the recovery images
@@ -108 +72 @@
         printf("time to process component %d: %f\n", sfiles->extnum, psTimerClear("PROCESS_COMPONENT"));
     } while (streakFilesNextExtension(sfiles));
-
-    printf("pixels:        %10ld\n", totalPixels);
-    printf("masked pixels: %10ld %4.2f%%\n", totalMaskedPixels, 100. * totalMaskedPixels / totalPixels);
-    printf("nand pixels:   %10ld %4.2f%%\n", totalNANdPixels, 100. * totalNANdPixels / totalPixels);
-    printf("nand weights:  %10ld %4.2f%%\n", totalNANdWeights, 100. * totalNANdWeights / totalPixels);
 
     psTimerStart("CLOSE_IMAGES");
@@ -292 +251 @@
     } else {
         // image data directly from psFits
-        readImage(sf->inImage, sf->extnum, sf->stage);
+        readImage(sf->inImage, sf->extnum, sf->stage, false);
     }
     sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header);
@@ -325 +284 @@
 
     if (sf->inMask) {
-        readImage(sf->inMask, sf->extnum, sf->stage);
+        readImage(sf->inMask, sf->extnum, sf->stage, true);
         sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
         if (sf->recMask) {
@@ -358 +317 @@
 
     if (sf->inWeight) {
-        readImage(sf->inWeight, sf->extnum, sf->stage);
+        readImage(sf->inWeight, sf->extnum, sf->stage, false);
         sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);
         if (sf->recWeight) {

trunk/magic/remove/src/streaksremove.c

@@ -9 +9 @@
 static void writeImages(streakFiles *sf, bool exciseImageCube);
 static bool replicateOutputs(streakFiles *sfiles);
+static void updateAstrometry(streakFiles *sfiles);
 
 int
@@ -21 +22 @@
     if (!config) {
         psError(PS_ERR_UNKNOWN, false, "failed to parse arguments\n");
-        return PS_EXIT_CONFIG_ERROR;
+        streaksExit("", PS_EXIT_CONFIG_ERROR);
     }
 
@@ -27 +28 @@
     if (!status) {
         psError(PM_ERR_CONFIG, false, "failed to lookup MASKS in recipes\n");
-        return PS_EXIT_CONFIG_ERROR;
+        streaksExit("", PS_EXIT_CONFIG_ERROR);
     }
     double maskStreak = (double) psMetadataLookupU8(&status, masks, "STREAK");
     if (!status) {
         psError(PM_ERR_CONFIG, false, "failed to lookup mask value for STREAK in recipes\n");
-        return PS_EXIT_CONFIG_ERROR;
-    }
+        streaksExit("", PS_EXIT_CONFIG_ERROR);
+    }
+
+    // optionally setting pixels with any mask bits execpt POOR.WARP to NAN
+    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");
+        streaksExit("", PS_EXIT_CONFIG_ERROR);
+    }
+    psU8 maskMask = ~poorWarp;
 
     psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS");
@@ -40 +49 @@
     if (!streaks) {
         psErrorStackPrint(stderr, "failed to read streaks file: %s", streaksFileName);
-        exit (PS_EXIT_PROG_ERROR);
+        streaksExit("", PS_EXIT_PROG_ERROR);
     }
 
     streakFiles *sfiles = openFiles(config, true);
     setupAstrometry(sfiles);
+
+    bool nanForRelease = psMetadataLookupBool(&status, config->arguments, "NAN_FOR_RELEASE");
+    if (nanForRelease && (sfiles->inMask == NULL)) {
+        psError(PS_ERR_UNKNOWN, true, "mask image must be supplied with -release.");
+        streaksExit("", PS_EXIT_PROG_ERROR);
+    }
 
     bool exciseAll = false;
@@ -113 +128 @@
             printf("time to get streak pixels: %f\n", psTimerClear("GET_STREAK_PIXELS"));
 
-            if (checkNonWarpedPixels) {
-                psTimerStart("EXCISE_NON_WARPED");
-
-                exciseNonWarpedPixels(sfiles, maskStreak);
-
-                printf("time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED"));
-            }
             
             if (sfiles->inImage->image) {
+                if (checkNonWarpedPixels) {
+                    psTimerStart("EXCISE_NON_WARPED");
+
+                    exciseNonWarpedPixels(sfiles, maskStreak);
+
+                    printf("time to excise non warped pixels: %f\n", psTimerClear("EXCISE_NON_WARPED"));
+                }
                 totalStreakPixels +=  psArrayLength(pixels);
                 psTimerStart("REMOVE_STREAKS");
@@ -137 +152 @@
                 }
                 printf("time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS"));
+
+                if (nanForRelease) {
+                    setMaskedToNAN(sfiles, maskMask, true);
+                }
+
             } else { 
                 // this component contains an image cube, excise it completely
@@ -142 +162 @@
             }
             psArrayElementsFree (pixels);
+        }
+
+        if (sfiles->stage == IPP_STAGE_CHIP) {
+        // if (CHIP_LEVEL_INPUT(sfiles->stage)) {
+            updateAstrometry(sfiles);
         }
 
@@ -302 +327 @@
             true);
     }
+
+    if ((argnum = psArgumentGet(argc, argv, "-release"))) {
+        psArgumentRemove(argnum, &argc, argv);
+        psMetadataAddBool(config->arguments, PS_LIST_TAIL, "NAN_FOR_RELEASE", 0, "set masked pixels to NAN",
+            true);
+    }
     
     if ((argnum = psArgumentGet(argc, argv, "-keepnonwarped"))) {
@@ -424 +455 @@
 
     return config;
+}
+
+static void
+updateAstrometry(streakFiles *sf)
+{
+    if (sf->bilevelAstrometry) {
+        linearizeTransforms(sf->astrom);
+        // XXX: pmAstrometry.c should set the correct CTYPE values
+        if (!pmAstromWriteWCS(sf->outImage->header, sf->inAstrom->fpa, sf->chip, 0.001)) {
+            psError(PS_ERR_UNKNOWN, false, "failed to update astrometry for extension %d", sf->extnum);
+            streaksExit("", PS_EXIT_UNKNOWN_ERROR);
+        }
+        if (sf->outMask) {
+            pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
+        }
+        if (sf->outWeight) {
+            pmAstromWriteWCS(sf->outWeight->header, sf->inAstrom->fpa, sf->chip, 0.001);
+        }
+    }
 }
 
@@ -438 +488 @@
     } else {
         // image data directly from psFits
-        readImage(sf->inImage, sf->extnum, sf->stage);
+        readImage(sf->inImage, sf->extnum, sf->stage, false);
         if (SFILE_IS_IMAGE(sf->inImage)) {
             sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false,
@@ -493 +543 @@
 
     if (sf->inMask) {
-        readImage(sf->inMask, sf->extnum, sf->stage);
-        sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
-        if (sf->recMask) {
-            sf->recMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
-        }
-        if (updateAstrometry) {
-            pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
-        }
-        setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll);
-        if (sf->outChMask) {
-            sf->outChMask->header = (psMetadata *) psMemIncrRefCounter(sf->outMask->header);
-            sf->outChMask->image = (psImage *) psMemIncrRefCounter(sf->outMask->image);
-            if (sf->recChMask) {
-                sf->recChMask->header = (psMetadata *) psMemIncrRefCounter(sf->recMask->header);
-                sf->recChMask->image = (psImage *) psMemIncrRefCounter(sf->recMask->image);
-            }
-        }
+        readImage(sf->inMask, sf->extnum, sf->stage, true);
+        if (sf->outMask) {
+            sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
+            if (sf->recMask) {
+                sf->recMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
+            }
+            if (updateAstrometry) {
+                pmAstromWriteWCS(sf->outMask->header, sf->inAstrom->fpa, sf->chip, 0.001);
+            }
+            setupImageRefs(sf->outMask, sf->recMask, sf->inMask, sf->extnum, exciseAll);
+            if (sf->outChMask) {
+                sf->outChMask->header = (psMetadata *) psMemIncrRefCounter(sf->outMask->header);
+                sf->outChMask->image = (psImage *) psMemIncrRefCounter(sf->outMask->image);
+                if (sf->recChMask) {
+                    sf->recChMask->header = (psMetadata *) psMemIncrRefCounter(sf->recMask->header);
+                    sf->recChMask->image = (psImage *) psMemIncrRefCounter(sf->recMask->image);
+                }
+            }
 
 #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);
-        }
-        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);
-            }
-        }
+            // 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);
+            }
+            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
+        }
     }
 
     if (sf->inWeight) {
-        readImage(sf->inWeight, sf->extnum, sf->stage);
+        readImage(sf->inWeight, sf->extnum, sf->stage, false);
         sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);
         if (sf->recWeight) {
@@ -661 +713 @@
 excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, double newMaskValue)
 {
+    double exciseValue = sfiles->inImage->exciseValue;
+
     // we clip so that the streak calculation code doesn't have to
     if ((x < 0) || (x >= sfiles->inImage->numCols) || (y < 0) || (y >= sfiles->inImage->numRows)) {
@@ -667 +721 @@
 
     double imageValue  = psImageGet (sfiles->inImage->image,  x, y);
-    if (sfiles->recImage && !isnan(imageValue) ) {
+    if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
         psImageSet (sfiles->recImage->image,  x, y, imageValue);
     }
 
     if (sfiles->transparentStreaks == 0) {
-        psImageSet (sfiles->outImage->image,  x, y, NAN);
+        psImageSet (sfiles->outImage->image,  x, y, exciseValue);
     } else {
         if (streak) {
@@ -678 +732 @@
             psImageSet (sfiles->outImage->image,  x, y, imageValue + sfiles->transparentStreaks);
         } else {
-            psImageSet (sfiles->outImage->image,  x, y, NAN);
-        }
-    }
-
-    if (sfiles->inWeight) {
+            psImageSet (sfiles->outImage->image,  x, y, exciseValue);
+        }
+    }
+
+    if (sfiles->outWeight) {
         if (sfiles->recWeight) {
             double weightValue = psImageGet (sfiles->inWeight->image, x, y);
             psImageSet (sfiles->recWeight->image, x, y, weightValue);
         }
+        // Assume that weight images are always a floating point type
         psImageSet (sfiles->outWeight->image, x, y, NAN);
     }
-    if (sfiles->inMask) {
+    if (sfiles->outMask) {
         if (sfiles->recMask) {
             double maskValue   = psImageGet (sfiles->inMask->image,   x, y);
@@ -705 +760 @@
     int xParity = sfiles->astrom->xParity;
     int yParity = sfiles->astrom->yParity;
-    int numCols = sfiles->inImage->numCols;
-    int numRows = sfiles->inImage->numRows;
+    int numCols = sfiles->inImage->numCols; // for raw images this was calculated from the width of datasec
+    int numRows = sfiles->inImage->numRows; // for raw images this was calculated from the height of datasec
 
 //    printf("%2d x0: %4d y0: %4d xpar: %d ypar: %d\n", sfiles->extnum, cell_x0, cell_y0, xParity, yParity);
 
+    // XXX: in the case of raw images we are making an implicit assumption here
+    // that the DATASEC begins at (x, y) = (0, 0) in the array for this image
+    // While this is true for GPC1 it isn't necessarily true
     for (int yCell = 0; yCell < numRows; yCell++) {
         int yChip;
@@ -751 +809 @@
 
     // 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)) {
@@ -756 +815 @@
     }
 
-    cellToChip(&chipCoord, sfiles->astrom, cellCoord);
-
-#ifdef notdef
-    // Shouldn't need this but for raw images I'm getting 
-    // x = warpedPixels->numCols for cell xy70 when cell.x = 0
-    // Perhaps cell_x0 should be IMNPIX-1 since IMNPIX is in the 1 based system of the fits
-    // headers
+    cellToChipInt(&chipCoord.x, &chipCoord.y, sfiles->astrom, cellCoord->x, cellCoord->y);
+
     if (chipCoord.x < 0 || chipCoord.x >= sfiles->warpedPixels->numCols) {
         return false;
@@ -769 +823 @@
         return false;
     }
-#endif
 
     return psImageGet(sfiles->warpedPixels, chipCoord.x, chipCoord.y) ? true : false;

trunk/magic/remove/src/streaksremove.h

@@ -26 +26 @@
     psArray     *imagecube;
     pmFPAfile   *pmfile;
-    int         numCols;
-    int         numRows;
+    int         col0;
+    int         numCols;    // number of columns in data section of image (doesn't include bias section)
+    int         xParity;
+    int         row0;
+    int         numRows;    // number of rows in data section of image (does not include bias section)
+    int         yParity;
+    double      exciseValue;
 } sFile;
 
@@ -75 +80 @@
     int numCols, int numRows);
 // can't declare this in streaksastrom due to header file ordering
-extern void cellToChip(PixelPos *chip, strkAstrom *astrom, PixelPos *cell);
+extern void cellToChip(double *xChip, double *yChip, strkAstrom *astrom, double xCell, double yCell);
+extern void cellToChipInt(int *xChip, int *yChip, strkAstrom *astrom, int xCell, int yCell);
 
 extern bool computeWarpedPixels(streakFiles *sf);

trunk/magic/remove/src/streaksreplace.c

@@ -327 +327 @@
     } else {
         // image data directly from psFits
-        readImage(sf->inImage, sf->extnum, sf->stage);
+        readImage(sf->inImage, sf->extnum, sf->stage, false);
     }
     // read the recovery pixels
-    readImage(sf->recImage, sf->extnum, sf->stage);
+    readImage(sf->recImage, sf->extnum, sf->stage, false);
 
     sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header);
@@ -359 +359 @@
 #endif
     if (sf->inMask) {
-        readImage(sf->inMask, sf->extnum, sf->stage);
-        readImage(sf->recMask, sf->extnum, sf->stage);
+        readImage(sf->inMask, sf->extnum, sf->stage, true);
+        readImage(sf->recMask, sf->extnum, sf->stage, true);
         sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
 
@@ -372 +372 @@
 
     if (sf->inWeight) {
-        readImage(sf->inWeight, sf->extnum, sf->stage);
-        readImage(sf->recWeight, sf->extnum, sf->stage);
+        readImage(sf->inWeight, sf->extnum, sf->stage, false);
+        readImage(sf->recWeight, sf->extnum, sf->stage, false);
 
         sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);

trunk/magic/remove/src/warpedpixels.c

@@ -41 +41 @@
         // write out the warped pixels
         psMetadata * header = psMetadataAlloc();
-        pmAstromWriteWCS(header, sf->inAstrom->fpa, sf->chip, 0.001);
+        if (sf->inAstrom->fpa->toSky->type != PS_PROJ_DIS) {
+            pmAstromWriteWCS(header, sf->inAstrom->fpa, sf->chip, 0.001);
+        }
         psFits *fits = psFitsOpen("warpedpixels.fits", "w");
 
@@ -76 +78 @@
     sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false, NULL, 
         sf->warpedPixels->numCols, sf->warpedPixels->numRows);
-    
 
     // convert corners of skycell to sky coordinates

@@ -2 +2 @@
 streaksreplace
 streakscompare
+streaksrelease