Changeset 14318

Timestamp:

Jul 19, 2007, 11:26:12 AM (19 years ago)

Author:

Paul Price

Message:

Fixing MAJOR bug --- the rejection step wasn't doing the convolution
correctly. It was, in effect, adding each kernel component, instead
of using the solution. Fixed to look more like the code within
pmSubtractionConvolve.

File:

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

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-19 03:36:06 $
+ *  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-07-19 21:26:12 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -25 +25 @@
 #include "pmSubtraction.h"
 
-//#define USE_FUNCTIONS_INSTEAD_OF_MACROS          // Can I pass the address of a static inline function?
+#define USE_FUNCTIONS_INSTEAD_OF_MACROS          // Can I pass the address of a static inline function?
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -851 +851 @@
     }
 
-    int numKernels = solution->n - 1;   // Number of kernel basis functions
-
-    // Measure
+    // Measure deviations
     psVector *deviations = psVectorAlloc(stamps->n, PS_TYPE_F32); // Mean deviation for stamps
     double totalSquareDev = 0.0;        // Total square deviation from zero
@@ -859 +857 @@
     {
         int numCols = refImage->numCols, numRows = refImage->numRows; // Image dimensions
-
         psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics
-
-        // Convolved image of stamp
-        psImage *convolvedStamp = psImageAlloc(2 * footprint + 1, 2 * footprint + 1, PS_TYPE_F32);
+        psImage *convolvedStamp = psImageAlloc(2 * footprint + 1, 2 * footprint + 1,
+                                               PS_TYPE_F32); // Convolved image of stamp
+        psKernel *kernelImage = NULL;       // The kernel, with which to convolve the stamps
+        float background = solution->data.F64[solution->n-1]; // The difference in background
+        int size = kernels->size;       // Half-size of the kernel
 
         for (int i = 0; i < stamps->n; i++) {
@@ -871 +870 @@
             }
             int xStamp = stamp->x, yStamp = stamp->y; // Coordinates of stamp
-            psImageInit(convolvedStamp, 0.0);
 
             // Precalulate polynomial values
@@ -878 +876 @@
                                                     2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows);
 
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = footprint; x <= footprint; x++) {
-                    for (int k = 0; k < numKernels; k++) {
 #ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS
-                        convolvedStamp->data.F32[y + footprint][x + footprint] +=
-                            convolvePixel(kernels, k, xStamp + x, yStamp + y, refImage, polyValues,
-                                          imageWeighting);
+            kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting);
 #else
-                        CONVOLVE_PIXEL(convolvedStamp->data.F32[y + footprint][x + footprint], kernels,
-                                       k, xStamp + x, yStamp + y, refImage, polyValues, imageWeighting);
+            SOLVED_KERNEL(kernelImage, solution, kernels, polyValues, imageWeighting);
 #endif
+            for (int y = yStamp - footprint, j = 0; y <= yStamp + footprint; y++, j++) {
+                for (int x = xStamp - footprint, i = 0; x <= xStamp + footprint; x++, i++) {
+                    convolvedStamp->data.F32[j][i] = background;
+                    for (int v = -size; v <= size; v++) {
+                        for (int u = -size; u <= size; u++) {
+                            convolvedStamp->data.F32[j][i] += kernelImage->kernel[v][u] *
+                                refImage->data.F32[y + v][x + u];
+                        }
                     }
                 }
             }
+
             psFree(polyValues);
 
@@ -900 +901 @@
                    inStamp->numRows == convolvedStamp->numRows);
             (void)psBinaryOp(convolvedStamp, inStamp, "-", convolvedStamp);
+            (void)psBinaryOp(convolvedStamp, convolvedStamp, "*", convolvedStamp);
             (void)psBinaryOp(convolvedStamp, convolvedStamp, "/", inStamp);
             psFree(inStamp);
-            (void)psBinaryOp(convolvedStamp, convolvedStamp, "*", convolvedStamp);
             if (!psImageStats(stats, convolvedStamp, NULL, 0)) {
                 psError(PS_ERR_UNKNOWN, false,
@@ -911 +912 @@
                 return -1;
             }
-            deviations->data.F32[i] = stats->sampleMean;
+            deviations->data.F32[i] = sqrt(stats->sampleMean / 2.0);
+            psTrace("psModules.imcombine", 1, "Deviation for stamp %d: %f\n", i, deviations->data.F32[i]);
             totalSquareDev += PS_SQR(deviations->data.F32[i]);
             numStamps++;
         }
 
+        psFree(kernelImage);
         psFree(convolvedStamp);
         psFree(stats);
     }
 
-    float limit = sigmaRej * sqrt(totalSquareDev / (numStamps - 1)); // Limit on maximum deviation
+    float limit = sigmaRej * sqrt(totalSquareDev / numStamps); // Limit on maximum deviation
+    psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit);
 
     int numRejected = 0;