Changeset 25406 for branches/eam_branches/20090715/magic/remove/src/streaksremove.c

Timestamp:

Sep 15, 2009, 4:02:42 PM (17 years ago)

Author:

eugene

Message:

updates from trunk

File:

• branches/eam_branches/20090715/magic/remove/src/streaksremove.c (modified) (19 diffs)

branches/eam_branches/20090715/magic/remove/src/streaksremove.c

@@ -12 +12 @@
 static pmConfig *parseArguments(int argc, char **argv);
 static bool readAndCopyToOutput(streakFiles *sf, bool exciseAll);
-static void exciseNonWarpedPixels(streakFiles *sfiles, double newMaskValue);
+static void exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue);
 static bool warpedPixel(streakFiles *sfiles, int x, int y);
-static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, double newMaskValue);
+static void excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, psImageMaskType newMaskValue);
 static void writeImages(streakFiles *sf, bool exciseImageCube);
 static void updateAstrometry(streakFiles *sfiles);
-static void censorSources(streakFiles *sfiles, psU32 maskStreak);
+static void censorSources(streakFiles *sfiles, psImageMaskType maskStreak);
+static long censorPixels(streakFiles *sfiles, psImage * pixels, bool checkNonWarpedPixels, psU16 maskStreak);
 
 int
@@ -35 +36 @@
     // 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 
+    psU32 maskMask = 0;             // value looked up for MASK.STREAK
 
     psString streaksFileName = psMetadataLookupStr(NULL, config->arguments, "STREAKS");
@@ -42 +43 @@
     Streaks *streaks = readStreaksFile(streaksFileName);
     if (!streaks) {
-        psError(PS_ERR_UNKNOWN, "failed to read streaks file: %s", streaksFileName);
+        psError(PS_ERR_UNKNOWN, false, "failed to read streaks file: %s", streaksFileName);
         streaksExit("", PS_EXIT_PROG_ERROR);
     }
@@ -83 +84 @@
         psLogMsg("streaksremove", PS_LOG_INFO, "time to compute warped pixels: %f\n", cwp_t);
     }
-    
+
     if (sfiles->stage == IPP_STAGE_RAW) {
         // Except for raw stage, all of our (GPC1) files have one image extension.
@@ -97 +98 @@
     }
 
-    int totalPixels = 0;
-    int totalStreakPixels = 0;
+    long totalPixels = 0;
+    long totalStreakPixels = 0;
 
     // Iterate through each component of the input (except for raw images there is only one)
@@ -132 +133 @@
                                                         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
@@ -146 +147 @@
                 }
 
-
                 psTimerStart("REMOVE_STREAKS");
 
-                for (int y=0 ; y < sfiles->inImage->numRows; y++) {
-                    for (int x = 0; x < sfiles->inImage->numCols; x++) {
-                        if (psImageGet(pixels, x, y)) {
-                            ++totalStreakPixels;
-                            if (!checkNonWarpedPixels || warpedPixel(sfiles, x, y)) {
-
-                                excisePixel(sfiles, x, y, true, maskStreak);
-
-                            } else {
-                                // This pixel was not included in any warp and has thus already excised
-                                // by exciseNonWarpedPixels
-                            }
-                        }
-                    }
-                }
+                totalStreakPixels += censorPixels(sfiles, pixels, checkNonWarpedPixels, maskStreak);
 
                 psLogMsg("streaksremove", PS_LOG_INFO, "time to remove streak pixels: %f\n", psTimerClear("REMOVE_STREAKS"));
@@ -172 +158 @@
                 }
 
-            } else { 
+            } else {
                 // this component contains an image cube
                 // For now excise it completely
@@ -215 +201 @@
 
     if (!replicateOutputs(sfiles)) {
-        psError(PS_ERR_UNKNOWN, false, "failed to replicate output files");
+        psErrorStackPrint(stderr, "failed to replicate output files");
         deleteTemps(sfiles);
-        psErrorStackPrint(stderr, "");
         exit(PS_EXIT_UNKNOWN_ERROR);
     }
@@ -247 +232 @@
     pmConceptsDone();
     pmModelClassCleanup();
-    streaksNebulousCleanup(); 
+    streaksNebulousCleanup();
     pmConfigDone();
     psLogMsg("streaksremove", PS_LOG_INFO, "time to run streaksremove: %f\n", psTimerClear("STREAKSREMOVE"));
@@ -255 +240 @@
 
     return 0;
+}
+
+static long
+censorPixels(streakFiles *sfiles, psImage *pixels, bool checkNonWarpedPixels, psU16 maskStreak)
+{
+    long streakPixels = 0;
+
+    for (int y=0 ; y < sfiles->inImage->numRows; y++) {
+        for (int x = 0; x < sfiles->inImage->numCols; x++) {
+            if (psImageGet(pixels, x, y)) {
+                ++streakPixels;
+                if (!checkNonWarpedPixels || warpedPixel(sfiles, x, y)) {
+
+                    excisePixel(sfiles, x, y, true, maskStreak);
+
+                } else {
+                    // This pixel was not included in any warp and has thus already excised
+                    // by exciseNonWarpedPixels
+                }
+            }
+        }
+    }
+    return streakPixels;
 }
 
@@ -356 +364 @@
             true);
     }
-    
+
     if ((argnum = psArgumentGet(argc, argv, "-keepnonwarped"))) {
         psArgumentRemove(argnum, &argc, argv);
@@ -374 +382 @@
         psArgumentRemove(argnum, &argc, argv);
     }
-        
+
     // if skycells are not provided then we have to execise all pixels  unless -keepnonwarped
     pmConfigFileSetsMD(config->arguments, &argc, argv, "SKYCELLS", "-skycell", "-skycelllist");
@@ -628 +636 @@
             initValue = NAN;
         } else {
-            // otherwise write it to the output 
+            // otherwise write it to the output
             writeImageCube(sf->outImage, sf->inImage->imagecube, extname, sf->extnum);
             initValue = 0;
@@ -666 +674 @@
 
 static void
-excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, double newMaskValue)
+excisePixel(streakFiles *sfiles, unsigned int x, unsigned int y, bool streak, psImageMaskType newMaskValue)
 {
     double exciseValue = sfiles->inImage->exciseValue;
@@ -675 +683 @@
     }
 
-    double imageValue  = psImageGet (sfiles->inImage->image,  x, y);
+    float imageValue  = sfiles->inImage->image->data.F32[y][x];
     if (sfiles->recImage && !isExciseValue(imageValue, sfiles->inImage->exciseValue) ) {
-        psImageSet (sfiles->recImage->image,  x, y, imageValue);
+        sfiles->recImage->image->data.F32[y][x] = imageValue;
     }
 
     if (sfiles->transparentStreaks == 0) {
-        psImageSet (sfiles->outImage->image,  x, y, exciseValue);
+        sfiles->outImage->image->data.F32[y][x] = exciseValue;
     } else {
         if (streak) {
             // as a visualization aid don't mask the pixel, just change the intensity
-            psImageSet (sfiles->outImage->image,  x, y, imageValue + sfiles->transparentStreaks);
+            sfiles->outImage->image->data.F32[y][x] = imageValue + sfiles->transparentStreaks;
         } else {
-            psImageSet (sfiles->outImage->image,  x, y, exciseValue);
+            sfiles->outImage->image->data.F32[y][x] = exciseValue;
         }
     }
@@ -693 +701 @@
     if (sfiles->outWeight) {
         if (sfiles->recWeight) {
-            double weightValue = psImageGet (sfiles->inWeight->image, x, y);
-            psImageSet (sfiles->recWeight->image, x, y, weightValue);
+            sfiles->recWeight->image->data.F32[y][x] = sfiles->inWeight->image->data.F32[y][x];
         }
         // Assume that weight images are always a floating point type
-        psImageSet (sfiles->outWeight->image, x, y, NAN);
+        sfiles->outWeight->image->data.F32[y][x] = NAN;
     }
     if (sfiles->outMask) {
         if (sfiles->recMask) {
-            double maskValue   = psImageGet (sfiles->inMask->image,   x, y);
-            psImageSet (sfiles->recMask->image,   x, y, maskValue);
-        }
-        psImageSet (sfiles->outMask->image,   x, y, newMaskValue);
+            sfiles->recMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] =
+                sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x];
+        }
+        sfiles->outMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] =
+            sfiles->inMask->image->data.PS_TYPE_IMAGE_MASK_DATA[y][x] | newMaskValue;
     }
 }
 
 static void
-exciseNonWarpedPixels(streakFiles *sfiles, double newMaskValue)
+exciseNonWarpedPixels(streakFiles *sfiles, psImageMaskType newMaskValue)
 {
     int cell_x0 = sfiles->astrom->cell_x0;
@@ -785 +793 @@
 // streak mask
 static void
-censorSources(streakFiles *sfiles, psU32 maskStreak)
+censorSources(streakFiles *sfiles, psImageMaskType maskStreak)
 {
     if ((!sfiles->inSources) || (!sfiles->outMask)) {
@@ -799 +807 @@
     sFile *out = sfiles->outSources;
 
-    in->header = psFitsReadHeader(NULL, in->fits);
-    if (!in->header) {
-        psError(PS_ERR_IO, false, "failed to read header from %s", in->resolved_name);
+
+    // Primary header, should be "something.hdr"
+    {
+        psMetadata *header = psFitsReadHeader(NULL, in->fits);
+        if (!header) {
+            psError(PS_ERR_IO, false, "failed to read header from %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+
+        bool status;
+        psString extname = psMetadataLookupStr(&status, header, "EXTNAME");
+        if (!extname) {
+            psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        addDestreakKeyword(header);
+
+        if (!psFitsWriteBlank(out->fits, header, extname)) {
+            psError(PS_ERR_IO, false, "failed to write blank in header of %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        psFree(header);
+    }
+
+    // Extension with PSF fits, should be "something.psf"
+    {
+        if (!psFitsMoveExtNum(in->fits, 1, true)) {
+            psErrorStackPrint(stderr, "failed to read header from %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+
+        psMetadata *header = psFitsReadHeader(NULL, in->fits);
+        if (!header) {
+            psErrorStackPrint(stderr, "failed to read header from %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+        psString extname = psMetadataLookupStr(NULL, header, "EXTNAME");
+        if (!extname) {
+            psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+
+        psArray *inTable = psFitsReadTable(in->fits);
+        if (!inTable->n) {
+            psErrorStackPrint(stderr, "table in %s is empty", in->resolved_name);
+            streaksExit("", PS_EXIT_DATA_ERROR);
+        }
+
+        psArray *outTable = psArrayAllocEmpty(inTable->n);
+        int j = 0;
+        int numCensored = 0;
+        for (int i = 0 ; i < inTable->n; i++) {
+            psMetadata *row = inTable->data[i];
+
+            psF32 x = psMetadataLookupF32(NULL, row, "X_PSF");
+            psF32 y = psMetadataLookupF32(NULL, row, "Y_PSF");
+
+            psImageMaskType mask = maskImage->data.PS_TYPE_IMAGE_MASK_DATA[(int)y][(int)x];
+
+            // Key the source if the center pixel is not masked with maskStreak
+            if (!(mask & maskStreak) ) {
+                psArraySet(outTable, j++, row);
+            } else {
+                numCensored++;
+            }
+        }
+
+        // get rid of unused elements (don't know if this is necessary)
+        psArrayRealloc(outTable, j);
+
+        addDestreakKeyword(header);
+        if (psArrayLength(outTable) > 0) {
+            printf("Censored %d sources\n", numCensored);
+            if (! psFitsWriteTable(out->fits, header, outTable, extname)) {
+                psErrorStackPrint(stderr, "failed to write table to %s", out->resolved_name);
+                streaksExit("", PS_EXIT_DATA_ERROR);
+            }
+        } else {
+            printf("Censored ALL %d sources\n", numCensored);
+            if (! psFitsWriteTableEmpty(out->fits, header, inTable->data[0], extname)) {
+                psErrorStackPrint(stderr, "failed to write empty table to %s", out->resolved_name);
+                streaksExit("", PS_EXIT_DATA_ERROR);
+            }
+        }
+        psFree(header);
+        psFree(outTable);
+        psFree(inTable);
+    }
+
+    // XXX Will need to update to handle extension with extended sources, etc.
+
+    if (!psFitsClose(out->fits)) {
+        psErrorStackPrint(stderr, "failed to close table %s", out->resolved_name);
         streaksExit("", PS_EXIT_DATA_ERROR);
     }
-
-    bool status;
-    psString extname = psMetadataLookupStr(&status, in->header, "EXTNAME");
-    if (!extname) {
-        psError(PS_ERR_IO, false, "failed to find extname in header of %s", in->resolved_name);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-
-    psArray *inTable = psFitsReadTable(in->fits);
-    if (!inTable->n) {
-        psError(PS_ERR_IO, false, "table in %s is empty", in->resolved_name);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-
-    psArray *outTable = psArrayAllocEmpty(inTable->n);
-    int j = 0;
-    int numCensored = 0;
-    for (int i = 0 ; i < inTable->n; i++) {
-        psMetadata *row = inTable->data[i];
-
-        psF32 x = psMetadataLookupF32 (&status, row, "X_PSF");
-        psF32 y = psMetadataLookupF32 (&status, row, "Y_PSF");
-        
-        psU32 mask = psImageGet(maskImage, x, y);
-
-        // Key the source if the center pixel is not masked with maskStreak
-        if (! (mask & maskStreak) ) {
-            psArraySet(outTable, j++, row);
-        } else {
-            numCensored++;
-        }
-    }
-
-    // get rid of unused elements (don't know if this is necessary)
-    psArrayRealloc(outTable, j);
-
-    addDestreakKeyword(in->header);
-    if (psArrayLength(outTable) > 0) {
-        printf("Censored %d sources\n", numCensored);
-        if (! psFitsWriteTable(out->fits, in->header, outTable, extname)) {
-            psError(PS_ERR_IO, false, "failed to write table to %s", out->resolved_name);
-            streaksExit("", PS_EXIT_DATA_ERROR);
-        }
-    } else {
-        printf("Censored ALL %d sources\n", numCensored);
-        if (! psFitsWriteTableEmpty(out->fits, in->header, inTable->data[0], extname)) {
-            psError(PS_ERR_IO, false, "failed to write empty table to %s", out->resolved_name);
-            streaksExit("", PS_EXIT_DATA_ERROR);
-        }
-    }
-
-    if (!psFitsClose(out->fits)) {
-        psError(PS_ERR_IO, false, "failed to close table %s", out->resolved_name);
-        streaksExit("", PS_EXIT_DATA_ERROR);
-    }
-}
+}