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

Timestamp:

Dec 20, 2004, 6:15:28 PM (22 years ago)

Author:

Paul Price

Message:

Running faster now by only transforming rejection masks where we're interested. To do: only transform and combine in the second iteration for areas where we're interested

File:

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

trunk/stac/src/stacRejection.c

@@ -42 +42 @@
 }
 
-#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
-
-
 psArray *stacRejection(psArray *inputs, // Input images
-                       psArray *transformed, // Transformed images
-                       psImage *combined, // Combined image
+                       psArray *rejected, // Rejected images
+                       psArray *regions, // Regions of interest
                        psArray *maps,   // Maps from input to transformed image
                        psArray *inverseMaps, // Maps from transformed to input image
-                       psArray *rejected, // Rejected images
                        stacConfig *config // Configuration
     )
@@ -85 +53 @@
 
     // Check inputs
-    if (inputs->n != transformed->n) {
-        psError("stac.rejection", "Number of input (%d) and transformed (%d) images does not match.\n",
-                inputs->n, transformed->n);
+    if (inputs->n != rejected->n) {
+        psError("stac.rejection", "Number of input images (%d) and rejection masks (%d) does not match.\n",
+                inputs->n, rejected->n);
+        return NULL;
+    }
+    if (inputs->n != regions->n) {
+        psError("stac.rejection", "Number of input images (%d) and region masks (%d) does not match.\n",
+                inputs->n, regions->n);
+        return NULL;
+    }
+    if (inputs->n != maps->n) {
+        psError("stac.rejection", "Number of input images (%d) and transformation maps (%d) does "
+                "not match.\n", inputs->n, maps->n);
+        return NULL;
+    }
+    if (inputs->n != inverseMaps->n) {
+        psError("stac.rejection", "Number of input images (%d) and inverse transformation maps (%d) does "
+                "not match.\n", inputs->n, inverseMaps->n);
         return NULL;
     }
     for (int i = 0; i < nImages; i++) {
-        if ((((psImage*)(transformed->data[i]))->numRows != combined->numRows) ||
-            (((psImage*)(transformed->data[i]))->numCols != combined->numCols)) {
+        if ((((psImage*)(inputs->data[i]))->numRows != ((psImage*)(regions->data[i]))->numRows) ||
+            (((psImage*)(inputs->data[i]))->numCols != ((psImage*)(regions->data[i]))->numCols)) {
             psError("stac.rejection",
-                    "Sizes of transformed (%dx%d) and combined (%dx%d) images do not match.\n",
-                    ((psImage*)(transformed->data[i]))->numCols, ((psImage*)(transformed->data[i]))->numRows,
-                    combined->numCols, combined->numRows);
+                    "Sizes of input image (%dx%d) and region mask (%dx%d) for input %d do not match.\n",
+                    ((psImage*)(inputs->data[i]))->numCols, ((psImage*)(inputs->data[i]))->numRows,
+                    ((psImage*)(regions->data[i]))->numCols, ((psImage*)(regions->data[i]))->numRows, i);
             return NULL;
         }
-        if ((((psImage*)(transformed->data[i]))->numRows != ((psImage*)(rejected->data[i]))->numRows) ||
-            (((psImage*)(transformed->data[i]))->numCols != ((psImage*)(rejected->data[i]))->numCols)) {
-            psError("stac.rejection",
-                    "Sizes of transformed (%dx%d) and rejected (%dx%d) images do not match.\n",
-                    ((psImage*)(transformed->data[i]))->numCols, ((psImage*)(transformed->data[i]))->numRows,
-                    ((psImage*)(rejected->data[i]))->numCols, ((psImage*)(rejected->data[i]))->numRows);
-            return NULL;
-        }
-    }
-
-    // Size of output images
-    int nxOutput = combined->numCols;
-    int nyOutput = combined->numRows;
+    }
 
     // Stuff for the transformations
@@ -132 +103 @@
         psImage *reject = rejected->data[i]; // Pull out the mask in the output frame
         psPlaneTransform *map = maps->data[i]; // The map from input to output
+        psImage *region = regions->data[i]; // The region of interest for this image
 
         psTrace("stac.rejection", 3, "Transforming rejection mask %d....\n", i);
@@ -153 +125 @@
         for (int y = 0; y < nyInput; y++) {
             for (int x = 0; x < nxInput; x++) {
-                inCoords->x = (double)x;
-                inCoords->y = (double)y;
-                (void)psPlaneTransformApply(outCoords, map, inCoords);
-                float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
-                                                               NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
+
+                // Only transform pixels of interest
+                if (region->data.U8[y][x]) {
+
+                    inCoords->x = (double)x;
+                    inCoords->y = (double)y;
+                    (void)psPlaneTransformApply(outCoords, map, inCoords);
+                    float maskVal = (float)psImagePixelInterpolate(reject, outCoords->x, outCoords->y,
+                                                                   NULL, 0, 0.0, PS_INTERPOLATE_BILINEAR);
 #ifdef TESTING
-                rejmap->data.F32[y][x] = maskVal;
-#endif
-
-#ifdef TESTING
-                // 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++;
+                    rejmap->data.F32[y][x] = maskVal;
+
+                    // Calculate 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;
-                }
-#else
-                float meanGrads = 10000.0;
-#endif
-
-                
-#ifdef TESTING
-                grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
-#endif
-                if (maskVal > config->frac) {
+                    if (numGrads > 0) {
+                        meanGrads /= (float)numGrads;
+                    } else {
+                        meanGrads = 0;
+                    }
+                    grad->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
+#endif
 
                     if ((maskVal > config->frac) &&
@@ -200 +167 @@
                         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
+                        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;
                     }
-                } else {
-                    mask->data.U8[y][x] = 0;
-                }
+
+                } // Only touching pixels of interest
 
             }