Changeset 2661 for trunk/stac/src/stacRejection.c

Timestamp:

Dec 7, 2004, 3:38:05 PM (22 years ago)

Author:

Paul Price

Message:

Working, but not getting faint CRs

File:

• trunk/stac/src/stacRejection.c (modified) (3 diffs)

trunk/stac/src/stacRejection.c

@@ -2 +2 @@
 #include "pslib.h"
 #include "stac.h"
+
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+
+float stacGradient(psImage *image,      // Input for which to measure the gradient
+                   int x, int y         // Coordinates at which to measure the gradient
+    )
+{
+    float sum = 0.0;                    // The sum of surrounding pixels
+    int num = 0;
+    psVector *pixels = psVectorAlloc(9, PS_TYPE_F32); // Array of pixels
+    psVector *mask = psVectorAlloc(9, PS_TYPE_U8); // Corresponding mask
+
+    // Get limits
+    int xMin = MAX(x - 1, 0);
+    int xMax = MIN(x + 1, image->numCols - 1);
+    int yMin = MAX(y - 1, 0);
+    int yMax = MIN(y + 1, image->numRows - 1);
+    for (int j = yMin; j <= yMax; j++) {
+        for (int i = xMin; i <= xMax; i++) {
+            sum += image->data.F32[j][i];
+            pixels->data.F32[num] = image->data.F32[j][i];
+            mask->data.U8[num] = 1;
+            num++;
+        }
+    }
+#if 0
+    sum -= image->data.F32[y][x];
+    sum /= (float)(num-1);
+//    sum -= image->data.F32[y][x];
+    return sum / image->data.F32[y][x];
+#endif 
+
+    // Fill out the array
+    for (int i = num; i < 9; i++) {
+        pixels->data.F32[i] = 0;
+        mask->data.U8[i] = 0;
+    }
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+    (void)psVectorStats(stats, pixels, mask, 0);
+    float median = stats->sampleMedian;
+    psFree(stats);
+    psFree(pixels);
+    psFree(mask);
+    return median / image->data.F32[y][x];
+}
+
+#if 0
+float stacGradient(psImage *image,      // Input for which to measure the gradient
+                   int x, int y         // Coordinates at which to measure the gradient
+                   )
+{
+    float sum = 0.0;                    // The sum of surrounding pixels
+    float maxDiff = 0.0;
+    int num = 0;
+    // Get limits
+    int xMin = MAX(x - 1, 0);
+    int xMax = MIN(x + 1, image->numCols - 1);
+    int yMin = MAX(y - 1, 0);
+    int yMax = MIN(y + 1, image->numRows - 1);
+    for (int j = yMin; j <= yMax; j++) {
+        for (int i = xMin; i <= xMax; i++) {
+            if (image->data.F32[j][i] - image->data.F32[y][x] < maxDiff) {
+                maxDiff = image->data.F32[j][i] - image->data.F32[y][x];
+            }
+            sum += image->data.F32[j][i];
+            num++;
+        }
+    }
+    sum -= image->data.F32[y][x];
+    sum /= (float)(num-1);
+    sum -= image->data.F32[y][x];
+//    return sum / image->data.F32[y][x];
+    return maxDiff / image->data.F32[y][x];
+}   
+#endif
 
 
@@ -55 +132 @@
         int nxInput = ((psImage*)(inputs->data[i]))->numCols;
         int nyInput = ((psImage*)(inputs->data[i]))->numRows;
-        psImage *mask = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The mask in the source frame
+        psImage *mask = psImageAlloc(nxInput, nyInput, PS_TYPE_U8); // The pixel mask for the input
+#ifdef TESTING
+        psImage *rejmap = psImageAlloc(nxInput, nyInput, PS_TYPE_F32); // The rejections in the source frame
+        psImage *grad = psImageAlloc(nxInput, nyInput, PS_TYPE_F32);    // The gradient image
+#endif
         psImage *reject = rejected->data[i]; // Pull out the mask in the output frame
         psPlaneTransform *map = maps->data[i]; // The map from input to output
 
         psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i);
+        int nBad = 0;                   // Number of bad pixels
+
+#ifdef CRFLUX
+        // Set up CR output
+        FILE *crs = NULL;               // File for outputting details of rejected pixels
+        char crfile[MAXCHAR];   // Filename
+        sprintf(crfile,"%s.cr",config->inputs->data[i]);
+        if ((crs = fopen(crfile, "w")) == NULL) {
+            psError("stac.rejection","Unable to open file for detailed output, %s\n");
+            return NULL;
+        }
+#endif
 
         // Transform the mask
@@ -67 +160 @@
         for (int y = 0; y < nyInput; y++) {
             for (int x = 0; x < nxInput; x++) {
-                inCoords->x = (double)x + 0.5;
-                inCoords->y = (double)y + 0.5;
+                inCoords->x = (double)x;
+                inCoords->y = (double)y;
                 (void)psPlaneTransformApply(outCoords, map, inCoords);
-                mask->data.F32[y][x] = psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
+                float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
                                                                NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
+#ifdef TESTING
+                rejmap->data.F32[y][x] = maskVal;
+                grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y);
+#endif
+                // Check threshold, then check gradient
+                if ((maskVal > config->frac) && (stacGradient(inputs->data[i], x, y) < config->grad)) {
+                    mask->data.U8[y][x] = 1;
+                    nBad++;
+#ifdef CRFLUX
+                    fprintf(crs, "%d %d --> %f %f %f\n", x, y, ((psImage*)(inputs->data[i]))->data.F32[y][x],
+                            maskVal, stacGradient(inputs->data[i], x, y));
+#endif
+                } else {
+                    mask->data.U8[y][x] = 0;
+                }
 
             }
-        }
+        } // Iterating over pixels
+
+#ifdef CRFLUX
+        // Close file used for detailed output
+        fclose(crs);
+#endif
+
+        psTrace("stac.rejection", 1, "%d pixels masked in image %d.\n", nBad, i);
+        // Clip the image, and convert to suitable mask format
 
 #ifdef TESTING
         // Write error image out to check
-        char maskfile[MAXCHAR]; // Filename of mask image
+        char maskfile[MAXCHAR];         // Filename of mask image
+        char rejmapfile[MAXCHAR];       // Filename of rejection image
+        char gradfile[MAXCHAR];         // Filename of gradient image
         sprintf(maskfile,"%s.mask",config->inputs->data[i]);
+        sprintf(rejmapfile,"%s.rejmap",config->inputs->data[i]);
+        sprintf(gradfile,"%s.grad",config->inputs->data[i]);
         psImageWriteSection(mask,0,0,0,NULL,0,maskfile);
-#endif
-
-        // Clip the image, and convert to suitable mask format
-        (void)psImageClip(mask, config->frac, 0.0, config->frac, 1.0);
-        psImage *maskU8 = psImageCopy(NULL, mask, PS_TYPE_U8);
-        psFree(mask);
-    
+        psImageWriteSection(rejmap,0,0,0,NULL,0,rejmapfile);
+        psImageWriteSection(grad,0,0,0,NULL,0,gradfile);
+        psFree(rejmap);
+        psFree(grad);
+#endif
+
         // Stuff into the array
-        inputRej->data[i] = maskU8;
-
-    }
+        inputRej->data[i] = mask;
+
+    }
+
 
     psFree(inCoords);