Changeset 14738

Timestamp:

Sep 4, 2007, 1:56:07 PM (19 years ago)

Author:

Paul Price

Message:

Consolidating some code. Making it so that the convolutions of the reference don't have to be performed twice if the optimum kernel is computed (before, the convolutions were calculated twice --- once in getting the optimum kernel set, and once in calculating the least-squares matrix).

Location:

trunk/psModules/src/imcombine

Files:

• pmSubtraction.c (modified) (4 diffs)
• pmSubtraction.h (modified) (3 diffs)
• pmSubtractionMatch.c (modified) (7 diffs)
• pmSubtractionParams.c (modified) (5 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-30 23:09:52 $
+ *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-04 23:56:07 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -548 +548 @@
 }
 
+bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint)
+{
+    PS_ASSERT_PTR_NON_NULL(stamp, false);
+    PS_ASSERT_KERNEL_NON_NULL(stamp->reference, false);
+    PS_ASSERT_PTR_NON_NULL(kernels, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
+    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
+
+    if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Stamp not marked for calculation.");
+        return false;
+    }
+
+    if (stamp->convolutions) {
+        // Already done
+        return true;
+    }
+
+    stamp->convolutions = psArrayAlloc(kernels->num);
+    psKernel *reference = stamp->reference; // Reference postage stamp
+
+    for (int i = 0; i < kernels->num; i++) {
+        stamp->convolutions->data[i] = convolveRef(kernels, i, reference, footprint);
+    }
+
+    return true;
+}
+
 
 bool pmSubtractionCalculateEquation(psArray *stamps, const pmSubtractionKernels *kernels, int footprint)
@@ -591 +619 @@
         psVectorInit(stampVector, 0.0);
 
-        psKernel *reference = stamp->reference; // Reference postage stamp
         psKernel *input = stamp->input; // Input postage stamp
         psKernel *weight = stamp->weight; // Weight map postage stamp
 
         // Generate convolutions of the reference
+        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
+            return NULL;
+        }
+
         psArray *convolutions = stamp->convolutions; // Convolutions of the reference for each kernel
-        if (!convolutions) {
-            stamp->convolutions = convolutions = psArrayAlloc(numKernels);
-        }
-        for (int j = 0; j < numKernels; j++) {
-            if (convolutions->data[j]) {
-                psFree(convolutions->data[j]);
-            }
-            convolutions->data[j] = convolveRef(kernels, j, reference, footprint);
-#ifdef TESTING
-            {
-                psKernel *conv = convolutions->data[j]; // Convolution of interest
-                psString filename = NULL;
-                psStringAppend(&filename, "conv_%03d_%03d.fits", i, j);
-                psFits *fits = psFitsOpen(filename, "w");
-                psFree(filename);
-                psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
-                psFitsClose(fits);
-            }
-#endif
-       }
-
         psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms
 
@@ -900 +911 @@
             stamp->y = 0;
             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
+            // Recalculate convolutions
+            psFree(stamp->convolutions);
+            stamp->convolutions = NULL;
         }
     }

trunk/psModules/src/imcombine/pmSubtraction.h

@@ -6 +6 @@
  * @author GLG, MHPCC
  *
- * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-08-30 21:45:26 $
+ * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-09-04 23:56:07 $
  *
  * Copyright 2004-207 Institute for Astronomy, University of Hawaii
@@ -17 +17 @@
 #include <pslib.h>
 #include "pmSubtractionKernels.h"
+#include "pmSubtractionStamps.h"
 
 /// @addtogroup imcombine Image Combinations
@@ -37 +38 @@
                            int size, ///< Half-size of the kernel (pmSubtractionKernels.size)
                            int footprint ///< Half-size of the kernel footprint
+    );
+
+/// Convolve the reference stamp with the kernel components
+bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, ///< Stamp to convolve
+                                const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                int footprint ///< Half-size of region over which to calculate equation
     );
 

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -214 +214 @@
     // Putting important variable declarations here, since they are freed after a "goto" if there is an error.
     psImage *subMask = NULL;            // Mask for subtraction
-    psArray *stamps = NULL;             // Stamps for matching PSF
     psRegion *region = NULL;            // Iso-kernel region
     psString regionString = NULL;       // String for region
+    psArray *stamps = NULL;             // Stamps for matching PSF
     psVector *solution = NULL;          // Solution to match PSF
     pmSubtractionKernels *kernels = NULL; // Kernel basis functions
@@ -224 +224 @@
     if (stampsName && strlen(stampsName) > 0) {
         psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName);
-        psVector *xStamp = NULL, *yStamp = NULL; // Stamp positions and fluxes
         if (input) {
             // We have x, y because the target is provided by the input image
@@ -234 +233 @@
         if (!stampsData) {
             psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName);
-            goto ERROR;
+            return false;
         }
 
         // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
-        xStamp = stampsData->data[0];
-        yStamp = stampsData->data[1];
+        psVector *xStamp = stampsData->data[0]; // x positions
+        psVector *yStamp = stampsData->data[1]; // y positions
         psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
         psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
@@ -251 +250 @@
 
     memCheck("mask");
-
-    if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
-        if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, NULL,
-                       threshold, stampSpacing, targetWidth, size, footprint)) {
-            goto ERROR;
-        }
-        kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
-                                                  stamps, footprint, optThreshold);
-        // XXX This is not quite the best thing to do here--- we'll find the same set of stamps again later,
-        // but since we may not be interested in the entire image on the first pass through, we have to blow
-        // the whole lot away.  The alternative is to mark stamps that are outside the region of interest, and
-        // ignore them when generating sums and rejecting.
-        psFree(stamps);
-        stamps = NULL;
-
-        if (!kernels) {
-            psErrorClear();
-            psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
-        }
-    }
-
-    if (kernels == NULL) {
-        // Not an ISIS/GUNK kernel, or the optimum kernel search failed
-        kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
-                                               inner, binning, ringsOrder);
-    }
-
-    psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL", PS_DATA_UNKNOWN,
-                     "Subtraction kernels", kernels);
-
-    memCheck("kernels");
 
     // Get region of interest
@@ -307 +275 @@
                         regionString, numCols, numRows);
             }
+
+            // Define kernel basis functions
+            if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
+                if (!getStamps(&stamps, stampsData, reference, input, subMask, weight, NULL,
+                               threshold, stampSpacing, targetWidth, size, footprint)) {
+                    goto ERROR;
+                }
+                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
+                                                          stamps, footprint, optThreshold);
+                if (!kernels) {
+                    psErrorClear();
+                    psWarning("Unable to derive optimum ISIS kernel --- switching to default.");
+                }
+            }
+            if (kernels == NULL) {
+                // Not an ISIS/GUNK kernel, or the optimum kernel search failed
+                kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
+                                                       inner, binning, ringsOrder);
+            }
+            psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL",
+                             PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
+
+            memCheck("kernels");
 
             int numRejected = -1;               // Number of rejected stamps in each iteration
@@ -422 +413 @@
                 goto ERROR;
             }
+            psFree(kernels);
+            kernels = NULL;
 
             // Put the solution on the metadata
@@ -466 +459 @@
         goto ERROR;
     }
-    psFree(kernels);
 
     memCheck("convolution");

trunk/psModules/src/imcombine/pmSubtractionParams.c

@@ -221 +221 @@
     }
 
-    // Generate the convolutions
-    for (int i = 0; i < numStamps; i++) {
-        pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
-        if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || stamp->status == PM_SUBTRACTION_STAMP_NONE) {
-            continue;
-        }
-        if (!stamp->convolutions) {
-            stamp->convolutions = psArrayAlloc(numKernels);
-        }
-
-        for (int j = 0; j < numKernels; j++) {
-            psKernel *kernel = kernels->preCalc->data[j]; // Precalculated kernel
-            stamp->convolutions->data[j] = p_pmSubtractionConvolveStampPrecalc(stamp->reference, kernel);
-        }
-    }
-
-    // Calculate sums invariant with the input image
+    // Generate the convolutions, accumulate sums, and measure initial chi^2
     double sum1 = 0.0;                  // sum of 1/sigma(x,y)^2
-    for (int i = 0; i < numStamps; i++) {
-        pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
-        if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED || stamp->status == PM_SUBTRACTION_STAMP_NONE) {
-            continue;
-        }
-        psKernel *weight = stamp->weight; // Weight map for stamp
-
-        for (int v = -footprint; v <= footprint; v++) {
-            psF32 *wt = &weight->kernel[v][-footprint]; // Dereference weight map
-            for (int u = -footprint; u <= footprint; u++, wt++) {
-                sum1 += 1.0 / *wt;
-            }
-        }
-        if (!isfinite(sum1)) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                    "Sum of 1/sigma^2 is non-finite for stamp %d (%d,%d)\n",
-                    i, (int)stamp->x, (int)stamp->y);
-            psFree(kernels);
-            return NULL;
-        }
-    }
     psVector *sumC = psVectorAlloc(numKernels, PS_TYPE_F64); // sum of R(x)*k(u)/sigma(x)^2
     psVector *sumCC = psVectorAlloc(numKernels, PS_TYPE_F64); // sum of [R(x)*k(u)]^2/sigma(x)^2
     psVectorInit(sumC, 0.0);
     psVectorInit(sumCC, 0.0);
-    for (int i = 0; i < numKernels; i++) {
-        for (int j = 0; j < numStamps; j++) {
-            pmSubtractionStamp *stamp = stamps->data[j]; // Stamp of interest
-            if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED ||
-                stamp->status == PM_SUBTRACTION_STAMP_NONE) {
-                continue;
-            }
-            accumulateConvolutions(&sumC->data.F64[i], &sumCC->data.F64[i], stamp, i, footprint);
-        }
-    }
-
-    // Initial chi^2
     double lastChi2 = 0.0;              // Chi^2 from last iteration
     int numPixels = 0;                  // Number of pixels contributing to chi^2
@@ -283 +234 @@
             continue;
         }
+        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
+            psFree(inputs);
+            psFree(kernels);
+            return NULL;
+        }
+
+        // This sum is invariant to the kernel
+        psKernel *weight = stamp->weight; // Weight map for stamp
+        for (int v = -footprint; v <= footprint; v++) {
+            psF32 *wt = &weight->kernel[v][-footprint]; // Dereference weight map
+            for (int u = -footprint; u <= footprint; u++, wt++) {
+                sum1 += 1.0 / *wt;
+            }
+        }
+        if (!isfinite(sum1)) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "Sum of 1/sigma^2 is non-finite for stamp %d (%d,%d)\n",
+                    i, (int)stamp->x, (int)stamp->y);
+            psFree(inputs);
+            psFree(kernels);
+            return NULL;
+        }
+
+        for (int j = 0; j < numKernels; j++) {
+            accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint);
+        }
+
         lastChi2 += initialChi2(inputs->data[i], stamp, footprint);
         numPixels += PS_SQR(2 * footprint + 1);
@@ -395 +374 @@
     psVector *widthsNew = psVectorAlloc(newSize, PS_TYPE_F32);
     psArray *preCalcNew = psArrayAlloc(newSize);
+    psArray *convNew = psArrayAlloc(numStamps);
+    for (int i = 0; i < numStamps; i++) {
+        convNew->data[i] = psArrayAlloc(newSize);
+    }
 
     for (int i = 0; i < numKernels; i++) {
@@ -403 +386 @@
             widthsNew->data.F32[rank] = kernels->widths->data.F32[i];
             preCalcNew->data[rank] = psMemIncrRefCounter(kernels->preCalc->data[i]);
+
+            for (int j = 0; j < numStamps; j++) {
+                psArray *convolutions = convNew->data[j]; // Convolutions for this stamp
+                pmSubtractionStamp *stamp = stamps->data[j]; // Stamp of interest
+                convolutions->data[rank] = psMemIncrRefCounter(stamp->convolutions->data[i]);
+            }
         }
     }
@@ -415 +404 @@
     kernels->num = newSize;
 
+    for (int i = 0; i < numStamps; i++) {
+        pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
+        psFree(stamp->convolutions);
+        stamp->convolutions = convNew->data[i];
+    }
+
     psFree(ranking);