Index: trunk/magic/remove/src/streaksio.c
===================================================================
--- trunk/magic/remove/src/streaksio.c	(revision 21085)
+++ trunk/magic/remove/src/streaksio.c	(revision 21156)
@@ -41,5 +41,5 @@
 
     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,5 +295,5 @@
 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,12 +315,8 @@
     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,5 +350,5 @@
 
     sFile *inMask = sfiles->inMask;
-    if (inMask) {
+    if (inMask && sfiles->outMask) {
         psMetadata *maskHeader = psFitsReadHeader(NULL, inMask->fits);
         if (!maskHeader) {
@@ -469,8 +465,8 @@
 }
 
-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,4 +476,19 @@
         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,5 +500,5 @@
 
 void
-readImage(sFile *in, int extnum, ippStage stage)
+readImage(sFile *in, int extnum, ippStage stage, bool isMask)
 {
     psRegion region = {0, 0, 0, 0};
@@ -526,4 +537,11 @@
             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,5 +557,5 @@
         psImage *image = (psImage *) (in->imagecube->data[0]);
     }
-    setDataExtent(stage, in);
+    setDataExtent(stage, in, (stage == IPP_STAGE_RAW) && !isMask);
 }
 
@@ -676,5 +694,5 @@
 writeImageCube(sFile *sfile, psArray *imagecube, psString extname, int extnum)
 {
-    if (!imagecube) {
+    if (!sfile || !imagecube) {
         return;
     }
@@ -687,58 +705,4 @@
     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,2 +916,71 @@
 }
 
+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);
+            }
+        }
+}