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

Timestamp:

Jul 10, 2007, 1:54:26 PM (19 years ago)

Author:

Paul Price

Message:

Adding additional convolution kernels: SPAM (binned pixels around a central single-pixel core) and FRIES (binning changes as a (Fibonacci-like) function of radius).

File:

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

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-09 01:06:48 $
+ *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-07-10 23:54:26 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -24 +24 @@
 
 #include "pmSubtraction.h"
+
+//#define CAN_PASS_STATIC_INLINE          // Can I pass the address of a static inline function?
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -53 +55 @@
 }
 
+#ifndef CAN_PASS_STATIC_INLINE
 // Generate an image of the solved kernel
 // Meant to replace the following function declaration:
@@ -61 +64 @@
 //                                      double (*weightFunc)(double value) // Function for weighting
 //     );
-// because even if the weightFunc can't be "static inline", they still appear in the assembly.
+// because even though the weightFunc is "static inline", it still appears in the assembly.
 // TARGET: Kernel, to return (psKernel*)
 // SOLUTION: Solution to the least-squares problem (psVector*)
@@ -75 +78 @@
     \
     /* Ensure the subIndex for POIS kernels is what is expected */ \
-    assert((KERNELS)->type != PM_SUBTRACTION_KERNEL_POIS || \
+    assert(((KERNELS)->type != PM_SUBTRACTION_KERNEL_POIS && \
+           (KERNELS)->type != PM_SUBTRACTION_KERNEL_SPAM && \
+           (KERNELS)->type != PM_SUBTRACTION_KERNEL_FRIES) || \
            ((KERNELS)->u->data.S32[(KERNELS)->subIndex] == 0 && \
             (KERNELS)->v->data.S32[(KERNELS)->subIndex] == 0 && \
@@ -103 +108 @@
               break; \
           } \
+          /* SPAM and FRIES use the same method */ \
+          case PM_SUBTRACTION_KERNEL_SPAM: \
+          case PM_SUBTRACTION_KERNEL_FRIES: { \
+            int uStart = (KERNELS)->u->data.S32[i]; \
+            int uStop = (KERNELS)->uStop->data.S32[i]; \
+            int vStart = (KERNELS)->v->data.S32[i]; \
+            int vStop = (KERNELS)->vStop->data.S32[i]; \
+            for (int v = vStart; v <= vStop; v++) { \
+                for (int u = uStart; u <= uStop; u++) { \
+                    (TARGET)->kernel[v][u] += value; \
+                } \
+            } \
+            if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \
+                /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \
+                (TARGET)->kernel[0][0] += subValue; \
+            } \
+            break; \
+          } \
           case PM_SUBTRACTION_KERNEL_ISIS: { \
               psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \
@@ -122 +145 @@
     } \
 }
-
-#if 0
+#else
 // Generate an image of the solved kernel
 static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return
@@ -139 +161 @@
 
     // Ensure the subIndex for POIS kernels is what is expected
-    assert(kernels->type != PM_SUBTRACTION_KERNEL_POIS ||
+    assert((kernels->type != PM_SUBTRACTION_KERNEL_POIS &&
+            kernels->type != PM_SUBTRACTION_KERNEL_SPAM &&
+            kernels->type != PM_SUBTRACTION_KERNEL_FRIES) ||
            (kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0 &&
             kernels->xOrder->data.S32[kernels->subIndex] == 0 &&
@@ -167 +191 @@
               break;
           }
+          /* SPAM and FRIES use the same method */
+          case PM_SUBTRACTION_KERNEL_SPAM:
+          case PM_SUBTRACTION_KERNEL_FRIES:
+            int uStart = kernels->u->data.S32[i];
+            int uStop = kernels->uStop->data.S32[i];
+            int vStart = kernels->v->data.S32[i];
+            int vStop = kernels->vStop->data.S32[i];
+            for (int v = vStart; v <= vStop; v++) {
+                for (int u = uStart; u <= uStop; u++) {
+                    kernel->kernel[v][u] += value;
+                }
+            }
+            if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
+                /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */
+                kernel->kernel[0][0] += subValue;
+            }
+            break;
           case PM_SUBTRACTION_KERNEL_ISIS: {
               psKernel *preCalc = kernels->preCalc->data[i];// Precalculated values
@@ -190 +231 @@
 #endif
 
+#ifndef CAN_PASS_STATIC_INLINE
 // Generate the convolved pixel value
 // Meant to replace the following function declaration:
@@ -199 +241 @@
 //                                    double (*weightFunc)(double value) // Function for weighting
 //     );
-// because even if the weightFunc can't be "static inline", they still appear in the assembly.
+// because even though weightFunc is "static inline", it still appears in the assembly.
 //
 // TARGET: The value to 'return' (double)
@@ -225 +267 @@
           break; \
       } \
-      case PM_SUBTRACTION_KERNEL_ISIS: { \
+      /* Method for SPAM and FRIES is the same */ \
+      case PM_SUBTRACTION_KERNEL_SPAM: \
+      case PM_SUBTRACTION_KERNEL_FRIES: { \
+          int uStart = (KERNELS)->u->data.S32[(INDEX)]; \
+          int uStop = (KERNELS)->uStop->data.S32[(INDEX)]; \
+          int vStart = (KERNELS)->v->data.S32[(INDEX)]; \
+          int vStop = (KERNELS)->vStop->data.S32[(INDEX)]; \
+          for (int v = vStart; v <= vStop; v++) { \
+              for (int u = uStart; u <= uStop; u++) { \
+                  (TARGET) += FUNC(polyValue) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \
+              } \
+          } \
+          if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \
+              /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */ \
+              (TARGET) += FUNC(-1.0) * (IMAGE)->data.F32[(Y)][(X)]; \
+          } \
+          break; \
+      } \
+            case PM_SUBTRACTION_KERNEL_ISIS: { \
           psKernel *kernel = (KERNELS)->preCalc->data[(INDEX)]; /* The convolution kernel */ \
           psKernel *subKernel = (KERNELS)->preCalc->data[(KERNELS)->subIndex]; /* Kernel to subtract */ \
@@ -231 +291 @@
           double sum = 0.0;             /* Accumulated sum from convolution */ \
           double sub = 0.0;             /* Accumulated sum to subtract */ \
-          for (int v = -size; v <= size; v++) { \
-              for (int u = -size; u <= size; u++) { \
-                  sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \
+          for (int v = -size, y = (Y) - size; v <= size; v++, y++) { \
+              for (int u = -size, x = (X) - size; u <= size; u++, x++) { \
+                  sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[y][x]; \
                   /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \
                   if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \
-                      sub += FUNC(subKernel->kernel[v][u]) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \
+                      sub += FUNC(subKernel->kernel[v][u]) * (IMAGE)->data.F32[y][x]; \
                   } \
               } \
@@ -247 +307 @@
     } \
 }
-
-#if 0
+#else
 // Generate the convolved pixel value
 static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions
@@ -273 +332 @@
           }
           return value;
+      }
+      /* Method for SPAM and FRIES is the same */
+      case PM_SUBTRACTION_KERNEL_SPAM:
+      case PM_SUBTRACTION_KERNEL_FRIES: {
+          int uStart = kernels->u->data.S32[index];
+          int uStop = kernels->uStop->data.S32[index];
+          int vStart = kernels->v->data.S32[index];
+          int vStop = kernels->vStop->data.S32[index];
+          for (int v = vStart; v <= vStop; v++) {
+              for (int u = uStart; u <= uStop; u++) {
+                  value += weightFunc(polyValue) * image->data.F32[y + v][x + u];
+              }
+          }
+          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
+              /* The (0,0) element is subtracted from most kernels to preserve photometric scaling */
+              value += weightFunc(-1.0) * image->data.F32[y][x];
+          }
+          break;
       }
       case PM_SUBTRACTION_KERNEL_ISIS: {
@@ -333 +410 @@
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-psImage *pmSubtractionMask(const psImage *inMask, const psImage *refMask, psMaskType maskVal,
-                           int size, int footprint)
-{
-    PS_ASSERT_IMAGE_NON_NULL(inMask, NULL);
-    PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, NULL);
-    PS_ASSERT_IMAGE_NON_NULL(refMask, NULL);
-    PS_ASSERT_IMAGE_TYPE(refMask, PS_TYPE_MASK, NULL);
-    PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, refMask, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(size, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);
-
-    // Size of the images
-    int numCols = inMask->numCols;
-    int numRows = inMask->numRows;
-
-    // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask
-    psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask
-    psImageInit(mask, 0);
-
-    // Dereference for convenience
-    psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA;
-    psMaskType **inData = inMask->data.PS_TYPE_MASK_DATA;
-    psMaskType **refData = refMask->data.PS_TYPE_MASK_DATA;
-
-    // Block out a border around the edge of the image
-
-    // Bottom stripe
-    for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) {
-        for (int x = 0; x < numCols; x++) {
-            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
-        }
-    }
-    // Either side
-    for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) {
-        for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) {
-            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
-        }
-        for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) {
-            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
-        }
-    }
-    // Top stripe
-    for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) {
-        for (int x = 0; x < numCols; x++) {
-            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
-        }
-    }
-
-    // Mask the bad pixels, and around them
-    for (int y = footprint + size; y < numRows - footprint - size; y++) {
-        for (int x = 0; x < numCols; x++) {
-            if (inData[y][x] & maskVal) {
-                maskData[y][x] |= PM_SUBTRACTION_MASK_INPUT;
-                // Block out the entire stamp footprint around this pixel
-                for (int v = y - footprint; v <= y + footprint; v++) {
-                    for (int u = x - footprint; u <= x + footprint; u++) {
-                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
-                    }
-                }
-            }
-            if (refData[y][x] & maskVal) {
-                maskData[y][x] |= PM_SUBTRACTION_MASK_REF;
-
-                // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with
-                // this bad pixel (within 'size').  Then we want to block out with the FOOTPRINT mask
-                // everything within a footprint's distance of those (within 'footprint').
-
-                // Bottom stripe
-                for (int v = PS_MAX(y - footprint - size, 0); v < y - size; v++) {
-                    for (int u = PS_MAX(x - footprint - size, 0);
-                         u <= PS_MIN(x + footprint + size, numCols - 1); u++) {
-                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
-                    }
-                }
-                // Middle
-                for (int v = PS_MAX(y - size, 0); v <= PS_MIN(y + size, numRows - 1); v++) {
-                    // Left side
-                    for (int u = PS_MAX(x - footprint - size, 0); u < x - size; u++) {
-                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
-                    }
-                    // Centre
-                    for (int u = PS_MAX(x - size, 0); u <= PS_MIN(x + size, numCols - 1); u++) {
-                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_CONVOLVE;
-                    }
-                    // Right side
-                    for (int u = x + size + 1; u <= PS_MIN(x + footprint + size, numCols - 1); u++) {
-                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
-                    }
-                }
-                // Top stripe
-                for (int v = y + size + 1; v <= PS_MIN(y + footprint + size, numRows - 1); v++) {
-                    for (int u = PS_MAX(x - footprint - size, 0);
-                         u <= PS_MIN(x + footprint + size, numCols - 1); u++) {
-                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
-                    }
-                }
-            }
-        }
-    }
-
-    return mask;
-}
-
 bool pmSubtractionCalculateEquation(psArray *stamps, const psImage *reference, const psImage *input,
                                     const psImage *weight, const pmSubtractionKernels *kernels, int footprint)
@@ -504 +478 @@
                     // Generate the convolutions
                     for (int i = 0; i < numKernels; i++) {
-#if 0
+#ifdef CAN_PASS_STATIC_INLINE
                         convolutions->data.F64[i] = convolvePixel(kernels, i, x, y, reference, polyValues,
                                                                   imageWeighting);
@@ -533 +507 @@
             psFree(polyValues);
 
-            // Fill in lower diagonal of symmetric matrix, while checking for bad values
-            bool bad = false;           // Are there bad values?
+            // Fill in lower diagonal of symmetric matrix
             for (int i = 0; i < numKernels; i++) {
                 for (int j = 0; j < i; j++) {
-                    if (!isfinite(stampMatrix->data.F64[j][i])) {
-                        bad = true;
-                    }
                     stampMatrix->data.F64[i][j] = stampMatrix->data.F64[j][i];
                 }
-                if (!isfinite(stampMatrix->data.F64[i][bgIndex]) ||
-                    !isfinite(stampMatrix->data.F64[i][i]) ||
-                    !isfinite(stampVector->data.F64[i])) {
-                    bad = true;
-                }
                 stampMatrix->data.F64[bgIndex][i] = stampMatrix->data.F64[i][bgIndex];
             }
-            if (!isfinite(stampVector->data.F64[bgIndex])) {
-                bad = true;
-            }
-
-            if (bad) {
-                stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-                psTrace("psModules.imcombine", 3, "Rejecting stamp %d because of bad equation\n", i);
-            } else {
-                stamp->status = PM_SUBTRACTION_STAMP_USED;
-            }
+
+            stamp->status = PM_SUBTRACTION_STAMP_USED;
 
             if (psTraceGetLevel("psModules.imcombine.equation") >= 10) {
@@ -653 +610 @@
 
 int pmSubtractionRejectStamps(psArray *stamps, const psImage *refImage, psImage *inImage,
-                              psImage *subMask, const psVector *solution, int footprint, float sigmaRej,
-                              const pmSubtractionKernels *kernels)
+                              psImage *mask, psMaskType badStampMaskVal, const psVector *solution,
+                              int footprint, float sigmaRej, const pmSubtractionKernels *kernels)
 {
     PS_ASSERT_ARRAY_NON_NULL(stamps, -1);
@@ -662 +619 @@
     PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, -1);
     PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, -1);
-    PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1);
-    PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1);
-    PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, subMask, -1);
+    PS_ASSERT_IMAGE_NON_EMPTY(mask, -1);
+    PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, -1);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, mask, -1);
     PS_ASSERT_VECTOR_NON_NULL(solution, -1);
     PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1);
@@ -707 +664 @@
                 for (int x = footprint; x <= footprint; x++) {
                     for (int k = 0; k < numKernels; k++) {
-#if 0
+#ifdef CAN_PASS_STATIC_INLINE
                         convolvedStamp->data.F32[y + footprint][x + footprint] +=
                             convolvePixel(kernels, k, xStamp + x, yStamp + y, refImage, polyValues,
@@ -753 +710 @@
         if (stamp->status == PM_SUBTRACTION_STAMP_USED && deviations->data.F32[i] > limit) {
             // Mask out the stamp in the image so you it's not found again
-            psTrace("psModules.imcombine", 3, "Clipping stamp %d\n", i);
+            psTrace("psModules.imcombine", 3, "Rejecting stamp %d\n", i);
             numRejected++;
             for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
                 for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
-                    subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
+                    mask->data.PS_TYPE_MASK_DATA[y][x] |= badStampMaskVal;
                 }
             }
@@ -794 +751 @@
 
     // The appropriate kernel
-#if 0
+#ifdef CAN_PASS_STATIC_INLINE
     psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, imageWeighting);
 #else
@@ -810 +767 @@
 
 bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask,
-                           const psImage *inImage, const psImage *inWeight, const psImage *subMask,
-                           psMaskType blank, const psVector *solution, const pmSubtractionKernels *kernels)
+                           const psImage *inImage, const psImage *inWeight, const psImage *inMask,
+                           psMaskType maskVal, psMaskType blank,
+                           const psVector *solution, const pmSubtractionKernels *kernels)
 {
     PS_ASSERT_IMAGE_NON_NULL(inImage, false);
     PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, false);
-    if (subMask) {
-        PS_ASSERT_IMAGE_NON_NULL(subMask, false);
-        PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, inImage, false);
+    if (inMask) {
+        PS_ASSERT_IMAGE_NON_NULL(inMask, false);
+        PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, inImage, false);
     }
     if (inWeight) {
@@ -835 +793 @@
         PS_ASSERT_IMAGES_SIZE_EQUAL(*outMask, inImage, false);
         PS_ASSERT_IMAGE_TYPE(*outMask, PS_TYPE_MASK, false);
-        PS_ASSERT_IMAGE_NON_NULL(subMask, false);
+        PS_ASSERT_IMAGE_NON_NULL(inMask, false);
     }
     if (outWeight && *outWeight) {
@@ -865 +823 @@
     }
     psImage *convMask = NULL;           // Convolved mask image
-    if (outMask && subMask) {
+    if (outMask && inMask) {
         if (!*outMask) {
             *outMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
@@ -896 +854 @@
                                            2.0 * (float)(i + size - numCols/2.0) / (float)numCols,
                                            2.0 * (float)(j + size - numRows/2.0) / (float)numRows);
-#if 0
+#ifdef CAN_PASS_STATIC_INLINE
             kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting);
             if (inWeight) {
@@ -914 +872 @@
 
                     // Check and propagate the kernel footprint, if required
-                    if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] &
-                                    (PM_SUBTRACTION_MASK_INPUT | PM_SUBTRACTION_MASK_CONVOLVE))) {
-                        convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
-                        convImage->data.F32[y][x] = NAN;
-                        if (inWeight) {
-                            convWeight->data.F32[y][x] = NAN;
+                    if (inMask) {
+                        for (int v = -size; v <= size; v++) {
+                            for (int u = -size; u <= size; u++) {
+                                convMask->data.PS_TYPE_MASK_DATA[y][x] |=
+                                    inMask->data.PS_TYPE_MASK_DATA[y + v][x + u];
+                            }
                         }
-                        continue;
+                        if (convMask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) {
+                            convImage->data.F32[y][x] = NAN;
+                            if (inWeight) {
+                                convWeight->data.F32[y][x] = NAN;
+                            }
+                            continue;
+                        }
                     }