Changeset 14701 for trunk/psModules/src/imcombine/pmSubtractionKernels.c

Timestamp:

Aug 29, 2007, 5:50:28 PM (19 years ago)

Author:

Paul Price

Message:

Extensive changes to image subtraction codes. Discovered that I
wasn't (again) implementing the algorithm correctly --- I was
accumulating

sum_x,y C_i(x,y) sum_x,y C_j(x,y)

instead of

sum_x,y C_i(x,y) C_j(x,y)

i.e., I was doing the sums separately. I'm amazed that it was
actually getting fairly decent results before, but now it's doing it
properly. As part of this fix, added in, for each stamp, images
corresponding to the convolution of the reference with each of the
kernel components. This allows fast rejection (instead of generating
the kernel and convolving, can simply sum the individual convolutions
with the correct scaling).
Added use of FFT to do convolutions. This involved changing around
the convolvePixel function (now convolveRef, since it generates the
entire image for the convolution of the reference stamp with a kernel
component, rather than a single pixel as previously). Tested with
POIS and ISIS kernels, but no others yet. FFT is nice and fast for
the large (size = 10) kernels I've been using, so I also moved it in
to the convolveRef function for appropriate kernel types.
As part of the FFT move, discovered that my convolutions were around
the wrong way --- supposed to be convolved = image[y][x] * conv[y - v][x - u]
instead of conv[y + v][x + u]. However, this problem didn't manifest
itself previously because everything was consistently around the wrong
way. The FFT does it the Correct Way, so it revealed the problem. This
is now fixed throughout the image subtraction code.

File:

• trunk/psModules/src/imcombine/pmSubtractionKernels.c (modified) (2 diffs)

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -216 +216 @@
 
     // Subtract a reference kernel from all the others to maintain flux scaling
-    int subIndex = 0;                   // Index of kernel to subtract from others
-    psKernel *subtract = kernels->preCalc->data[subIndex]; // Kernel to subtract from others
-    int numGaussians = fwhms->n;       // Number of Gaussians
-    for (int i = 0, index = 0; i < numGaussians; i++) {
-        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
-            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
-                if (uOrder % 2 == 0 && vOrder % 2 == 0 && index != subIndex) {
-                    psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract
-                    psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
+//    if (spatialOrder > 0) {
+        int subIndex = 0;                   // Index of kernel to subtract from others
+        psKernel *subtract = kernels->preCalc->data[subIndex]; // Kernel to subtract from others
+        int numGaussians = fwhms->n;       // Number of Gaussians
+        for (int i = 0, index = 0; i < numGaussians; i++) {
+            for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
+                for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
+                    if (uOrder % 2 == 0 && vOrder % 2 == 0 && index != subIndex) {
+                        psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract
+                        psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
+                    }
                 }
             }
         }
-    }
+//    }
 
     if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) {
@@ -452 +454 @@
     // Subtract unity from the kernels to maintain photometric flux scaling
     int numGaussians = fwhms->n;       // Number of Gaussians
-    for (int i = 0, index = 0; i < numGaussians; i++) {
-        for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
-            for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
-                if (uOrder % 2 == 0 && vOrder % 2 == 0) {
-                    psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract
-                    kernel->kernel[0][0] -= 1.0;
+//    if (spatialOrder > 0) {
+        for (int i = 0, index = 0; i < numGaussians; i++) {
+            for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
+                for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
+                    if (uOrder % 2 == 0 && vOrder % 2 == 0) {
+                        psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract
+                        kernel->kernel[0][0] -= 1.0;
+                    }
                 }
             }
         }
-    }
+//    }
 
     int numISIS = kernels->num;         // Number of ISIS kernels
     int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements
 
-    if (!p_pmSubtractionKernelsAddGrid(kernels, numGaussians, inner)) {
+    if (!p_pmSubtractionKernelsAddGrid(kernels, numISIS, inner)) {
         psAbort("Should never get here.");
     }