Changeset 20573 for trunk/magic/remove/src/streaksremove.c

Timestamp:

Nov 6, 2008, 6:42:41 PM (18 years ago)

Author:

bills

Message:

completed excising of non-wapred pixels and various other changes.

File:

• trunk/magic/remove/src/streaksremove.c (modified) (26 diffs)

trunk/magic/remove/src/streaksremove.c

@@ -1 @@
-#include "streaks.h"
-#include "streaksextern.h"
+#include "streaksremove.h"
 #include "libgen.h"
 #include "unistd.h"
 
+// XXX: I don't think we need the lock and unlock functions
 extern bool  sFileLock(sFile * sfile);
 extern bool  sFileUnlock(sFile * sfile);
@@ -10 +10 @@
 static nebServer *ourNebServer = NULL;
 
-// for clarity in this program we don't propagate errors up the stack
+// to enhance clarity we don't propagate errors up the stack
 // we just bail out
 void streaksremoveExit(psString str, int exitCode) {
@@ -70 +70 @@
 }
 
-psString
+static psString
 resolveFilename(pmConfig *config, sFile *sfile, bool create)
 {
@@ -166 +166 @@
         sfile->resolved_name = resolveFilename(config, sfile, true);
         if (!sfile->resolved_name) {
-            psError(PS_ERR_IO, false, "Failed to create file %s", sfile->name);
+            psError(PS_ERR_IO, false, "Failed to resolve filename for %s", sfile->name);
             sFileFree(sfile);
             streaksremoveExit("", 1);
         }
         sfile->fits = psFitsOpen(sfile->resolved_name, "w");
-        sfile->fits->options = psFitsOptionsAlloc();
+        if (sfile->fits) {
+            sfile->fits->options = psFitsOptionsAlloc();
+        }
     } else {
         sfile->name = psStringCopy(name);
@@ -197 +199 @@
 
 
-bool
-astrometry_read(streakFiles *sf)
+static bool
+readAstrometry(streakFiles *sf)
 {
     pmHDU *phu = pmFPAviewThisPHU(sf->view, sf->inAstrom->fpa);
     if (phu) {
-        char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1");
+        bool status;
+        char *ctype = psMetadataLookupStr(&status, phu->header, "CTYPE1");
         if (ctype) {
             sf->bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
@@ -229 +232 @@
 }
 
+// load the fpa containing the astrometry, find the appropriate chip and process the data
 static void
 setupAstromFromFPA(streakFiles *sf)
@@ -259 +263 @@
     }
 
-    // read in the astrometry
-    astrometry_read(sf);
+    readAstrometry(sf);
 }
 
@@ -312 +315 @@
     }
     setupAstromFromFPA(sf);
-    if (CHIP_LEVEL_INPUT(sf->stage)) {
-        computeWarpedPixels(sf);
-    }
      
     psElemType tileType;                // Type corresponding to "long"
@@ -373 +373 @@
 {
     // unless stage is raw or chip when we get here we're done
-    if (!CHIP_LEVEL_INPUT(sf->stage)) {
+    if (sf->stage != IPP_STAGE_RAW) {
         sf->view->readout = -1;
         pmFPAfileIOChecks(sf->config, sf->view, PM_FPA_AFTER);
@@ -442 +442 @@
     }
 
-    // bool status;
     int argnum;
     ippStage stage = IPP_STAGE_NONE;
@@ -667 +666 @@
 
 // Determine whether the current header for this file is an image or not
+// Find a cleaner way to do this
 static bool
 isImage(sFile *in)
@@ -686 +686 @@
             }
         }
+    } else if (psMetadataLookupBool(&status, in->header, "ZIMAGE")) {
+        return true;
     } else if (in->nHDU == 1) {
         // no extensions in the file, can just return true? For now require 
@@ -721 +723 @@
 }
 
+void
+setDataExtent(ippStage stage, sFile *in)
+{
+    if (stage == IPP_STAGE_RAW) {
+        psString datasec = psMetadataLookupStr(NULL, in->header, "DATASEC");
+        if (!datasec) {
+            psError(PS_ERR_IO, false, "failed to find DATASEC in header");
+            streaksremoveExit("", PS_EXIT_DATA_ERROR);
+        }
+        int xmin, xmax, ymin, ymax;
+        sscanf(datasec, "[%d:%d,%d:%d]", &xmin, &xmax, &ymin, &ymax);
+        in->numCols = xmax - xmin + 1;
+        in->numRows = ymax - ymin + 1;
+    } else {
+        in->numCols = in->image->numCols;
+        in->numRows = in->image->numRows;
+    }
+}
+
 static void
-readImage(sFile *in, int extnum)
+readImage(sFile *in, int extnum, ippStage stage)
 {
     psRegion region = {0, 0, 0, 0};
@@ -756 +777 @@
             streaksremoveExit("", PS_EXIT_DATA_ERROR);
         }
-        in->numCols = in->image->numCols;
         in->numRows = in->image->numRows;
     }  else {
+        if (stage != IPP_STAGE_RAW) {
+            psError(PS_ERR_PROGRAMMING, true, "unexpected image cube found in non-raw image");
+            streaksremoveExit("", PS_EXIT_PROG_ERROR);
+        }
         in->imagecube = psFitsReadImageCube(in->fits, region);
         if (!in->imagecube) {
@@ -766 +790 @@
         }
         psImage *image = (psImage *) (in->imagecube->data[0]);
-        in->numCols = image->numCols;
-        in->numRows = image->numRows;
-    }
-
+    }
+    setDataExtent(stage, in);
 }
 
 static void
-setFitsOptions(sFile *sfile, psString extname, int bitpix, float bscale, float bzero)
+setFitsOptions(sFile *sfile, int bitpix, float bscale, float bzero)
 {
     if (sfile->fits->options) {
@@ -812 +834 @@
         bzero = psMetadataItemParseF32(bzeroItem);
     }
-    bool mdok;
-
-    psString extname = psMetadataLookupStr(&mdok, in->header, "EXTNAME");
-    setFitsOptions(out, extname, bitpix, bscale, bzero);
-    setFitsOptions(rec, extname, bitpix, bscale, bzero);
+
+//    setFitsOptions(out, bitpix, bscale, bzero);
+//    setFitsOptions(rec, bitpix, bscale, bzero);
 }
 
@@ -869 +889 @@
         // image data from pmFPAfile (diff or warp file)
         readImageFrom_pmFile(sf);
-        sf->astrom = streakSetAstrometry(sf->astrom, sf->inAstrom->fpa, sf->chip, true, sf->cell->concepts,
-            sf->inImage->numCols, sf->inImage->numRows);
+        sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa,
+                sf->chip, true, sf->cell->concepts, sf->inImage->numCols, sf->inImage->numRows);
     } else {
         // image data directly from psFits
-        readImage(sf->inImage, sf->extnum);
+        readImage(sf->inImage, sf->extnum, sf->stage);
         if (SFILE_IS_IMAGE(sf->inImage)) {
-            sf->astrom = streakSetAstrometry(sf->astrom, sf->inAstrom->fpa, sf->chip, false,
+            sf->astrom = streakSetAstrometry(sf->astrom, sf->stage, sf->inAstrom->fpa, sf->chip, false,
                 sf->inImage->header, sf->inImage->numCols, sf->inImage->numRows);
         }
+    }
+    if (!sf->astrom) {
+        psError(PS_ERR_UNKNOWN, false, "failed to set up astrometry for extnsion %d", sf->extnum);
+        streaksremoveExit("", PS_EXIT_UNKNOWN_ERROR);
     }
     sf->outImage->header = (psMetadata*) psMemIncrRefCounter(sf->inImage->header);
@@ -906 +930 @@
 
     if (sf->inMask) {
-        readImage(sf->inMask, sf->extnum);
+        readImage(sf->inMask, sf->extnum, sf->stage);
         sf->outMask->header = (psMetadata*) psMemIncrRefCounter(sf->inMask->header);
         sf->outMask->image = (psImage*) psMemIncrRefCounter(sf->inMask->image);
@@ -920 +944 @@
 
     if (sf->inWeight) {
-        readImage(sf->inWeight, sf->extnum);
+        readImage(sf->inWeight, sf->extnum, sf->stage);
         sf->outWeight->header = (psMetadata*) psMemIncrRefCounter(sf->inWeight->header);
         sf->outWeight->image = (psImage*) psMemIncrRefCounter(sf->inWeight->image);
@@ -943 +967 @@
 writeImage(sFile *sfile, psString extname, int extnum)
 {
+    if (!sfile->image) {
+        return;
+    }
     if (!psFitsWriteImage(sfile->fits, sfile->header, sfile->image, 0, extname)) {
         psError(PS_ERR_IO, false, "failed to write image to %s extnum: %d", 
@@ -953 +980 @@
     sfile->header = NULL;
 }
+
 static void
 writeImageCube(sFile *sfile, psArray *imagecube, psString extname, int extnum)
@@ -966 +994 @@
 
 static void
-writeImages(streakFiles *sf)
+writeImages(streakFiles *sf, bool exciseAll)
 {
     psString extname = NULL;
     if (sf->nHDU > 1) {
-        extname = psMetadataLookupStr(NULL, sf->inImage->header, "EXTNAME");
+        bool mdok;
+        extname = psMetadataLookupStr(&mdok, sf->inImage->header, "EXTNAME");
     }
     if (sf->inImage->numZPlanes == 0)  {
@@ -976 +1005 @@
         writeImage(sf->recImage, extname, sf->extnum);
     } else {
-        // XXX: TODO: if a streak touched this cell then
-        // write the image cube to the recovery image if no streaks write it
-        // to the output image
-        writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
+        double initValue;
+        if (exciseAll) {
+            writeImageCube(sf->recImage, sf->inImage->imagecube, extname, sf->extnum);
+            initValue = NAN;
+        } else {
+            writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
+            initValue = 0;
+        }
+
+        // borrow one of the images from the imagecube and set it to NAN
+        psImage *image = psArrayGet (sf->inImage->imagecube, 0);
+        psMemIncrRefCounter(image);
+        psImageInit(image, initValue);
+        if (exciseAll) {
+            sf->outImage->image = image;
+            writeImage(sf->outImage, extname, sf->extnum);
+        } else {
+            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) {
@@ -1194 +1241 @@
 }
 
+static void
+excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, double newMaskValue)
+{
+    double imageValue  = psImageGet (sfiles->inImage->image,  x, y);
+    psImageSet (sfiles->recImage->image,  x, y, imageValue);
+    psImageSet (sfiles->outImage->image,  x, y, NAN);
+
+    // TODO:
+    // Need to get mask weight for PM_MASK_DETECTOR.  Is this stored
+    // in a config somewhere?  Not sure how to properly set the maskValue.
+ 
+    if (sfiles->inWeight) {
+        double weightValue = psImageGet (sfiles->inWeight->image, x, y);
+        psImageSet (sfiles->recWeight->image, x, y, weightValue);
+        psImageSet (sfiles->outWeight->image, x, y, NAN);
+    }
+    if (sfiles->inMask) {
+        double maskValue   = psImageGet (sfiles->inMask->image,   x, y);
+        psImageSet (sfiles->recMask->image,   x, y, maskValue);
+        psImageSet (sfiles->outMask->image,   x, y, newMaskValue);
+    }
+}
+
+static void
+exciseNonWarpedPixels(streakFiles *sfiles, double newMaskValue)
+{
+    int cell_x0 = sfiles->astrom->cell_x0;
+    int cell_y0 = sfiles->astrom->cell_y0;
+    int xParity = sfiles->astrom->xParity;
+    int yParity = sfiles->astrom->yParity;
+    int numCols = sfiles->inImage->numCols;
+    int numRows = sfiles->inImage->numRows;
+
+//    printf("%2d x0: %4d y0: %4d xpar: %d ypar: %d\n", sfiles->extnum, cell_x0, cell_y0, xParity, yParity);
+
+    for (int yCell = 0; yCell < numRows; yCell++) {
+        int yChip;
+        if (yParity == 1) {
+            yChip = cell_y0 + yCell;
+        } else {
+            yChip = cell_y0 - yCell;
+        }
+
+        psU8 *pixels = sfiles->warpedPixels->data.U8[yChip];
+
+        if (xParity == 1) {
+            pixels += cell_x0;
+            for (int xCell = 0; xCell < numCols; xCell++, pixels++) {
+                if (! *pixels ) {
+                    excisePixel(sfiles, xCell, yCell, newMaskValue);
+                }
+            }
+        } else {
+            // negative parity
+            // point to the right most pixel in chip coordinates (lowest in cell coords)
+            pixels += cell_x0 - numCols;
+            for (int xCell = numCols - 1; xCell >= 0 ; xCell--, pixels++) {
+                if (!*pixels) {
+                    excisePixel(sfiles, xCell, yCell, newMaskValue);
+                }
+            }
+        }
+    }
+}
+
+bool
+warpedPixel(streakFiles *sfiles, PixelPos *cellCoord)
+{
+    PixelPos chipCoord;
+
+    if (!CHIP_LEVEL_INPUT(sfiles->stage)) {
+        // if we're here on a skycell image by definition this pixel was warped
+        return true;
+    }
+
+    cellToChip(&chipCoord, sfiles->astrom, cellCoord);
+
+    return psImageGet(sfiles->warpedPixels, chipCoord.x, chipCoord.y) ? true : false;
+}
+
 int
 main(int argc, char *argv[])
@@ -1220 +1347 @@
     streakFiles *sfiles = openFiles(config);
 
-    if (sfiles->stage == IPP_STAGE_RAW) {
-        // copy PHU to output files
-        copyPHU(sfiles);
-        // advance to the first image extension
-        if (!streakFilesNextExtension(sfiles)) {
-            psErrorStackPrint(stderr, "failed to advance to first extension of input images");
-            exit (PS_EXIT_PROG_ERROR);
-        }
-    }
-
-    // Iterate through components of image
-    do {
-
-        // read the images and copy the pixels image from the inputs to the outputs
-        // This also sets up the astrometry and initializes the recovery images
-        if (!readAndCopyToOutput(sfiles)) {
-            // this extension is not an image skip (video descriptor? 
-            // XXX is there anything else that could be in an input file that we need to handle?
-            continue;
-        }
-
-        // Identify pixels to mask because of streaks
-
-        pixels = streak_on_component (streaks, sfiles->astrom,
-                                      sfiles->inImage->numCols, sfiles->inImage->numRows);
-
-#ifdef notyet
-        // XXX: How is this formatted?
-        // PFS: Aren't these chips that are never warped?  They should be cleared
-        //      entirely or copy them to archive location.  I don't know how
-        //      to find these chips.  The fact that a warp is not created for
-        //      a chip must be recorded somewhere.  I don't like this approach
-        //      of including everything (ALL) and then removing them if they
-        //      are warped.  I don't know if it is practical.
+    if (CHIP_LEVEL_INPUT(sfiles->stage)) {
         // From ICD:
         // In the raw and detrended images, the pixels which were not
@@ -1260 +1354 @@
         // than a small fraction of the pixels are lit by the input image.
 
-        // Identify pixels to mask because not in a warp
-        warp_pixels = pixel_list_initialise(ALL); // List of pixels to be excised because not in a warp
-        foreach warp (warps) {
-            // Calculate warp region on image
-            image_warp = warp_on_component(warp, astrom, numCols, numRows);
-            foreach pixel specified by image_warp {
-                remove_pixel(warp_pixels, pixel);
+        computeWarpedPixels(sfiles);
+    }
+
+    if (sfiles->stage == IPP_STAGE_RAW) {
+        // copy PHU to output files
+        copyPHU(sfiles);
+
+        // advance to the first image extension
+        if (!streakFilesNextExtension(sfiles)) {
+            psErrorStackPrint(stderr, "failed to advance to first extension of input images");
+            exit (PS_EXIT_PROG_ERROR);
+        }
+    }
+
+    // Iterate through components of image
+    do {
+        bool exciseAll = false;
+
+        // read the images and copy the pixels image from the inputs to the outputs
+        // This also sets up the astrometry and initializes the recovery images
+        if (!readAndCopyToOutput(sfiles)) {
+            // this extension is not an image skip (video descriptor?) 
+            // XXX is there anything else that could be in an input file that we need to handle?
+            continue;
+        }
+
+        // Identify pixels to mask because of streaks
+
+        pixels = streak_on_component (streaks, sfiles->astrom,
+                                      sfiles->inImage->numCols, sfiles->inImage->numRows);
+        // XXX: 
+        // 
+        // PFS: Need to get mask weight for PM_MASK_DETECTOR.  Is this stored
+        // in a config somewhere?  Not sure how to properly set the maskValue.
+        // WES just pick some numbers
+#define MASK_STREAK     42
+#define MASK_NOT_WARPED 51
+
+        exciseNonWarpedPixels(sfiles, MASK_NOT_WARPED);
+
+        if (sfiles->inImage->image) {
+            for (i = 0; i < psArrayLength (pixels); ++i) {
+                pixelPos = psArrayGet (pixels, i);
+
+                if (warpedPixel(sfiles, pixelPos)) {
+
+                    excisePixel(sfiles, pixelPos->x, pixelPos->y, MASK_STREAK);
+
+                } else {
+                    // This pixel was not included in any warp and has thus already excised above
+                }
             }
-        }
-#endif
-
-        // XXX: This loop following will SEGV for raw video Cells becuase the image is null
-        for (i = 0; i < psArrayLength (pixels); ++i) {
-            pixelPos = psArrayGet (pixels, i);
-            imageValue  = psImageGet (sfiles->inImage->image,  pixelPos->x, pixelPos->y);
-            psImageSet (sfiles->recImage->image,  pixelPos->x, pixelPos->y, imageValue);
-            psImageSet (sfiles->outImage->image,  pixelPos->x, pixelPos->y, NAN);
-            if (sfiles->inWeight) {
-                weightValue = psImageGet (sfiles->inWeight->image, pixelPos->x, pixelPos->y);
-                psImageSet (sfiles->recWeight->image, pixelPos->x, pixelPos->y, weightValue);
-                psImageSet (sfiles->outWeight->image, pixelPos->x, pixelPos->y, NAN);
-            }
-            if (sfiles->inMask) {
-                maskValue   = psImageGet (sfiles->inMask->image,   pixelPos->x, pixelPos->y);
-                psImageSet (sfiles->recMask->image,   pixelPos->x, pixelPos->y, maskValue);
-                // TODO:
-                // Need to get mask weight for PM_MASK_DETECTOR.  Is this stored
-                // in a config somewhere?  Not sure how to properly set the maskValue.
-                psImageSet (sfiles->outMask->image,   pixelPos->x, pixelPos->y, maskValue);
-            }
+        } else { 
+            // This video cell is touched by a streak excise
+            // For now we excise all video cells.
+            // if (psArrayLength(pixels))
+            exciseAll = true;
         }
         psArrayElementsFree (pixels);
 
         // write out the destreaked temporary images and the recovery images
-        writeImages(sfiles);
+        writeImages(sfiles, exciseAll);
+
     } while (streakFilesNextExtension(sfiles));