Changeset 25027 for branches/pap/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Aug 7, 2009, 4:08:25 PM (17 years ago)

Author:

Paul Price

Message:

Merging trunk (r25026) to get up-to-date on old branch.

Location:

branches/pap

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmSubtraction.c (modified) (8 diffs)

branches/pap/psModules/src/imcombine/pmSubtraction.c

@@ -196 +196 @@
 {
     psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image
-    int numBytes = (2 * footprint + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
     for (int y = -footprint; y <= footprint; y++) {
-        // Convolution with a delta function is just the value specified by the offset
-        memcpy(&convolved->kernel[y][-footprint], &image->kernel[y - v][-footprint - u], numBytes);
+        for (int x = -footprint; x <= footprint; x++) {
+            convolved->kernel[y][x] = image->kernel[y - v][x - u];
+        }
     }
     return convolved;
 }
 
-// Take a subset of an image
-static inline psImage *convolveSubsetAlloc(psImage *image, // Image to be convolved
-                                           psRegion region, // Region of interest (with border)
-                                           bool threaded // Are we running threaded?
-                                           )
-{
-    if (!image) {
-        return NULL;
-    }
-    // XXX if (threaded) {
-    // XXX     psMutexLock(image);
-    // XXX }
-    psImage *subset = psImageSubset(image, region); // Subset image, to return
-    // XXX if (threaded) {
-    // XXX     psMutexUnlock(image);
-    // XXX }
-    return subset;
-}
-
-// Free the subset of an image
-static inline void convolveSubsetFree(psImage *parent, // Parent image
-                                      psImage *child, // Child (subset) image
-                                      bool threaded // Are we running threaded?
-                                      )
-{
-    if (!child || !parent) {
-        return;
-    }
-    // XXX if (threaded) {
-    // XXX     psMutexLock(parent);
-    // XXX }
-    psFree(child);
-    // XXX if (threaded) {
-    // XXX     psMutexUnlock(parent);
-    // XXX }
-    return;
-}
 
 // Convolve an image using FFT
@@ -256 +219 @@
                                   region.y0 - size, region.y1 + size); // Add a border
 
-    bool threaded = pmSubtractionThreaded(); // Are we running threaded?
-
-    psImage *subImage = convolveSubsetAlloc(image, border, threaded); // Subimage to convolve
-    psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Subimage mask
+    psImage *subImage = image ? psImageSubset(image, border) : NULL; // Subimage to convolve
+    psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Subimage mask
 
     psImage *convolved = psImageConvolveFFT(NULL, subImage, subMask, maskVal, kernel); // Convolution
 
-    convolveSubsetFree(image, subImage, threaded);
-    convolveSubsetFree(mask, subMask, threaded);
+    psFree(subImage);
+    psFree(subMask);
 
     // Now, we have to stick it in where it belongs
@@ -300 +261 @@
                                   region.y0 - size, region.y1 + size); // Add a border
 
-    bool threaded = pmSubtractionThreaded(); // Are we running threaded?
-
-    psImage *subVariance = convolveSubsetAlloc(variance, border, threaded); // Variance map
-    psImage *subSys = convolveSubsetAlloc(sys, border, threaded); // Systematic error image
-    psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Mask
+    psImage *subVariance = variance ? psImageSubset(variance, border) : NULL; // Variance map
+    psImage *subSys = sys ? psImageSubset(sys, border) : NULL; // Systematic error image
+    psImage *subMask = mask ? psImageSubset(mask, border) : NULL; // Mask
 
     // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once
@@ -310 +269 @@
     psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys
 
-    convolveSubsetFree(variance, subVariance, threaded);
-    convolveSubsetFree(sys, subSys, threaded);
-    convolveSubsetFree(mask, subMask, threaded);
+    psFree(subVariance);
+    psFree(subSys);
+    psFree(subMask);
 
     // Now, we have to stick it in where it belongs
@@ -420 +379 @@
         if (box > 0) {
             int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds
-            bool threaded = pmSubtractionThreaded(); // Are we running threaded?
-
             psRegion region = psRegionSet(colMin - box, colMax + box,
                                           rowMin - box, rowMax + box); // Region to convolve
 
-            psImage *image = convolveSubsetAlloc(subMask, region, threaded); // Mask to convolve
+            psImage *image = subMask ? psImageSubset(subMask, region) : NULL; // Mask to convolve
 
             psImage *convolved = psImageConvolveMask(NULL, image, subBad, subConvBad,
                                                      -box, box, -box, box); // Convolved subtraction mask
 
-            convolveSubsetFree(subMask, image, threaded);
+            psFree(image);
 
             psAssert(convolved->numCols - 2 * box == colMax - colMin, "Bad number of columns");
@@ -653 +610 @@
                   float value = 0.0;    // Value of convolved pixel
                   int uMin = x - size, uMax = x + size; // Range for u
-                  psF32 *xKernelData = xKernel->data.F32; // Kernel values
+                  psF32 *xKernelData = &xKernel->data.F32[xKernel->n - 1]; // Kernel values
                   psF32 *imageData = &image->kernel[y][uMin]; // Image values
-                  for (int u = uMin; u <= uMax; u++, xKernelData++, imageData++) {
+                  for (int u = uMin; u <= uMax; u++, xKernelData--, imageData++) {
                       value += *xKernelData * *imageData;
                   }
@@ -668 +625 @@
                   float value = 0.0;    // Value of convolved pixel
                   int vMin = y - size, vMax = y + size; // Range for v
-                  psF32 *yKernelData = yKernel->data.F32; // Kernel values
+                  psF32 *yKernelData = &yKernel->data.F32[yKernel->n - 1]; // Kernel values
                   psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way!
-                  for (int v = vMin; v <= vMax; v++, yKernelData++, imageData++) {
+                  for (int v = vMin; v <= vMax; v++, yKernelData--, imageData++) {
                       value += *yKernelData * *imageData;
                   }
@@ -821 +778 @@
     }
     psFree(mask);
+
+    // XXX raise an error?
+    if (isnan(stats->sampleMean)) {
+        psFree(stats);
+        return -1;
+    }
 
     double mean, rms;                 // Mean and RMS of deviations