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

Timestamp:

Nov 2, 2009, 10:38:23 AM (17 years ago)

Author:

Paul Price

Message:

Merging PSF-matching code from branches/pap/ for dual convolution.

Location:

trunk/psModules/src/imcombine

Files:

• . (modified) (1 prop)
• pmSubtractionKernels.c (modified) (8 diffs)

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -81 +81 @@
     kernels->penalties = psVectorRealloc(kernels->penalties, start + numNew);
     kernels->inner = start;
+    kernels->num += numNew;
 
     // Generate a set of kernels for each (u,v)
@@ -94 +95 @@
             kernels->v->data.S32[index] = v;
             kernels->preCalc->data[index] = NULL;
-            kernels->penalties->data.F32[index] = kernels->penalty * (PS_SQR(u) + PS_SQR(v));
+            kernels->penalties->data.F32[index] = kernels->penalty * PS_SQR(PS_SQR(u) + PS_SQR(v));
 
             psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v);
         }
     }
+
+    kernels->widths->n = start + numNew;
+    kernels->u->n = start + numNew;
+    kernels->v->n = start + numNew;
+    kernels->preCalc->n = start + numNew;
+    kernels->penalties->n = start + numNew;
 
     return true;
@@ -140 +147 @@
         for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
             for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
-                psArray *preCalc = psArrayAlloc(2); // Array to hold precalculated values
+                psArray *preCalc = psArrayAlloc(3); // Array to hold precalculated values
                 psVector *xKernel = preCalc->data[0] = subtractionKernelISIS(sigma, uOrder, size); // x Kernel
                 psVector *yKernel = preCalc->data[1] = subtractionKernelISIS(sigma, vOrder, size); // y Kernel
+                psKernel *kernel = preCalc->data[2] = psKernelAlloc(-size, size, -size, size);      // Kernel
 
                 // Calculate moments
@@ -149 +157 @@
                     for (int u = -size, x = 0; u <= size; u++, x++) {
                         double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
-                        moment += value * (PS_SQR(u) + PS_SQR(v));
+                        kernel->kernel[v][u] = value;
+                        moment += value * PS_SQR((PS_SQR(u) + PS_SQR(v)));
                     }
                 }
@@ -164 +173 @@
                     psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
                     psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
+                    psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(PS_SQR(sum), PS_TYPE_F32));
+                    kernel->kernel[0][0] -= 1.0;
                     moment *= PS_SQR(sum);
                 }
+
+
+#if 0
+                double sum = 0.0;   // Sum of kernel component
+                for (int v = -size; v <= size; v++) {
+                    for (int u = -size; u <= size; u++) {
+                        sum += kernel->kernel[v][u];
+                    }
+                }
+                fprintf(stderr, "%d sum: %lf\n", index, sum);
+#endif
 
                 kernels->widths->data.F32[index] = fwhms->data.F32[i];
@@ -456 +478 @@
     PS_ASSERT_INT_LESS_THAN(inner, size, NULL);
 
-    // XXX GUNK doesn't seem to work --- doesn't add the POIS components, or at least, they're not noticed
-
     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
                                                                   penalty, mode); // Kernels
+    kernels->type = PM_SUBTRACTION_KERNEL_GUNK;
     psStringPrepend(&kernels->description, "GUNK=");
     psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder);
@@ -577 +598 @@
                                     poly->data.F32[j] = polyVal;
                                     norm += polyVal;
-                                    moment += polyVal * (PS_SQR(u) + PS_SQR(v));
+                                    moment += polyVal * PS_SQR(PS_SQR(u) + PS_SQR(v));
 
                                     psVectorExtend(uCoords, RINGS_BUFFER, 1);
@@ -603 +624 @@
                         psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32));
                     }
-//                    moment /= norm;
+                    moment /= norm;
                 }