Changeset 15756 for trunk/psModules/src/imcombine/pmSubtractionMatch.c

Timestamp:

Dec 6, 2007, 3:57:15 PM (18 years ago)

Author:

Paul Price

Message:

Implementing dual-convolution. This required reworking those
functions involved with calculating and solving the equations; moved
these into a separate file and made various other cleanups (e.g.,
assertions for pmSubtractionKernels). Removed the normalisation
(central pixel) term from all kernels, because this shouldn't be in
the second kernel (there's only one normalisation term between the two
kernels); the normalisation is treated explicitly in the equations,
along with the background (still only a constant background is
supported for now, but there's a lot of support for a polynomial
background for when I get around to putting it in the equations).
Tested the single convolution, and it's working; same results as
before, apparently. Haven't tested dual convolution yet.

File:

• trunk/psModules/src/imcombine/pmSubtractionMatch.c (modified) (19 diffs)

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -13 +13 @@
 #include "pmSubtractionKernels.h"
 #include "pmSubtractionStamps.h"
+#include "pmSubtractionEquation.h"
 #include "pmSubtraction.h"
+#include "pmSubtractionMask.h"
 #include "pmSubtractionMatch.h"
 
+#define TOL 1.0e-3                      // Tolerance for subtraction solution
+#define KERNEL_MOSAIC 2                 // Half-number of kernel instances in the mosaic image
 
 static bool useFFT = true;              // Do convolutions using FFT
@@ -63 +67 @@
 {
     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
-    *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, stampSpacing, mode);
+    *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, footprint,
+                                      stampSpacing, mode);
     if (!*stamps) {
         psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
@@ -72 +77 @@
 
     psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
-    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, footprint, size)) {
+    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, size)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
         return false;
@@ -86 +91 @@
 
 
-bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *ro1, const pmReadout *ro2,
+bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
                         int footprint, float regionSize, float stampSpacing, float threshold,
                         const psArray *sources, const char *stampsName,
@@ -95 +100 @@
                         psMaskType maskBlank, float badFrac, pmSubtractionMode mode)
 {
-    PS_ASSERT_PTR_NON_NULL(convolved, false);
+    PS_ASSERT_PTR_NON_NULL(conv1, false);
+    if (mode == PM_SUBTRACTION_MODE_DUAL) {
+        PS_ASSERT_PTR_NON_NULL(conv2, false);
+    }
     PS_ASSERT_PTR_NON_NULL(ro1, false);
     PS_ASSERT_IMAGE_NON_NULL(ro1->image, false);
@@ -174 +182 @@
 
     // Reset the output readout, just in case
-    if (convolved->image) {
-        psFree(convolved->image);
-        convolved->image = NULL;
-    }
-    if (convolved->mask) {
-        psFree(convolved->mask);
-        convolved->mask = NULL;
-    }
-    if (convolved->weight) {
-        psFree(convolved->weight);
-        convolved->weight = NULL;
+    if (conv1->image) {
+        psFree(conv1->image);
+        conv1->image = NULL;
+    }
+    if (conv1->mask) {
+        psFree(conv1->mask);
+        conv1->mask = NULL;
+    }
+    if (conv1->weight) {
+        psFree(conv1->weight);
+        conv1->weight = NULL;
     }
 
@@ -250 +258 @@
 
             if (sources) {
-                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, mode);
+                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, footprint,
+                                                           stampSpacing, mode);
             } else if (stampsName && strlen(stampsName) > 0) {
-                stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, stampSpacing, mode);
+                stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, footprint,
+                                                        stampSpacing, mode);
             }
 
@@ -295 +305 @@
                 // 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",
+                                                       inner, binning, ringsOrder, mode);
+            }
+            psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL",
                              PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
 
@@ -312 +322 @@
 
                 psTrace("psModules.imcombine", 3, "Calculating equation...\n");
-                if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) {
+                if (!pmSubtractionCalculateEquation(stamps, kernels)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
@@ -320 +330 @@
 
                 psTrace("psModules.imcombine", 3, "Solving equation...\n");
-                solution = pmSubtractionSolveEquation(solution, stamps);
-                if (!solution) {
+
+                if (!pmSubtractionSolveEquation(kernels, stamps, TOL)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
@@ -328 +338 @@
                 memCheck("  solve equation");
 
-                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
-                                                                        mode); // Deviations for each stamp
+                psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
                 if (!deviations) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
@@ -351 +360 @@
             if (numRejected > 0) {
                 psTrace("psModules.imcombine", 3, "Solving equation...\n");
-                solution = pmSubtractionSolveEquation(solution, stamps);
-                if (!solution) {
+                if (!pmSubtractionSolveEquation(kernels, stamps, TOL)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
                 }
-                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
-                                                                        mode); // Deviations for each stamp
+                psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
                 if (!deviations) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
@@ -376 +383 @@
                 psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32);
                 psImageInit(convKernels, NAN);
-                for (int j = -2; j <= 2; j++) {
-                    for (int i = -2; i <= 2; i++) {
-                        psImage *kernel = pmSubtractionKernelImage(solution, kernels, (float)i / 2.0,
-                                                                   (float)j / 2.0); // Image of the kernel
+                for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) {
+                    for (int i = -KERNEL_MOSAIC; i <= KERNEL_MOSAIC; i++) {
+                        psImage *kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC,
+                                                                   (float)j / (float)KERNEL_MOSAIC,
+                                                                   false); // Image of the kernel
                         if (!kernel) {
                             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
@@ -386 +394 @@
                         }
 
-                        if (psImageOverlaySection(convKernels, kernel, (i + 2) * fullSize,
-                                                  (j + 2) * fullSize, "=") == 0) {
+                        if (psImageOverlaySection(convKernels, kernel, (i + KERNEL_MOSAIC) * fullSize,
+                                                  (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
                             psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
                             psFree(kernel);
@@ -399 +407 @@
                 psString comment = NULL; // Comment for metadata
                 psStringAppend(&comment, "Subtraction kernel for region %s", regionString);
-                psMetadataAddImage(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE",
+                psMetadataAddImage(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE",
                                    PS_META_DUPLICATE_OK, comment, convKernels);
                 psFree(comment);
@@ -427 +435 @@
 
             psTrace("psModules.imcombine", 2, "Convolving...\n");
-            if (!pmSubtractionConvolve(convolved, ro1, ro2, subMask, maskBlank, region,
-                                       solution, kernels, mode, useFFT)) {
+            if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBlank, region, kernels, useFFT)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
                 goto MATCH_ERROR;
@@ -439 +446 @@
                 psString comment = NULL; // Comment for metadata
                 psStringAppend(&comment, "Subtraction solution for region %s", regionString);
-                psMetadataAddVector(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.SOLUTION",
+                psMetadataAddVector(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.SOLUTION",
                                     PS_META_DUPLICATE_OK, comment, solution);
                 psFree(comment);
                 if (region) {
-                    psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
+                    psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
                                      PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region);
                 } else {
                     region = psRegionAlloc(0, numCols, 0, numRows);
-                    psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
+                    psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
                                      PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region);
                     psFree(region);
@@ -458 +465 @@
 
             // There is data in the readout now
-            convolved->data_exists = true;
-            if (convolved->parent) {
-                convolved->parent->data_exists = true;
-                convolved->parent->parent->data_exists = true;
+            conv1->data_exists = true;
+            if (conv1->parent) {
+                conv1->parent->data_exists = true;
+                conv1->parent->parent->data_exists = true;
+            }
+            if (mode == PM_SUBTRACTION_MODE_DUAL) {
+                conv2->data_exists = true;
+                if (conv2->parent) {
+                    conv2->parent->data_exists = true;
+                    conv2->parent->parent->data_exists = true;
+                }
             }
         }
@@ -474 +488 @@
     weight = NULL;
 
-    if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) {
+    if (!pmSubtractionBorder(conv1->image, conv1->weight, conv1->mask, size, maskBlank)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
         goto MATCH_ERROR;