Changeset 13389 for trunk/psModules/src/imcombine/pmSubtraction.c

Timestamp:

May 15, 2007, 2:01:03 PM (19 years ago)

Author:

Paul Price

Message:

Fixing subtraction of nominated kernel. This is used to maintain photometric constancy across the image when using a spatially variable kernel (see Alard 2000). Works for POIS kernels, but not yet for ISIS.

File:

• trunk/psModules/src/imcombine/pmSubtraction.c (modified) (5 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.6 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-15 04:33:58 $
+ *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-05-16 00:01:03 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -66 +66 @@
     assert(kernels->xOrder->n == numKernels);
     assert(kernels->yOrder->n == numKernels);
+
+    // Ensure the subIndex for POIS kernels is what is expected
+    assert(kernels->type != PM_SUBTRACTION_KERNEL_POIS ||
+           (kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0 &&
+            kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
+            kernels->yOrder->data.S32[kernels->subIndex] == 0));
 
     int size = kernels->size;           // Kernel half-size
@@ -82 +88 @@
               int v = kernels->v->data.S32[i]; // Offset in y
               kernel->kernel[v][u] += weightFunc(polyValue * solution->data.F64[i]);
+              if (i != kernels->subIndex) {
+                  // The (0,0) element is subtracted from most kernels to preserve photometric scaling
+                  kernel->kernel[0][0] -= weightFunc(solution->data.F64[i]);
+              }
               break;
           }
@@ -119 +129 @@
           int u = kernels->u->data.S32[index]; // Offset in x
           int v = kernels->v->data.S32[index]; // Offset in y
-          return weightFunc(polyValue) * image->data.F32[y + v][x + u];
+          double value = weightFunc(polyValue) * image->data.F32[y + v][x + u]; // Value of convolution
+          if (index != kernels->subIndex) {
+              // The (0,0) element is subtracted from most kernels to preserve photometric scaling
+              value -= image->data.F32[y][x];
+          }
+          return value;
       }
       case PM_SUBTRACTION_KERNEL_ISIS: {
@@ -243 +258 @@
                         convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues,
                                                                   imageWeighting);
-                    }
-
-                    if (spatialOrder > 0) {
-                        // Subtract a nominated convolution, to force flux conservation
-                        int subIndex = kernels->subIndex; // Index of convolution to subtract
-                        for (int i = 0; i < subIndex; i++) {
-                            convolutions->data.F64[i] -= convolutions->data.F64[subIndex];
-                        }
-                        for (int i = subIndex + 1; i < numKernels; i++) {
-                            convolutions->data.F64[i] -= convolutions->data.F64[subIndex];
-                        }
                     }