Changeset 14738 for trunk/psModules/src/imcombine/pmSubtractionParams.c

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).

File:

• trunk/psModules/src/imcombine/pmSubtractionParams.c (modified) (5 diffs)

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