Changeset 26479

Timestamp:

Dec 23, 2009, 8:08:30 AM (16 years ago)

Author:

eugene

Message:

adding deconvolved hermite functions

Location:

branches/eam_branches/20091201/psModules/src/imcombine

Files:

• Makefile.am (modified) (2 diffs)
• pmSubtractionDeconvolve.c (added)
• pmSubtractionDeconvolve.h (added)
• pmSubtractionKernels.c (modified) (1 diff)
• pmSubtractionKernels.h (modified) (4 diffs)

branches/eam_branches/20091201/psModules/src/imcombine/Makefile.am

@@ -14 +14 @@
         pmSubtractionKernels.c  \
         pmSubtractionHermitian.c        \
+        pmSubtractionDeconvolve.c       \
         pmSubtractionMask.c     \
         pmSubtractionMatch.c    \
@@ -34 +35 @@
         pmSubtractionKernels.h  \
         pmSubtractionHermitian.h        \
+        pmSubtractionDeconvolve.h       \
         pmSubtractionMask.h     \
         pmSubtractionMatch.h    \

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c

@@ -346 +346 @@
     return kernels;
 }
+
+# if (0)
+pmSubtractionKernels *pmSubtractionKernelsDECON_HERM(int size, int spatialOrder,
+                                                     const psVector *fwhmsIN, const psVector *ordersIN,
+                                                     float penalty, pmSubtractionMode mode)
+{
+    PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(fwhmsIN, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(ordersIN, NULL);
+    PS_ASSERT_VECTOR_TYPE(ordersIN, PS_TYPE_S32, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(fwhmsIN, ordersIN, NULL);
+    PS_ASSERT_INT_POSITIVE(size, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
+
+    // check the requested fwhm values: any values <= 0.0 should be dropped
+    psVector *fwhms  = psVectorAllocEmpty (fwhmsIN->n, PS_TYPE_F32);
+    psVector *orders = psVectorAllocEmpty (ordersIN->n, PS_TYPE_S32);
+    for (int i = 0; i < fwhmsIN->n; i++) {
+        if (fwhmsIN->data.F32[i] <= FLT_EPSILON) continue;
+        psVectorAppend(fwhms, fwhmsIN->data.F32[i]);
+        psVectorAppend(orders, ordersIN->data.S32[i]);
+    }
+
+    int numGaussians = fwhms->n;       // Number of Gaussians
+
+    int num = 0;                        // Number of basis functions
+    psString params = NULL;             // List of parameters
+    for (int i = 0; i < numGaussians; i++) {
+        int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
+        psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
+        num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
+    }
+
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECON_HERM, size, spatialOrder, penalty, mode); // The kernels
+    psStringAppend(&kernels->description, "DECON_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "DECONVOLVED HERM kernel: %s,%d --> %d elements", params, spatialOrder, num);
+    psFree(params);
+
+    // XXXXX hard-wired reference sigma for now of 1.7 pix
+    // generate the Gaussian deconvolution kernel
+    # define DECON_SIGMA 1.7
+    psKernel *kernelGauss = pmSubtractionDeconvolveGauss (size, DECON_SIGMA);
+
+    // Set the kernel parameters
+    int fullSize = 2 * size + 1;        // Full size of kernels
+    for (int i = 0, index = 0; i < numGaussians; i++) {
+        float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
+        // Iterate over (u,v) order
+        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(3); // Array to hold precalculated values
+                psVector *xKernel = preCalc->data[0] = subtractionKernelHERM(sigma, uOrder, size); // x Kernel
+                psVector *yKernel = preCalc->data[1] = subtractionKernelHERM(sigma, vOrder, size); // y Kernel
+                psKernel *kernelTarget = psKernelAlloc(-size, size, -size, size);       // Kernel
+
+                // generate 2D kernel, calculate moments
+                for (int v = -size, y = 0; v <= size; v++, y++) {
+                    for (int u = -size, x = 0; u <= size; u++, x++) {
+                        double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
+                        kernelTarget->kernel[v][u] = value;
+                    }
+                }
+
+                // deconvolve the target by the gaussian:
+                psKernel *kernel = pmSubtractionDeconvolveKernel(kernelTarget, kernelGauss); // Kernel
+                preCalc->data[2] = kernel;
+
+                psImage *kernelConv = psImageConvolveFFT(NULL, kernel, NULL, 0, kernelGauss);
+                pmSubtractionVisualShowSubtraction (kernelTarget->image, kernel->image, kernelConv->image);
+
+                // Normalise sum of kernel component to unity for even functions
+                if (uOrder % 2 == 0 && vOrder % 2 == 0) {
+                    double sum = 0.0;   // Sum of kernel component
+                    for (int v = 0; v < fullSize; v++) {
+                        for (int u = 0; u < fullSize; u++) {
+                            sum += xKernel->data.F32[u] * yKernel->data.F32[v];
+                        }
+                    }
+                    sum = 1.0 / sqrt(sum);
+                    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));
+
+#if 1
+                    fprintf(stderr, "%d norm: %e, null: %e\n", index, sum, kernel->kernel[0][0]);
+#endif
+                    
+                    kernel->kernel[0][0] -= 1.0;
+                    moment *= PS_SQR(sum);
+                }
+
+
+#if 1
+                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: %e\n", index, sum);
+#endif
+
+                kernels->widths->data.F32[index] = fwhms->data.F32[i];
+                kernels->u->data.S32[index] = uOrder;
+                kernels->v->data.S32[index] = vOrder;
+                if (kernels->preCalc->data[index]) {
+                    psFree(kernels->preCalc->data[index]);
+                }
+                kernels->preCalc->data[index] = preCalc;
+                kernels->penalties->data.F32[index] = kernels->penalty * fabsf(moment);
+
+                psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d %f\n", index,
+                        fwhms->data.F32[i], uOrder, vOrder, fabsf(moment));
+            }
+        }
+    }
+
+    return kernels;
+}
+# endif
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h

@@ -11 +11 @@
     PM_SUBTRACTION_KERNEL_ISIS,         ///< Traditional kernel --- gaussians modified by polynomials
     PM_SUBTRACTION_KERNEL_HERM,         ///< Hermitian polynomial kernels
+    PM_SUBTRACTION_KERNEL_DECON_HERM,   ///< Deconvolved Hermitian polynomial kernels
     PM_SUBTRACTION_KERNEL_SPAM,         ///< Summed Pixels for Advanced Matching --- summed delta functions
     PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
@@ -31 +32 @@
     psString description;               ///< Description of the kernel parameters
     long num;                           ///< Number of kernel components (not including the spatial ones)
-    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS or HERM)
-    psVector *widths;                   ///< Gaussian FWHMs (ISIS or HERM)
+    psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECON_HERM)
+    psVector *widths;                   ///< Gaussian FWHMs (ISIS, HERM or DECON_HERM)
     psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
-    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS or HERM)
+    psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS, HERM or DECON_HERM)
     float penalty;                      ///< Penalty for wideness
     psVector *penalties;                ///< Penalty for each kernel component
@@ -70 +71 @@
         PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
     } \
+    if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_DECON_HERM) { \
+        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
+        PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
+        PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
+    } \
     if ((KERNELS)->uStop || (KERNELS)->vStop) { \
         PS_ASSERT_VECTOR_NON_NULL((KERNELS)->uStop, RETURNVALUE); \
@@ -157 +163 @@
                                                );
 
+/// Generate DECON_HERM kernels
+pmSubtractionKernels *pmSubtractionKernelsDECON_HERM(int size, ///< Half-size of the kernel
+                                                     int spatialOrder, ///< Order of spatial variations
+                                                     const psVector *fwhms, ///< Gaussian FWHMs
+                                                     const psVector *orders, ///< order of hermitian polynomials
+                                                     float penalty, ///< Penalty for wideness
+                                                     pmSubtractionMode mode ///< Mode for subtraction
+    );
+
 /// Generate SPAM kernels
 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel