Changeset 26739 for branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c

Timestamp:

Jan 29, 2010, 6:02:38 PM (16 years ago)

Author:

Paul Price

Message:

Fairly extensive change to make subtractions more stable on images with limited lit area in the presence of spatial variation. Took the polynomials which were scaled over the entire image and scale them only over the lit area. Removed some fprintf statements. Fixed the problem causing the background term to always come out 1.0.

File:

• branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c (modified) (5 diffs)

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

@@ -424 +424 @@
                 } else if (*source & subConvPoor) {
                     *target |= maskPoor;
+                } else {
+                    *target &= ~maskBad & ~maskPoor;
                 }
             }
@@ -574 +576 @@
 }
 
+void p_pmSubtractionPolynomialNormCoords(float *xOut, float *yOut, float xIn, float yIn,
+                                         int xMin, int xMax, int yMin, int yMax)
+{
+    float xNormSize = xMax - xMin, yNormSize = yMax - yMin; // Size to use for normalisation
+    *xOut = 2.0 * (float)(xIn - xMin - xNormSize/2.0) / xNormSize;
+    *yOut = 2.0 * (float)(yIn - yMin - yNormSize/2.0) / yNormSize;
+    return;
+}
+
 psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, const pmSubtractionKernels *kernels,
-                                             int numCols, int numRows, int x, int y)
+                                             int x, int y)
 {
     assert(kernels);
-    assert(numCols > 0 && numRows > 0);
-
-    // Size to use when calculating normalised coordinates (different from actual size when convolving
-    // subimage)
-    int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
-    int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
-
-    // Normalised coordinates
-    float yNorm = 2.0 * (float)(y - yNormSize/2.0) / (float)yNormSize;
-    float xNorm = 2.0 * (float)(x - xNormSize/2.0) / (float)xNormSize;
-
+
+    float xNorm, yNorm;                 // Normalised coordinates
+    p_pmSubtractionPolynomialNormCoords(&xNorm, &yNorm, x, y,
+                                        kernels->xMin, kernels->xMax, kernels->yMin, kernels->yMax);
     return p_pmSubtractionPolynomial(output, kernels->spatialOrder, xNorm, yNorm);
 }
@@ -1072 +1076 @@
     // Only generate polynomial values every kernel footprint, since we have already assumed
     // (with the stamps) that it does not vary rapidly on this scale.
-    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
-                                                              xMin + x0 + size + 1,
-                                                              yMin + y0 + size + 1);
+    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, xMin + x0 + size + 1,
+                                                              yMin + y0 + size + 1);        // Polynomial
     float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term
 
@@ -1200 +1203 @@
     bool threaded = pmSubtractionThreaded(); // Running threaded?
 
-    // Outputs
-    if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        if (!out1->image) {
-            out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-        }
-        if (ro1->variance) {
-            if (!out1->variance) {
-                out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            }
-            psImageInit(out1->variance, 0.0);
-        }
-    }
-    if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        if (!out2->image) {
-            out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-        }
-        if (ro2->variance) {
-            if (!out2->variance) {
-                out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            }
-            psImageInit(out2->variance, 0.0);
-        }
-    }
     psImage *convMask = NULL;           // Convolved mask image (common to inputs 1 and 2)
     if (subMask) {
         if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-            if (!out1->mask) {
-                out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-            }
-            psImageInit(out1->mask, 0);
             convMask = out1->mask;
         }
         if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-            if (!out2->mask) {
-                out2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-            }
-            psImageInit(out2->mask, 0);
             if (!convMask) {
                 convMask = out2->mask;
@@ -1256 +1228 @@
 
     // Get region for convolution: [xMin:xMax,yMin:yMax]
-    int xMin = size, xMax = numCols - size;
-    int yMin = size, yMax = numRows - size;
+    int xMin = kernels->xMin + size, xMax = kernels->xMax - size;
+    int yMin = kernels->yMin + size, yMax = kernels->yMax - size;
     if (region) {
         xMin = PS_MAX(region->x0, xMin);