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

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)
• streaksremove.c (modified) (18 diffs)

trunk/magic/remove/src

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

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

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;

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