Changeset 2783 for trunk/stac/src/stacCombine.c

Timestamp:

Dec 21, 2004, 4:07:04 PM (22 years ago)

Author:

Paul Price

Message:

Running fast now. Only optimisation tweaks to be made now.

File:

• trunk/stac/src/stacCombine.c (modified) (6 diffs)

trunk/stac/src/stacCombine.c

@@ -26 +26 @@
     for (int i = 0; i < num; i++) {
         if (! masks->data.U8[i]) {
-            sum += values->data.F32[i] / SQUARE(errors->data.F32[i]);
-            weights += 1.0 / SQUARE(errors->data.F32[i]);
+            // "error" here is the variance
+            sum += values->data.F32[i] / errors->data.F32[i];
+            weights += 1.0 / errors->data.F32[i];
         }
     }
@@ -53 +54 @@
 
 
-psImage *stacCombine(psArray *images,   // Array of transformed images
-                     psArray *errors,   // Array of transformed error images
-                     int nReject,       // Number of rejection iterations
-                     psArray **rejected, // Array of rejection masks
-                     stacConfig *config // Configuration
+
+bool stacCombine(psImage **combined,    // The combined image for output
+                 psArray **rejected,    // Array of rejection masks
+                 psArray *images,       // Array of transformed images
+                 psArray *errors,       // Array of transformed error images
+                 int nReject,           // Number of rejection iterations
+                 psImage *region,       // Region to combine
+                 stacConfig *config     // Configuration
     )
 {
@@ -76 +80 @@
                     "Image size mismatch on error image %d\nExpected %dx%d, got %dx%d\n",
                     i, numCols, numRows, image->numCols, image->numRows);
-            return NULL;
+            return false;
         }
         if ((error->numCols != numCols) || (error->numRows != numRows)) {
@@ -82 +86 @@
                     "Image size mismatch on error image %d\nExpected %dx%d, got %dx%d\n",
                     i, numCols, numRows, error->numCols, error->numRows);
-            return NULL;
-        }
+            return false;
+        }
+    }
+
+    // Check combined image
+    if (*combined == NULL) {
+        *combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
+    } else if (((*combined)->numRows != numRows) || ((*combined)->numCols != numCols)) {
+        psError("stac.combine", "Size of combined image (%dx%d) does not match inputs (%dx%d)\n",
+                (*combined)->numCols, (*combined)->numRows, numCols, numRows);
+        return false;
+    }
+
+    // Check area of interest
+    if (region && ((region->numRows != numRows) || (region->numCols != numCols))) {
+        psError("stac.combine", "Size of area of interest (%dx%d) does not match inputs (%dx%d)\n",
+                region->numCols, region->numRows, numCols, numRows);
+        return false;
     }
 
     psTrace("stac.combine", 1, "Combining images....\n");
-    psTrace("stac.combine", 3, "Saturation: %f Bad: %f\n", saturated, bad);
 
     psVector *pixels = psVectorAlloc(nImages, PS_TYPE_F32); // Will hold the pixels in the statistics step
     psVector *deltas = psVectorAlloc(nImages, PS_TYPE_F32); // Will hold the errors in the statistics step
-    psVector *mask = psVectorAlloc(nImages, PS_TYPE_U8); // Mask bad pixels
-    psImage *combined = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Combined image
-
+    psVector *mask = psVectorAlloc(nImages, PS_TYPE_U8); // Will hold the mask in the statistics step
 
     // Set up rejection masks
@@ -105 +122 @@
             return NULL;
         }
+
+        // Create and initialise rejection masks
         for (int i = 0; i < nImages; i++) {
-            (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8); // Create rejection mask
-            // Zero out the mask
+            (*rejected)->data[i] = psImageAlloc(numCols, numRows, PS_TYPE_U8);
             for (int r = 0; r < numRows; r++) {
                 for (int c = 0; c < numCols; c++) {
@@ -124 +142 @@
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
+
+            // Only combine those pixels requested
+            if (!region || (region && region->data.U8[y][x])) {
             
-            // Export pixels into the vector and get stats
-            for (int i = 0; i < nImages; i++) {
-                float pixel = ((psImage*)images->data[i])->data.F32[y][x];
-                float delta = ((psImage*)errors->data[i])->data.F32[y][x];
-                pixels->data.F32[i] = pixel;
-                deltas->data.F32[i] = delta;
-                if ((pixel >= saturated) || (pixel <= bad) || (! isfinite(pixel)) || (! isfinite(delta))) {
-                    mask->data.U8[i] = (psU8)1; // Don't use!
-                } else {
-                    mask->data.U8[i] = (psU8)0; // Use.
-                }
-            }
+                // Export pixels into the vector and get stats
+                for (int i = 0; i < nImages; i++) {
+                    float pixel = ((psImage*)images->data[i])->data.F32[y][x];
+                    float delta = ((psImage*)errors->data[i])->data.F32[y][x]; // This is the variance!
+                    pixels->data.F32[i] = pixel;
+                    deltas->data.F32[i] = delta;
+                    if ((pixel >= saturated) || (pixel <= bad) || (! isfinite(pixel)) || (! isfinite(delta))) {
+                        mask->data.U8[i] = (psU8)1; // Don't use!
+                    } else {
+                        mask->data.U8[i] = (psU8)0; // Use.
+                    }
+                }
             
-            float average = stacCombineMean(pixels, deltas, mask); // Combined value
+                float average = stacCombineMean(pixels, deltas, mask); // Combined value
 
 #ifdef TESTING
-            // Calculate chi^2
-            chi2->data.F32[y][x] = 0.0;
-            int numGoodPix = 0;
-            for (int i = 0; i < nImages; i++) {
-                if (! mask->data.U8[i]) {
-                    chi2->data.F32[y][x] += SQUARE((pixels->data.F32[i] - average) / deltas->data.F32[i]);
-                    numGoodPix++;
-                }
-            }
-            chi2->data.F32[y][x] /= (float)numGoodPix;
-#endif
+                // Calculate chi^2
+                chi2->data.F32[y][x] = 0.0;
+                int numGoodPix = 0;
+                for (int i = 0; i < nImages; i++) {
+                    if (! mask->data.U8[i]) {
+                        chi2->data.F32[y][x] += SQUARE(pixels->data.F32[i] - average) / deltas->data.F32[i];
+                        numGoodPix++;
+                    }
+                }
+                chi2->data.F32[y][x] /= (float)numGoodPix;
+#endif
+                
+                // Rejection iterations
+                bool keepGoing = true;  // Keep going with rejection?
+                for (int rejNum = 0; (rejNum < nReject) && keepGoing; rejNum++) {
+                    float max = 0.0;    // Maximum deviation
+                    int maxIndex = -1;  // Index of the maximum deviation
+                    for (int i = 0; i < nImages; i++) {
+                        if (!mask->data.U8[i] &&
+                            ((pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]) > max)) {
+                            max = (pixels->data.F32[i] - average) / sqrtf(deltas->data.F32[i]);
+                            maxIndex = i;
+                        }
+                    }
+                    // Reject the pixel with the maximum deviation
+                    if (max > reject) {
+                        mask->data.U8[maxIndex] = 1;
+                        ((psImage*)((*rejected)->data[maxIndex]))->data.U8[y][x] += 1;
+                        // Re-do combination following rejection
+                        average = stacCombineMean(pixels, deltas, mask);
+                    } else {
+                        keepGoing = false;
+                    }
+                }
             
-            // Rejection iterations
-            bool keepGoing = true;      // Keep going with rejection?
-            for (int rejNum = 0; (rejNum < nReject) && keepGoing; rejNum++) {
-                float max = 0.0;        // Maximum deviation
-                int maxIndex = -1;      // Index of the maximum deviation
-                for (int i = 0; i < nImages; i++) {
-                    if (!mask->data.U8[i] && ((pixels->data.F32[i] - average) / deltas->data.F32[i] > max)) {
-                        max = (pixels->data.F32[i] - average) / deltas->data.F32[i];
-                        maxIndex = i;
-                    }
-                }
-                // Reject the pixel with the maximum deviation
-                if (max > reject) {
-                    mask->data.U8[maxIndex] = 1;
-                    ((psImage*)((*rejected)->data[maxIndex]))->data.U8[y][x] += 1.0;
-                    // Re-do combination following rejection
-                    average = stacCombineMean(pixels, deltas, mask);
-                } else {
-                    keepGoing = false;
-                }
-            }
-            
-            combined->data.F32[y][x] = average;
+                (*combined)->data.F32[y][x] = average;
+
+            } // Pixels of interest
+
         }
     } // Iterating over output pixels