Changeset 2783 for trunk/stac/src/stacTransform.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/stacTransform.c (modified) (7 diffs)

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;
 }