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

Timestamp:

Dec 8, 2004, 3:12:31 PM (22 years ago)

Author:

Paul Price

Message:

More robust gradient measurement

File:

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

trunk/stac/src/stacRejection.c

@@ -11 +11 @@
     )
 {
-    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
+    psVector *pixels = psVectorAlloc(8, PS_TYPE_F32); // Array of pixels
+    psVector *mask = psVectorAlloc(8, PS_TYPE_U8); // Corresponding mask
 
     // Get limits
@@ -23 +22 @@
     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;
-    }
+            if ((i != x) && (j != y)) {
+                pixels->data.F32[num] = image->data.F32[j][i];
+                mask->data.U8[num] = 0;
+                num++;
+            }
+        }
+    }
+    pixels->n = num;
+    mask->n = num;
+
+    // Get the median
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    (void)psVectorStats(stats, pixels, mask, 0);
+    (void)psVectorStats(stats, pixels, mask, 1);
     float median = stats->sampleMedian;
     psFree(stats);
@@ -85 +77 @@
                        psImage *combined, // Combined image
                        psArray *maps,   // Maps from input to transformed image
+                       psArray *inverseMaps, // Maps from transformed to input image
                        psArray *rejected, // Rejected images
                        stacConfig *config // Configuration
@@ -167 +160 @@
 #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++;
+#endif
+
+                // Check gradient
+                float meanGrads = 0.0;
+                int numGrads = 0;
+                for (int j = 0; j < nImages; j++) {
+                    if (i != j) {
+                        // Get coordinates for that image
+                        (void)psPlaneTransformApply(inCoords, inverseMaps->data[j], outCoords);
+                        int xPix = (int)(inCoords->x + 0.5);
+                        int yPix = (int)(inCoords->y + 0.5);
+                        if ((xPix >= 0) && (xPix <= ((psImage*)(inputs->data[j]))->numCols - 1) &&
+                            (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
+                            // Calculate the gradient
+                            meanGrads += stacGradient(inputs->data[j], xPix, yPix);
+                            numGrads++;
+                        }
+                    }
+                }
+                if (numGrads > 0) {
+                    meanGrads /= (float)numGrads;
+                } else {
+                    meanGrads = 0;
+                }
+                
+#ifdef TESTING
+                grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
+#endif
+                if (maskVal > config->frac) {
+
+                    if ((maskVal > config->frac) &&
+                        (stacGradient(inputs->data[i], x, y)) < config->grad * meanGrads) {
+                        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));
+                        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 {
@@ -189 +211 @@
 #endif
 
-        psTrace("stac.rejection", 1, "%d pixels masked in image %d.\n", nBad, i);
+        psTrace("stac.rejection", 1, "%d pixels masked in image %d\n", nBad, i);
         // Clip the image, and convert to suitable mask format