Changeset 2783

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.

Location:

trunk/stac/src

Files:

• stac.c (modified) (8 diffs)
• stac.h (modified) (3 diffs)
• stacAreaOfInterest.c (modified) (1 diff)
• stacCombine.c (modified) (6 diffs)
• stacConfig.c (modified) (5 diffs)
• stacTransform.c (modified) (7 diffs)

trunk/stac/src/stac.c

@@ -36 +36 @@
     psTraceSetLevel("stac.rejection",10);
     psTraceSetLevel("stac.time",10);
-    
+    psTraceSetLevel("stac.area",10);
 
     // Set logging level
@@ -50 +50 @@
     psArray *maps = stacReadMaps(config);
 
+    psTrace("stac.time", 1, "I/O completed at %f seconds\n", getTime() - startTime);
+
     // Invert maps
     psArray *inverseMaps = stacInvertMaps(maps);
@@ -57 +59 @@
 
     // Transform inputs and errors
+    psArray *transformed = NULL;
     psArray *transformedErrors = NULL;
-    psArray *transformed = stacTransform(inputs, inverseMaps, errors, &transformedErrors, NULL, config);
+    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, NULL, NULL, config);
 
     psTrace("stac.time",1,"Transformation completed at %f seconds\n", getTime() - startTime);
@@ -64 +67 @@
     // Combine with rejection
     psArray *rejected = NULL;
-    psImage *combined = stacCombine(transformed, transformedErrors, 2,  &rejected, config);
+    psImage *combined = NULL;
+    (void)stacCombine(&combined, &rejected, transformed, transformedErrors, config->nReject, NULL, config);
+
+    psTrace("stac.time",1,"First combination completed at %f seconds\n", getTime() - startTime);
 
 #ifdef TESTING
@@ -75 +81 @@
     psArray *regions = psArrayAlloc(inputs->n); // Array of images denoting regions of interest
     for (int i = 0; i < inputs->n; i++) {
-        fprintf(stderr, "Finding area of interest for image %d\n", i);
         regions->data[i] = stacAreaOfInterest(rejected->data[i], inverseMaps->data[i],
                                               ((psImage*)(inputs->data[i]))->numCols,
@@ -89 +94 @@
     psArray *rejectedSource = stacRejection(inputs, rejected, regions, maps, inverseMaps, config);
 
-    // Redo transformation with the mask
-    psFree(transformed);
-    psFree(transformedErrors);
-    psArray *transformedErrorsNew = NULL;
-    psArray *transformedNew = stacTransform(inputs, inverseMaps, errors, &transformedErrorsNew,
-                                            rejectedSource, config);
+    // Get regions of interest in the output frame
+    psImage *combineRegion = psImageAlloc(config->outnx, config->outny, PS_TYPE_U8);
+    for (int y = 0; y < config->outny; y++) {
+        for (int x = 0; x < config->outnx; x++) {
+            combineRegion->data.U8[y][x] = 0;
+        }
+    }
+    for (int i = 0; i < inputs->n; i++) {
+        psImage *region = stacAreaOfInterest(rejectedSource->data[i], maps->data[i], config->outnx,
+                                             config->outny);
+        psBinaryOp(combineRegion, combineRegion, "+", region);
+        psFree(region);
+    }
+
+    // Redo transformation with the masks
+    (void)stacTransform(&transformed, &transformedErrors, inputs, inverseMaps, errors, rejectedSource,
+                        combineRegion, config);
 
     // Combine the newly-transformed CR-free images, no rejection
-    psFree(combined);
     psFree(rejected);
     rejected = NULL;
-    combined = stacCombine(transformedNew, transformedErrorsNew, 0,  &rejected, config);
+    (void)stacCombine(&combined, &rejected, transformed, transformedErrors, 0, combineRegion, config);
 
     // Write output image
@@ -107 +122 @@
     // Free everything I've used
     stacConfigFree(config);
+    psFree(combineRegion);
     psFree(regions);
     psFree(inputs);
@@ -113 +129 @@
     psFree(errors);
     psFree(transformedErrors);
-    psFree(transformedNew);
+    psFree(transformed);
     psFree(rejectedSource);
     psFree(combined);
 
-    psTrace("stac.time",1,"Combination completed at %f seconds\n", getTime() - startTime);
+    psTrace("stac.time",1,"Final combination completed at %f seconds\n", getTime() - startTime);
 
    // Check memory for leaks, corruption

trunk/stac/src/stac.h

@@ -28 +28 @@
                                         // considered bad
     float grad;                         // Multiplier of the gradient
+    int nReject;                        // Number of rejection iterations
 } stacConfig;
 
@@ -67 +68 @@
 // stacTransform.c
 
-// Transform input images
-psArray *stacTransform(const psArray *images, // Array of images to be transformed
-                       const psArray *maps, // Array of polynomials that do the transformation
-                       const psArray *errors, // Array of error images
-                       psArray **outErrors, // Array of output errors
-                       const psArray *masks, // Masks of input images
-                       const stacConfig *config // Configuration
+// Transform input images, return success
+bool stacTransform(psArray **outputs,   // Transformed images for output
+                   psArray **outErrors, // Transformed error images for output
+                   const psArray *images, // Array of images to be transformed
+                   const psArray *maps, // Array of polynomials that do the transformation
+                   const psArray *errors, // Array of error images to be transformed
+                   const psArray *masks, // Masks of input images
+                   const psImage *region, // Region of interest for transformation
+                   const stacConfig *config // Configuration
     );
 
@@ -110 +113 @@
 
 // Combine the transformed images
-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
     );
 

trunk/stac/src/stacAreaOfInterest.c

@@ -89 +89 @@
     }
 
+    psTrace("stac.area", 3, "Finding area of interest....\n");
+
+
     // Coordinates for the transformations
     psPlane *inFrame = psAlloc(sizeof(psPlane));

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

trunk/stac/src/stacConfig.c

@@ -9 +9 @@
 {
     fprintf (stderr, "STAC: Simultaneous Telescope Array Combination\n"
-             "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-s SAT] [-b BAD] [-k REJ] [-k FRAC] [-G GRAD] OUT IN1 IN2...\n"
+             "Usage: %s [-h] [-v] [-g GAIN] [-r RN] [-o NX NY] [-s SAT] [-b BAD] [-k REJ] [-n NREJECT] [-k FRAC] [-G GRAD] OUT IN1 IN2...\n"
              "where\n"
              "\t-h           Help (this info)\n"
@@ -19 +19 @@
              "\t-b BAD       Bad level (0)\n"
              "\t-k REJ       Rejection level (k-sigma; 3.5)\n"
+             "\t-n NREJECT   Number of rejection iterations (2)\n"
              "\t-f FRAC      Fraction of pixel to be marked before considered bad (0.5)\n"
              "\t-d GRAD      Gradient threshold for pixel to be marked (0.0)\n"
@@ -45 +46 @@
     config->frac = 0.5;
     config->grad = 0.4;
+    config->nReject = 2;
 
     return config;
@@ -74 +76 @@
 
     /* Parse command-line arguments using getopt */
-    while ((opt = getopt(argc, argv, "hvg:r:o:s:b:k:f:G:")) != -1) {
+    while ((opt = getopt(argc, argv, "hvg:r:o:s:b:k:n:f:G:")) != -1) {
         switch (opt) {
           case 'h':
@@ -115 +117 @@
           case 'k':
             if (sscanf(optarg, "%f", &config->reject) != 1) {
+                help(programName);
+            }
+            break;
+          case 'n':
+            if (sscanf(optarg, "%d", &config->nReject) != 1) {
                 help(programName);
             }

trunk/stac/src/stacTransform.c

@@ -105 +105 @@
 
 
-psArray *stacTransform(const psArray *images, // Array of images to be transformed
-                       const psArray *maps, // Array of polynomials that do the transformation
-                       const psArray *errors, // Array of error images to be transformed
-                       psArray **outErrors, // Transformed error images for output
-                       const psArray *masks, // Masks of input images
-                       const stacConfig *config // Configuration
+bool stacTransform(psArray **outputs,   // Transformed images for output
+                   psArray **outErrors, // Transformed error images for output
+                   const psArray *images, // Array of images to be transformed
+                   const psArray *maps, // Array of polynomials that do the transformation
+                   const psArray *errors, // Array of error images to be transformed
+                   const psArray *masks, // Masks of input images
+                   const psImage *region, // Region of interest for transformation
+                   const stacConfig *config // Configuration
     )
 {
@@ -121 +123 @@
         psError("stac.transform", "Number of maps (%d) does not match number of images (%d).\n",
                 maps->n, images->n);
-        return NULL;
+        return false;
     }
     if (errors && (images->n != errors->n)) {
         psError("stac.transform", "Number of error images (%d) does not match number of images (%d).\n",
                 errors->n, images->n);
-        return NULL;
-    }
-
+        return false;
+    }
+
+    // Allocate the output images if required, otherwise check the number
+    if (*outputs == NULL) {
+        *outputs = psArrayAlloc(nImages);
+        psTrace("stac.transform", 5, "Allocating space for transformed images, %dx%d\n", nx, ny);
+        for (int i = 0; i < nImages; i++) {
+            (*outputs)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32);
+        }
+    } else if ((*outputs)->n != nImages) {
+        psError("stac.transform", "Number of output images (%d) does not match number of input images "
+                "(%d).\n", (*outputs)->n, nImages);
+        return false;
+    }
+
+    // Allocate the output error images, if required, otherwise check the number
     if (errors && (*outErrors == NULL)) {
-        // Allocate the output error images, if required
         *outErrors = psArrayAlloc(errors->n);
+        psTrace("stac.transform", 5, "Allocating space for transformed error images, %dx%d\n", nx, ny);
+        for (int i = 0; i < nImages; i++) {
+            (*outErrors)->data[i] = psImageAlloc(nx, ny, PS_TYPE_F32);
+        }
     } else if (errors->n != (*outErrors)->n) {
         psError("stac.transform", "Number of error output images (%d) does not match number of input"
                 "images (%d).\n", (*outErrors)->n, errors->n);
-        return NULL;
-    }
-
+        return false;
+    }
+
+    // Check the masks, if specified
     if (masks != NULL) {
         if (masks->n != nImages) {
             psError("stac.transform", "Number of masks (%d) does not match number of input images (%d).\n",
                     masks->n, nImages);
-            return NULL;
+            return false;
         }
         for (int i = 0; i < nImages; i++) {
@@ -151 +171 @@
                          "Size of input mask (%dx%d) does not match that of input image (%dx%d).\n",
                          mask->numCols, mask->numRows, image->numCols, image->numRows);
-                return NULL;
+                return false;
             }
         }
     }
 
-    // Arrays of images for return
-    psArray *outImages = psArrayAlloc(nImages); // Output images
 
     // Stuff for the transformations
     psPlane *detector = psAlloc(sizeof(psPlane)); // Coordinates on the detector
     psPlane *sky = psAlloc(sizeof(psPlane)); // Coordinates on the sky
+
 
     // Iterate over the images
@@ -171 +190 @@
         psPlaneTransform *map = maps->data[n]; // The map
         psImage *error = errors->data[n]; // The error image
-
-        // Initialise the output images
-        psImage *outImage = psImageAlloc(nx, ny, PS_TYPE_F32);
-        psImage *outError = psImageAlloc(nx, ny, PS_TYPE_F32);
-        psTrace("stac.transform", 5, "Allocating space for transformed image, %dx%d\n", nx, ny);
+        psImage *outImage = (*outputs)->data[n]; // The output image
+        psImage *outError = (*outErrors)->data[n]; // The output error image
 
 #if 0
@@ -193 +209 @@
         }
 
-#if 0
-        // Calculate scales; could be adapted for higher order than linear by moving this into the
-        // pixel-dependent section and calculating derivatives
-        double xscale = 0.5*sqrt(SQUARE(map->x->coeff[0][1]) + SQUARE(map->x->coeff[1][0]));
-        double yscale = 0.5*sqrt(SQUARE(map->y->coeff[0][1]) + SQUARE(map->y->coeff[1][0]));
-        double areascale = 0.25 / xscale / yscale;
-#endif
-
         // Iterate over the output image pixels
         for (int y = 0; y < ny; y++) {
             for (int x = 0; x < nx; x++) {
-                // Transform!
-                sky->x = (double)x + 0.5;
-                sky->y = (double)y + 0.5;
-                (void)psPlaneTransformApply(detector, map, sky);
-
-                // Change PS_INTERPOLATE_BILINEAR to best available technique.
-                outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x + 0.5,
-                                                                          detector->y + 0.5, mask, 1, 0.0,
-                                                                          PS_INTERPOLATE_BILINEAR);
-                outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
-                                                                                        detector->x + 0.5,
-                                                                                        detector->y + 0.5,
-                                                                                        mask, 1, 0.0);
-                outError->data.F32[y][x] = sqrtf(outError->data.F32[y][x]);
+
+                // Only transform those pixels requested
+                if (!region || (region && region->data.U8[y][x])) {
+                    // Transform!
+                    sky->x = (double)x + 0.5;
+                    sky->y = (double)y + 0.5;
+                    (void)psPlaneTransformApply(detector, map, sky);
+                    
+                    // Change PS_INTERPOLATE_BILINEAR to best available technique.
+                    outImage->data.F32[y][x] = (psF32)psImagePixelInterpolate(image, detector->x + 0.5,
+                                                                              detector->y + 0.5, mask, 1, 0.0,
+                                                                              PS_INTERPOLATE_BILINEAR);
+                    // Error is actually the variance
+                    outError->data.F32[y][x] = (psF32)p_psImageErrorInterpolateBILINEAR_F32(error,
+                                                                                            detector->x + 0.5,
+                                                                                            detector->y + 0.5,
+                                                                                            mask, 1, 0.0);
+                } // Pixels of interest
+
             }
         } // Iterating over output pixels
@@ -233 +246 @@
 #endif
 
-        // Plug the new images into arrays
-        outImages->data[n] = outImage;
-        (*outErrors)->data[n] = outError;
-
     } // Iterating over images
 
@@ -243 +252 @@
     psFree(sky);
 
-    return outImages;
+    return true;
 }