Changeset 14420

Timestamp:

Aug 7, 2007, 9:02:55 AM (19 years ago)

Author:

Paul Price

Message:

Need to normalise sum of kernel components to unity in order to preserve photometric scaling.

Location:

trunk/psModules/src/imcombine

Files:

• pmSubtraction.c (modified) (7 diffs)
• pmSubtractionKernels.c (modified) (2 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-21 04:30:57 $
+ *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-07 19:02:55 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -124 +124 @@
             int vStart = (KERNELS)->v->data.S32[i]; \
             int vStop = (KERNELS)->vStop->data.S32[i]; \
+            /* Normalising sum of kernel component to unity */ \
+            value /= FUNC((uStop - uStart) * (vStop - vStart)); \
             for (int v = vStart; v <= vStop; v++) { \
                 for (int u = uStart; u <= uStop; u++) { \
@@ -180 +182 @@
               psVector *poly = preCalc->data[2]; /* Polynomial values */ \
               int num = uCoords->n;     /* Number of pixels */ \
+              value /= weightFunc(num); /* Normalising sum of kernel component to unity */ \
               for (int j = 0; j < num; j++) { \
                   int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; /* Kernel coordinates */ \
@@ -248 +251 @@
               int vStart = kernels->v->data.S32[i];
               int vStop = kernels->vStop->data.S32[i];
+
+              // Normalising sum of kernel component to unity
+              value /= weightFunc((uStop - uStart) * (vStop - vStart));
+
               for (int v = vStart; v <= vStop; v++) {
                   for (int u = uStart; u <= uStop; u++) {
@@ -304 +311 @@
               psVector *poly = preCalc->data[2]; // Polynomial values
               int num = uCoords->n;     // Number of pixels
-
+              value /= weightFunc(num); // Normalising sum of kernel component to unity
               for (int j = 0; j < num; j++) {
                   int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
@@ -362 +369 @@
               }
           }
+          sum /= (uStop - uStart) * (vStop - vStart); // Normalising sum of kernel component to unity
           if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
               // The (0,0) element is subtracted from most kernels to preserve photometric scaling
@@ -424 +432 @@
               sum += image->data.F32[y + v][x + u] * poly->data.F32[j];
           }
+          sum /= (double)num;           // Normalising sum of kernel component to unity
           // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
           return polyValue * sum - image->data.F32[0][0];

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -147 +147 @@
                 // Set the pre-calculated kernel
                 psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
+                double sum = 0.0;       // Normalisation
                 for (int v = -size; v <= size; v++) {
                     for (int u = -size; u <= size; u++) {
-                        preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *
+                        sum += preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *
                             expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i]));
                     }
                 }
+                // Normalise sum of kernel component to unity
+                psBinaryOp(preCalc->image, preCalc->image, "*", psScalarAlloc(1.0/sum, PS_TYPE_F32));
 
                 // Iterate over spatial order.  This loop creates the terms for
@@ -450 +453 @@
                 // Set the pre-calculated kernel
                 psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
+                double sum = 0.0;       // Normalisation
                 for (int v = -size; v <= size; v++) {
                     for (int u = -size; u <= size; u++) {
-                        preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *
+                        sum += preCalc->kernel[v][u] = power(u, uOrder) * power(v, vOrder) *
                             expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigmas->data.F32[i]));
                     }
                 }
+                // Normalise sum of kernel component to unity
+                psBinaryOp(preCalc->image, preCalc->image, "*", psScalarAlloc(1.0/sum, PS_TYPE_F32));
 
                 // Iterate over spatial order.  This loop creates the terms for