Changeset 30289

Timestamp:

Jan 17, 2011, 5:18:09 PM (15 years ago)

Author:

eugene

Message:

apply window only to the moments portions of the calculation; apply weights only to the chisq portion of the calculation; do not renormalize the chisq portion by the flux (keeps SINGLE and DUAL on same scale); move stamp convolution out of CalculateEquation; explicit tests if the convolutions are needed; remove the old unused matrix solution functions; calculate the coefficient errors; add function to calculate chisq and moments statistics, and to generate a subtraction quality statistic; pull model visualization out of CalculateDeviations; adjust the size of the normalization window (if too small, we likely miss some of the flux); exclude excessively faint sources from model fit

Location:

branches/eam_branches/ipp-20101205/psModules/src/imcombine

Files:

• pmSubtractionEquation.c (modified) (47 diffs)
• pmSubtractionEquation.h (modified) (3 diffs)
• pmSubtractionStamps.c (modified) (7 diffs)
• pmSubtractionVisual.c (modified) (11 diffs)
• pmSubtractionVisual.h (modified) (1 diff)

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c

@@ -20 +20 @@
 //# define USE_WEIGHT                      // Include weight (1/variance) in equation?
 # define USE_WINDOW                      // window to avoid neighbor contamination
-// XXX if we want to have a weight and window, we'll need to pass through to the chisq function
+
+/* I believe we want to apply the WEIGHT to the chisq portions of the calculation (but not the WINDOW),
+ * and the WINDOW to the moments portiosn of the calculations (but not the WEIGHT)
+ *
+ */
 
 # define PENALTY false
@@ -106 +110 @@
                         cc *= weight->kernel[y][x];
                     }
-                    if (window) {
+                    // XXX NOTE: do NOT apply the window to the chisq portions of the calculation
+                    if (false && window) {
                         cc *= window->kernel[y][x];
                     }
@@ -137 +142 @@
                     rc *= wtVal;
                 }
-                if (window) {
+                // XXX NOTE: do NOT apply the window to the chisq portions of the calculation
+                if (false && window) {
                     float winVal = window->kernel[y][x];
                     ic *= winVal;
@@ -172 +178 @@
 }
 
+# define RENORM_BY_FLUX 0
 
 // Calculate the least-squares matrix and vector for dual convolution
@@ -267 +274 @@
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
+
+                    // XXX NOTE: clipping low S/N pixels does not seem to work very well
+                    if (false && weight) {
+                        float i1 = image1->kernel[y][x];
+                        float i2 = image2->kernel[y][x];
+                        float sn = (i1 + i2) / sqrt (weight->kernel[y][x]);
+                        if (sn < 0.5) continue;
+                    }
+
                     double aa = iConv1->kernel[y][x] * jConv1->kernel[y][x] * PS_SQR(normValue);
                     double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x] * normValue;
-                    if (weight) {
-                        float wtVal = weight->kernel[y][x];
-                        aa *= wtVal;
-                        bb *= wtVal;
-                        ab *= wtVal;
-                    }
-                    if (window) {
-                        float wtVal = window->kernel[y][x];
-                        aa *= wtVal;
-                        bb *= wtVal;
-                        ab *= wtVal;
-                    }
-                    sumAA += aa;
-                    sumBB += bb;
-                    sumAB += ab;
+
+                    float wtVal = (weight) ? weight->kernel[y][x] : 1.0;
+                    sumAA += wtVal*aa;
+                    sumBB += wtVal*bb;
+                    sumAB += wtVal*ab;
 
                     if (MOMENTS) {
-                        MxxAA += x*x*aa;
-                        MyyAA += y*y*aa;
-                        MxxBB += x*x*bb;
-                        MyyBB += y*y*bb;
+                        float winVal = (window) ? window->kernel[y][x] : 1.0;
+                        MxxAA += winVal*x*x*aa;
+                        MyyAA += winVal*y*y*aa;
+                        MxxBB += winVal*x*x*bb;
+                        MyyBB += winVal*y*y*bb;
                     }
                 }
             }
 
+            // XXX does normSquare1,2 mess up the relative scaling?
+            // XXX no: normSquare1,2 is the sum of the flux^2 for the source
             if (MOMENTS) {
                 MxxAA /= stamp->normSquare1 * PS_SQR(normValue);
@@ -304 +313 @@
             // XXX this makes the Chisq portion independent of the normalization and star flux
             // but may be mis-scaling between stars of different fluxes
+# if (RENORM_BY_FLUX)       
             sumAA /= PS_SQR(stamp->normI2);
             sumAB /= PS_SQR(stamp->normI2);
             sumBB /= PS_SQR(stamp->normI2);
+# endif
 
             // fprintf (stderr, "i,j : %d %d : M(xx,yy)(AA,BB) : %f %f %f %f\n", i, j, MxxAA, MyyAA, MxxBB, MyyBB);
@@ -327 +338 @@
                         matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE;
                         matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE;
-
                         matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE;
                         matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE;
@@ -333 +343 @@
                         matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE;
                         matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE;
-
                         matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE;
                         matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE;
@@ -353 +362 @@
                         ab *= weight->kernel[y][x];
                     }
-                    if (window) {
+                    // XXX NOTE: do NOT apply the window to the chisq portions of the calculation
+                    if (false && window) {
                         ab *= window->kernel[y][x];
                     }
@@ -362 +372 @@
             // XXX this makes the Chisq portion independent of the normalization and star flux
             // but may be mis-scaling between stars of different fluxes
+# if (RENORM_BY_FLUX)
             sumAB /= PS_SQR(stamp->normI2);
+# endif
 
             // Spatial variation of kernel coeffs
@@ -390 +402 @@
                 float i2 = image2->kernel[y][x];
 
+                // XXX NOTE: clipping low S/N pixels does not seem to work very well
+                if (false && weight) {
+                    float sn = (i1 + i2) / sqrt (weight->kernel[y][x]);
+                    if (sn < 0.5) continue;
+                }
+
                 double ai2 = a * i2 * normValue;
                 double bi2 = b * i2;
@@ -395 +413 @@
                 double bi1 = b * i1 * normValue;
 
-                if (weight) {
-                    float wtVal = weight->kernel[y][x];
-                    ai2 *= wtVal;
-                    bi2 *= wtVal;
-                    ai1 *= wtVal;
-                    bi1 *= wtVal;
-                }
-                if (window) {
-                    float wtVal = window->kernel[y][x];
-                    ai2 *= wtVal;
-                    bi2 *= wtVal;
-                    ai1 *= wtVal;
-                    bi1 *= wtVal;
-                }
-                sumAI2 += ai2;
-                sumBI2 += bi2;
-                sumAI1 += ai1;
-                sumBI1 += bi1;
+                float wtVal = (weight) ? weight->kernel[y][x] : 1.0;
+                sumAI2 += wtVal*ai2;
+                sumBI2 += wtVal*bi2;
+                sumAI1 += wtVal*ai1;
+                sumBI1 += wtVal*bi1;
 
                 if (MOMENTS) {
-                    MxxAI1 += x*x*ai1;
-                    MyyAI1 += y*y*ai1;
-                    MxxBI2 += x*x*bi2;
-                    MyyBI2 += y*y*bi2;
+                    float winVal = (window) ? window->kernel[y][x] : 1.0;
+                    MxxAI1 += winVal*x*x*ai1;
+                    MyyAI1 += winVal*y*y*ai1;
+                    MxxBI2 += winVal*x*x*bi2;
+                    MyyBI2 += winVal*y*y*bi2;
                 }
             }
@@ -434 +440 @@
         // XXX this makes the Chisq portion independent of the normalization and star flux
         // but may be mis-scaling between stars of different fluxes
+# if (RENORM_BY_FLUX)
         sumAI1 /= PS_SQR(stamp->normI2);
         sumBI1 /= PS_SQR(stamp->normI2);
         sumAI2 /= PS_SQR(stamp->normI2);
         sumBI2 /= PS_SQR(stamp->normI2);
+# endif
 
         // Spatial variation
@@ -799 +807 @@
     pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
     int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
-    pmSubtractionEquationCalculationMode mode  = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model
-
-    return pmSubtractionCalculateEquationStamp(stamps, kernels, index, mode);
-}
-
-bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
-                                         int index, const pmSubtractionEquationCalculationMode mode)
+
+    return pmSubtractionCalculateEquationStamp(stamps, kernels, index);
+}
+
+bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels, int index)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -832 +838 @@
     psAssert(stamp->status == PM_SUBTRACTION_STAMP_CALCULATE, "We only operate on stamps with this state.");
 
-    // Generate convolutions: these are generated once and saved
-    if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
-        psError(psErrorCodeLast(), false, "Unable to convolve stamp %d.", index);
-        return NULL;
-    }
-
-#ifdef TESTING
-    for (int j = 0; j < numKernels; j++) {
-        if (stamp->convolutions1) {
-            psString convName = NULL;
-            psStringAppend(&convName, "conv1_%03d_%03d.fits", index, j);
-            psFits *fits = psFitsOpen(convName, "w");
-            psFree(convName);
-            psKernel *conv = stamp->convolutions1->data[j];
-            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
-            psFitsClose(fits);
-        }
-
-        if (stamp->convolutions2) {
-            psString convName = NULL;
-            psStringAppend(&convName, "conv2_%03d_%03d.fits", index, j);
-            psFits *fits = psFitsOpen(convName, "w");
-            psFree(convName);
-            psKernel *conv = stamp->convolutions2->data[j];
-            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
-            psFitsClose(fits);
-        }
-    }
-#endif
-
-    // XXX visualize the set of convolved stamps
-
-    psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder,
-                                                    stamp->xNorm, stamp->yNorm); // Polynomial terms
-
-    bool new = stamp->vector ? false : true; // Is this a new run?
+    psImage *polyValues = p_pmSubtractionPolynomial(NULL, spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms
+
+    // Is this a new run? Have we allocated the correct sized vector/matrix?
+    bool new = stamp->vector ? false : true;
+    if (!new && (stamp->vector->n != numParams)) {
+        psFree (stamp->vector);
+        psFree (stamp->matrix);
+        new = true;
+    }
+
     if (new) {
         stamp->matrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
@@ -940 +920 @@
 }
 
-bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels,
-                                    const pmSubtractionEquationCalculationMode mode)
+bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -952 +931 @@
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+
         if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
             continue;
@@ -968 +948 @@
             psArrayAdd(job->args, 1, (pmSubtractionKernels*)kernels); // Casting away const to put on array
             PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32);
             if (!psThreadJobAddPending(job)) {
                 return false;
             }
         } else {
-            pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode);
+            pmSubtractionCalculateEquationStamp(stamps, kernels, i);
         }
     }
@@ -986 +965 @@
     pmSubtractionVisualPlotLeastSquares(stamps);
 
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec",
-             psTimerClear("pmSubtractionCalculateEquation"));
-
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate equation: %f sec", psTimerClear("pmSubtractionCalculateEquation"));
 
     return true;
@@ -997 +974 @@
 bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header);
 
-psImage *p_pmSubSolve_wUt (psVector *w, psImage *U);
-psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt);
-
-bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask);
-
-bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B);
-bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB);
-bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB);
-
-bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x);
-bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y);
-bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
-
-psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w);
-
-double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps);
-
-bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels,
-                                const pmSubtractionStampList *stamps,
-                                const pmSubtractionEquationCalculationMode mode)
+bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps)
 {
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+
+    psTimerStart("pmSubtractionSolveEquation");
 
     // Check inputs
@@ -1041 +1001 @@
         }
 
+        if (stamp->vector->n != numParams) {
+            fprintf (stderr, "mismatch length\n");
+        }
         PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false);
         PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false);
@@ -1088 +1051 @@
 #endif
 
+        // XXX TEST : print the matrix & vector
+        if (0) {
+            for (int iy = 0; iy < sumMatrix->numRows; iy++) {
+                for (int ix = 0; ix < sumMatrix->numCols; ix++) {
+                    fprintf (stderr, "%e  ", sumMatrix->data.F64[iy][ix]);
+                }
+                fprintf (stderr, " : %e\n", sumVector->data.F64[iy]);
+            }
+        }
+
         psImage *invMatrix = NULL;
         psVector *solution = NULL;                       // Solution to equation!
@@ -1096 +1069 @@
         // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
         // SINGLE solution
-        if (1) {
-            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
-            invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
-        } else {
-            psVector *PERM = NULL;
-            psImage *LU = psMatrixLUDecomposition(NULL, &PERM, sumMatrix);
-            solution = psMatrixLUSolution(solution, LU, sumVector, PERM);
-            psFree (LU);
-            psFree (PERM);
-        }
-
 # if (1)
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+        invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
+# endif
+# if (0)
+        psMatrixLUSolve(sumMatrixLU, sumVector);
+        solution = psMemIncrRefCounter(sumVector);
+        invMatrix = psMemIncrRefCounter(sumMatrix);
+# endif
+# if (0)
+        psMatrixGJSolve(sumMatrix, sumVector);
+        invMatrix = psMemIncrRefCounter(sumMatrix);
+        solution = psMemIncrRefCounter(sumVector);
+# endif
+
+# if (0)
         for (int i = 0; i < solution->n; i++) {
-            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(invMatrix->data.F64[i][i]));
+            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Single solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(fabs(invMatrix->data.F64[i][i])));
         }
 # endif
 
-        if (!kernels->solution1) {
-            kernels->solution1 = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
-            psVectorInit(kernels->solution1, 0.0);
-            kernels->solution1err = psVectorAlloc(sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
-            psVectorInit(kernels->solution1err, 0.0);
-        }
+        // ensure we have a solution vector of the right size
+        kernels->solution1    = psVectorRecycle(kernels->solution1,    sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+        kernels->solution1err = psVectorRecycle(kernels->solution1err, sumVector->n + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+        psVectorInit(kernels->solution1, 0.0);
+        psVectorInit(kernels->solution1err, 0.0);
 
         int numKernels = kernels->num;
         int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
         int numPoly = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of polynomial terms
+
         for (int i = 0; i < numKernels * numPoly; i++) {
             kernels->solution1->data.F64[i] = solution->data.F64[i];
@@ -1180 +1157 @@
         }
 
+        // XXX TEST : print the matrix & vector
+        if (0) {
+            for (int iy = 0; iy < sumMatrix->numRows; iy++) {
+                for (int ix = 0; ix < sumMatrix->numCols; ix++) {
+                    fprintf (stderr, "%e  ", sumMatrix->data.F64[iy][ix]);
+                }
+                fprintf (stderr, " : %e\n", sumVector->data.F64[iy]);
+            }
+        }
+
         psImage *invMatrix = NULL;
         psVector *solution = NULL;                       // Solution to equation!
@@ -1187 +1174 @@
         // DUAL solution
 # if (1)
-        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-7);
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
         invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
 # endif
 # if (0)
         psMatrixLUSolve(sumMatrix, sumVector);
-        solution = sumVector;
-        invMatrix = sumMatrix;
+        solution = psMemIncrRefCounter(sumVector);
+        invMatrix = psMemIncrRefCounter(sumMatrix);
 # endif
 
-#if (1)
+#if (0)
         for (int i = 0; i < solution->n; i++) {
             fprintf(stderr, "Dual solution %d: %lf +/- %lf\n", i, solution->data.F64[i], sqrt(invMatrix->data.F64[i][i]));
@@ -1202 +1189 @@
 #endif
 
-        if (!kernels->solution1) {
-            kernels->solution1 = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64);
-            psVectorInit (kernels->solution1, 0.0);
-            kernels->solution1err = psVectorAlloc(numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
-            psVectorInit(kernels->solution1err, 0.0);
-        }
-        if (!kernels->solution2) {
-            kernels->solution2 = psVectorAlloc(numSolution2, PS_TYPE_F64);
-            psVectorInit (kernels->solution2, 0.0);
-            kernels->solution2err = psVectorAlloc(numSolution2, PS_TYPE_F64);
-            psVectorInit(kernels->solution2err, 0.0);
-        }
+        // XXX TEST: manually set the coeffs to a desired solution
+        // solution->data.F64[0] = +1.826;
+        // solution->data.F64[1] = -0.115;
+        // solution->data.F64[2] =  0.0;
+        // solution->data.F64[3] =  0.0;
+
+        // ensure we have solution vectors of the right size
+        kernels->solution1    = psVectorRecycle(kernels->solution1,    numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+        kernels->solution1err = psVectorRecycle(kernels->solution1err, numSolution1 + 2, PS_TYPE_F64); // 1 for norm, 1 for bg
+        kernels->solution2    = psVectorRecycle(kernels->solution2,    numSolution2,     PS_TYPE_F64); // 1 for norm, 1 for bg
+        kernels->solution2err = psVectorRecycle(kernels->solution2err, numSolution2,     PS_TYPE_F64); // 1 for norm, 1 for bg
+
+        psVectorInit(kernels->solution1, 0.0);
+        psVectorInit(kernels->solution1err, 0.0);
+        psVectorInit(kernels->solution2, 0.0);
+        psVectorInit(kernels->solution2err, 0.0);
 
         // for DUAL convolution analysis, we apply the normalization to I1 as follows:
@@ -1259 +1250 @@
         }
      }
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Solve equation: %f sec", psTimerClear("pmSubtractionSolveEquation"));
 
     // pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const
@@ -1311 +1304 @@
 
 // given the convolved image(s) and the residual image, calculate the second moment(s) and the chisq
-bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window) {
+bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window) {
+
+# ifndef USE_WEIGHT
+    psAssert(weight == NULL, "impossible!");
+# endif
+# ifndef USE_WINDOW
+    psAssert(window == NULL, "impossible!");
+# endif
 
     int npix = 0;
@@ -1319 +1319 @@
     for (int y = residual->yMin; y <= residual->yMax; y++) {
         for (int x = residual->xMin; x <= residual->xMax; x++) {
-            chisq  += PS_SQR(residual->kernel[y][x]);
+            float value = PS_SQR(residual->kernel[y][x]);
+            if (weight) {
+                value *= weight->kernel[y][x];
+            }
+            // XXX NOTE: do NOT apply the window to the chisq portions of the calculation
+            if (false && window) {
+                value *= window->kernel[y][x];
+            }
+            chisq += value;
             npix ++;
         }
@@ -1345 +1353 @@
                 value2  = PS_SQR(value1);
 
-                if (weight) {
+                if (window) {
+                    value1 *= window->kernel[y][x];
+                    value2 *= window->kernel[y][x];
+                }
+
+                fluxC1 += value1;
+                flux2  += value2;
+                fluxX  += x * value2;
+                fluxY  += y * value2;
+                // fluxX2 += PS_SQR(x) * value2;
+                // fluxY2 += PS_SQR(y) * value2;
+                fluxX2 += PS_SQR(x) * value1;
+                fluxY2 += PS_SQR(y) * value1;
+            }
+        }
+        // float Mx = fluxX / flux2;
+        // float My = fluxY / flux2;
+        // float Mxx = fluxX2 / flux2;
+        // float Myy = fluxY2 / flux2;
+        float Mxx = fluxX2 / fluxC1;
+        float Myy = fluxY2 / fluxC1;
+
+        // fprintf (stderr, "conv1, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix);
+        moment += Mxx + Myy;
+    }
+
+    // get the moments from convolved1
+    if (convolved2) {
+        for (int y = residual->yMin; y <= residual->yMax; y++) {
+            for (int x = residual->xMin; x <= residual->xMax; x++) {
+                value1  = convolved2->kernel[y][x];
+                value2  = PS_SQR(value1);
+
+                // XXX NOTE: do NOT apply the weight to the moments calculation
+                if (false && weight) {
                     value2 *= weight->kernel[y][x];
                 }
@@ -1353 +1395 @@
                 }
 
-                fluxC1 += value1;
+                fluxC2 += value1;
                 flux2  += value2;
                 fluxX  += x * value2;
                 fluxY  += y * value2;
-                fluxX2 += PS_SQR(x) * value2;
-                fluxY2 += PS_SQR(y) * value2;
+                // fluxX2 += PS_SQR(x) * value2;
+                // fluxY2 += PS_SQR(y) * value2;
+                fluxX2 += PS_SQR(x) * value1;
+                fluxY2 += PS_SQR(y) * value1;
             }
         }
         // float Mx = fluxX / flux2;
         // float My = fluxY / flux2;
-        float Mxx = fluxX2 / flux2;
-        float Myy = fluxY2 / flux2;
-
-        // fprintf (stderr, "conv1, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix);
-        moment += Mxx + Myy;
-    }
-
-    // get the moments from convolved1
-    if (convolved2) {
-        for (int y = residual->yMin; y <= residual->yMax; y++) {
-            for (int x = residual->xMin; x <= residual->xMax; x++) {
-                value1  = convolved2->kernel[y][x];
-                value2  = PS_SQR(value1);
-
-                if (weight) {
-                    value2 *= weight->kernel[y][x];
-                }
-                if (window) {
-                    value1 *= window->kernel[y][x];
-                    value2 *= window->kernel[y][x];
-                }
-
-                fluxC2 += value1;
-                flux2  += value2;
-                fluxX  += x * value2;
-                fluxY  += y * value2;
-                fluxX2 += PS_SQR(x) * value2;
-                fluxY2 += PS_SQR(y) * value2;
-            }
-        }
-        // float Mx = fluxX / flux2;
-        // float My = fluxY / flux2;
-        float Mxx = fluxX2 / flux2;
-        float Myy = fluxY2 / flux2;
+        // float Mxx = fluxX2 / flux2;
+        // float Myy = fluxY2 / flux2;
+        float Mxx = fluxX2 / fluxC2;
+        float Myy = fluxY2 / fluxC2;
 
         // fprintf (stderr, "conv2, flux2: %f, Mx: %f, My: %f, Mxx: %f, Myy: %f, chisq: %f, npix: %d\n", flux2, Mx, My, Mxx, Myy, chisq, npix);
@@ -1410 +1424 @@
     psVectorAppend(momentVector, moment);
     psVectorAppend(fluxesVector, flux);
+    psVectorAppend(stampMask, 0);
     return true;
 }
 
-// typedef struct {
-//     float moment;
-//     float chisq;
-// } diffStat;
-
-bool pmSubtractionCalculateChisqAndMoments(pmSubtractionStampList *stamps,
+bool pmSubtractionCalculateChisqAndMoments(pmSubtractionQuality **bestMatch, 
+                                           pmSubtractionStampList *stamps,
                                            pmSubtractionKernels *kernels)
 {
@@ -1425 +1436 @@
     PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
 
+    psTimerStart("pmSubtractionCalculateChisqAndMoments");
+
+    // XXX need to save these somewhere
     psVector *fluxes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     psVector *chisq = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     psVector *moments = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *stampMask = psVectorAllocEmpty(stamps->num, PS_TYPE_VECTOR_MASK);
 
     int footprint = stamps->footprint; // Half-size of stamps
@@ -1434 +1449 @@
     psImage *polyValues = NULL;         // Polynomial values
 
+    // storage for the image (convolved2 is not used in SINGLE mode)
     psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
-
-    // storage for the convolved source image (only convolved1 is used in SINGLE mode)
     psKernel *convolved1 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
     psKernel *convolved2 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
 
+    int nGood = 0;
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
         if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            // mark this stamp as unused (note that we have to append NANs to the other vectors to keep the lengths in sync)
+            psVectorAppend(moments, NAN);
+            psVectorAppend(fluxes, NAN);
+            psVectorAppend(chisq, NAN);
+            psVectorAppend(stampMask, 0x01);
             continue;
         }
+        nGood ++;
 
         // Calculate coefficients of the kernel basis functions
@@ -1454 +1475 @@
         psImageInit(residual->image, 0.0);
 
-    psKernel *weight = NULL;
-    psKernel *window = NULL;
-
+        psKernel *weight = NULL;
+        psKernel *window = NULL;
+    
 #ifdef USE_WEIGHT
     weight = stamp->weight;
@@ -1511 +1532 @@
             }
             // XXX if we want to have a weight and window, we'll need to pass through to here
-            pmSubtractionChisqStats(fluxes, chisq, moments, convolved1, NULL, residual, weight, window);
+            pmSubtractionChisqStats(fluxes, chisq, moments, stampMask, convolved1, NULL, residual, weight, window);
 
         } else {
@@ -1555 +1576 @@
             } 
 
-            pmSubtractionChisqStats(fluxes, chisq, moments, convolved1, convolved2, residual, weight, window);
-        }
-    }
-
-    pmSubtractionVisualPlotChisqAndMoments(fluxes, chisq, moments);
+            pmSubtractionChisqStats(fluxes, chisq, moments, stampMask, convolved1, convolved2, residual, weight, window);
+        }
+    }
 
     // find the mean chisq and mean moment
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
-    psVectorStats (stats, chisq, NULL, NULL, 0);
+    psVectorStats (stats, chisq, NULL, stampMask, 0xff);
     float chisqValue = stats->sampleMean;
 
     psStatsInit(stats);
-    psVectorStats (stats, moments, NULL, NULL, 0);
+    psVectorStats (stats, moments, NULL, stampMask, 0xff);
     float momentValue = stats->sampleMean;
 
-    fprintf (stderr, "chisq: %f, moment: %f\n", chisqValue, momentValue);
+    // XXX this is probably wrong : I want to score against higher moments **compared with the raw momements**
+    // score : chisq * dMoments * modeFactor
+    float modeFactor = kernels->mode == PM_SUBTRACTION_MODE_DUAL ? 2.0 : 1.0;
+    float orderFactor = 0.5*(kernels->spatialOrder + 1) * (kernels->spatialOrder + 2);
+    float score = chisqValue * momentValue * modeFactor * orderFactor;
+
+    fprintf (stderr, "chisq: %f, moment: %f, score: %f\n", chisqValue, momentValue, score);
+
+    // save this result if it is the first or the best (skip if bestMatch is NULL)
+    if (bestMatch) {
+        pmSubtractionQuality *match = *bestMatch;
+        bool keep = false;
+        if (match == NULL) {
+            *bestMatch = match = pmSubtractionQualityAlloc();
+            keep = true;
+        } else {
+            if (score < match->score) {
+                psFree(match->fluxes);
+                psFree(match->chisq);
+                psFree(match->moments);
+                psFree(match->stampMask);
+                keep = true;
+            }
+        }
+        if (keep) {
+            fprintf (stderr, "keeping order: %d, mode: %d, score: %f\n", kernels->spatialOrder, kernels->mode, score);
+            match->score        = score;
+            match->spatialOrder = kernels->spatialOrder;
+            match->mode         = kernels->mode;
+            match->nGood        = nGood;
+            match->fluxes       = psMemIncrRefCounter(fluxes);
+            match->chisq        = psMemIncrRefCounter(chisq);
+            match->moments      = psMemIncrRefCounter(moments);
+            match->stampMask    = psMemIncrRefCounter(stampMask);
+        }           
+    }
+
+    pmSubtractionVisualPlotChisqAndMoments(fluxes, chisq, moments);
 
     psFree(stats);
@@ -1576 +1632 @@
     psFree(fluxes);
     psFree(moments);
+    psFree(stampMask);
 
     psFree(residual);
@@ -1582 +1639 @@
     psFree(polyValues);
 
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate Chisq and Moments: %f sec", psTimerClear("pmSubtractionCalculateChisqAndMoments"));
+
     return true;
 }
 
-psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps,
-                                           pmSubtractionKernels *kernels)
+// XXX for now, let's not use this, and let's instead just use values from pmSubtractionCalculateChisqAndMoments
+psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
     PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
+
+    psTimerStart("pmSubtractionCalculateDeviations");
 
     psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps
@@ -1600 +1661 @@
     psImage *polyValues = NULL;         // Polynomial values
     psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
-
-    // set up holding images for the visualization
-    pmSubtractionVisualShowFitInit (stamps);
 
     psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
@@ -1703 +1761 @@
             }
 
-            // XXX visualize the target, source, convolution and residual
-            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
-
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
@@ -1740 +1795 @@
             }
 
-            // XXX visualize the target, source, convolution and residual
-            pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
-
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
@@ -1757 +1809 @@
         }
 
+        double flux = 0.0;
         double deviation = 0.0;         // Sum of differences
         for (int y = - footprint; y <= footprint; y++) {
@@ -1762 +1815 @@
                 double dev = PS_SQR(residual->kernel[y][x]) * weight->kernel[y][x];
                 deviation += dev;
-#ifdef TESTING
-                residual->kernel[y][x] = dev;
-#endif
+                flux += stamp->image1->kernel[y][x] + stamp->image2->kernel[y][x];
             }
         }
         deviations->data.F32[i] = devNorm * deviation;
-        psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
-                i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]);
+        psTrace("psModules.imcombine", 5, "Deviation and Flux for stamp %d (%d,%d): %f %f\n",
+                i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i], flux);
         psStringAppend(&log, "Stamp %d (%d,%d): %f\n",
                        i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5), deviations->data.F32[i]);
         if (!isfinite(deviations->data.F32[i])) {
             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-            psTrace("psModules.imcombine", 5,
-                    "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
-                    i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
+            psTrace("psModules.imcombine", 5, "Rejecting stamp %d (%d,%d) because of non-finite deviation\n", i, (int)(stamp->x - 0.5), (int)(stamp->y - 0.5));
             continue;
         }
-
-#ifdef TESTING
-        {
-            psString filename = NULL;
-            psStringAppend(&filename, "resid_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        if (stamp->image1) {
-            psString filename = NULL;
-            psStringAppend(&filename, "stamp_image1_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        if (stamp->image2) {
-            psString filename = NULL;
-            psStringAppend(&filename, "stamp_image2_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, stamp->image2->image, 0, NULL);
-            psFitsClose(fits);
-        }
-        if (stamp->weight) {
-            psString filename = NULL;
-            psStringAppend(&filename, "stamp_weight_%03d.fits", i);
-            psFits *fits = psFitsOpen(filename, "w");
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL);
-            psFitsClose(fits);
-        }
-#endif
-
     }
 
@@ -1829 +1842 @@
         psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
 
-        pmSubtractionVisualShowFit(norm);
-        pmSubtractionVisualPlotFit(kernels);
-
         psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
         psVectorStats (stats, fResSigma, NULL, NULL, 0);
@@ -1862 +1872 @@
     psFree(polyValues);
 
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Calculate Deviations: %f sec", psTimerClear("pmSubtractionCalculateDeviations"));
+
     return deviations;
-}
-
-// we are supplied U, not Ut; w represents a diagonal matrix (also, we apply 1/w instead of w)
-psImage *p_pmSubSolve_wUt (psVector *w, psImage *U) {
-
-    psAssert (w->n == U->numCols, "w and U dimensions do not match");
-
-    // wUt has dimensions transposed relative to Ut.
-    psImage *wUt = psImageAlloc (U->numRows, U->numCols, PS_TYPE_F64);
-    psImageInit (wUt, 0.0);
-
-    for (int i = 0; i < wUt->numCols; i++) {
-        for (int j = 0; j < wUt->numRows; j++) {
-            if (!isfinite(w->data.F64[j])) continue;
-            if (w->data.F64[j] == 0.0) continue;
-            wUt->data.F64[j][i] = U->data.F64[i][j] / w->data.F64[j];
-        }
-    }
-    return wUt;
-}
-
-// XXX this is just standard matrix multiplication: use psMatrixMultiply?
-psImage *p_pmSubSolve_VwUt (psImage *V, psImage *wUt) {
-
-    psAssert (V->numCols == wUt->numRows, "matrix dimensions do not match");
-
-    psImage *Ainv = psImageAlloc (wUt->numCols, V->numRows, PS_TYPE_F64);
-
-    for (int i = 0; i < Ainv->numCols; i++) {
-        for (int j = 0; j < Ainv->numRows; j++) {
-            double sum = 0.0;
-            for (int k = 0; k < V->numCols; k++) {
-                sum += V->data.F64[j][k] * wUt->data.F64[k][i];
-            }
-            Ainv->data.F64[j][i] = sum;
-        }
-    }
-    return Ainv;
-}
-
-// we are supplied U, not Ut
-bool p_pmSubSolve_UtB (psVector **UtB, psImage *U, psVector *B) {
-
-    psAssert (U->numRows == B->n, "U and B dimensions do not match");
-
-    UtB[0] = psVectorRecycle (UtB[0], U->numCols, PS_TYPE_F64);
-
-    for (int i = 0; i < U->numCols; i++) {
-        double sum = 0.0;
-        for (int j = 0; j < U->numRows; j++) {
-            sum += B->data.F64[j] * U->data.F64[j][i];
-        }
-        UtB[0]->data.F64[i] = sum;
-    }
-    return true;
-}
-
-// w is diagonal
-bool p_pmSubSolve_wUtB (psVector **wUtB, psVector *w, psVector *UtB) {
-
-    psAssert (w->n == UtB->n, "w and UtB dimensions do not match");
-
-    // wUt has dimensions transposed relative to Ut.
-    wUtB[0] = psVectorRecycle (wUtB[0], w->n, PS_TYPE_F64);
-    psVectorInit (wUtB[0], 0.0);
-
-    for (int i = 0; i < w->n; i++) {
-        if (!isfinite(w->data.F64[i])) continue;
-        if (w->data.F64[i] == 0.0) continue;
-        wUtB[0]->data.F64[i] = UtB->data.F64[i] / w->data.F64[i];
-    }
-    return true;
-}
-
-// this is basically matrix * vector
-bool p_pmSubSolve_VwUtB (psVector **VwUtB, psImage *V, psVector *wUtB) {
-
-    psAssert (V->numCols == wUtB->n, "V and wUtB dimensions do not match");
-
-    VwUtB[0] = psVectorRecycle (*VwUtB, V->numRows, PS_TYPE_F64);
-
-    for (int j = 0; j < V->numRows; j++) {
-        double sum = 0.0;
-        for (int i = 0; i < V->numCols; i++) {
-            sum += V->data.F64[j][i] * wUtB->data.F64[i];
-        }
-        VwUtB[0]->data.F64[j] = sum;
-    }
-    return true;
-}
-
-// this is basically matrix * vector
-bool p_pmSubSolve_Ax (psVector **B, psImage *A, psVector *x) {
-
-    psAssert (A->numCols == x->n, "A and x dimensions do not match");
-
-    B[0] = psVectorRecycle (*B, A->numRows, PS_TYPE_F64);
-
-    for (int j = 0; j < A->numRows; j++) {
-        double sum = 0.0;
-        for (int i = 0; i < A->numCols; i++) {
-            sum += A->data.F64[j][i] * x->data.F64[i];
-        }
-        B[0]->data.F64[j] = sum;
-    }
-    return true;
-}
-
-// this is basically Vector * vector
-bool p_pmSubSolve_VdV (double *value, psVector *x, psVector *y) {
-
-    psAssert (x->n == y->n, "x and y dimensions do not match");
-
-    double sum = 0.0;
-    for (int i = 0; i < x->n; i++) {
-        sum += x->data.F64[i] * y->data.F64[i];
-    }
-    *value = sum;
-    return true;
-}
-
-bool p_pmSubSolve_y2 (double *y2, pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
-
-    int footprint = stamps->footprint; // Half-size of stamps
-
-    double sum = 0.0;
-    for (int i = 0; i < stamps->num; i++) {
-
-        pmSubtractionStamp *stamp = stamps->stamps->data[i];
-        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
-
-        psKernel *weight = NULL;
-        psKernel *window = NULL;
-        psKernel *input = NULL;
-
-#ifdef USE_WEIGHT
-        weight = stamp->weight;
-#endif
-#ifdef USE_WINDOW
-        window = stamps->window;
-#endif
-
-        switch (kernels->mode) {
-            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
-          case PM_SUBTRACTION_MODE_1:
-            input = stamp->image2;
-            break;
-          case PM_SUBTRACTION_MODE_2:
-            input = stamp->image1;
-            break;
-          default:
-            psAbort ("programming error");
-        }
-
-        for (int y = - footprint; y <= footprint; y++) {
-            for (int x = - footprint; x <= footprint; x++) {
-                double in = input->kernel[y][x];
-                double value = in*in;
-                if (weight) {
-                    float wtVal = weight->kernel[y][x];
-                    value *= wtVal;
-                }
-                if (window) {
-                    float  winVal = window->kernel[y][x];
-                    value *= winVal;
-                }
-                sum += value;
-            }
-        }
-    }
-    *y2 = sum;
-    return true;
-}
-
-double p_pmSubSolve_ChiSquare (pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) {
-
-    int footprint = stamps->footprint; // Half-size of stamps
-    int numKernels = kernels->num;      // Number of kernels
-
-    double sum = 0.0;
-
-    psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
-    psImageInit(residual->image, 0.0);
-
-    psImage *polyValues = NULL;         // Polynomial values
-
-    for (int i = 0; i < stamps->num; i++) {
-
-        pmSubtractionStamp *stamp = stamps->stamps->data[i];
-        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
-
-        psKernel *weight = NULL;
-        psKernel *window = NULL;
-        psKernel *target = NULL;
-        psKernel *source = NULL;
-
-        psArray *convolutions = NULL;
-
-#ifdef USE_WEIGHT
-        weight = stamp->weight;
-#endif
-#ifdef USE_WINDOW
-        window = stamps->window;
-#endif
-
-        switch (kernels->mode) {
-            // MODE_1 : convolve image 1 to match image 2 (and vice versa)
-          case PM_SUBTRACTION_MODE_1:
-            target = stamp->image2;
-            source = stamp->image1;
-            convolutions = stamp->convolutions1;
-            break;
-          case PM_SUBTRACTION_MODE_2:
-            target = stamp->image1;
-            source = stamp->image2;
-            convolutions = stamp->convolutions2;
-            break;
-          default:
-            psAbort ("programming error");
-        }
-
-        // Calculate coefficients of the kernel basis functions
-        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
-        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
-        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
-
-        psImageInit(residual->image, 0.0);
-        for (int j = 0; j < numKernels; j++) {
-            psKernel *convolution = convolutions->data[j]; // Convolution
-            double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient;
-                }
-            }
-        }
-
-        for (int y = - footprint; y <= footprint; y++) {
-            for (int x = - footprint; x <= footprint; x++) {
-                double resid = target->kernel[y][x] - background - source->kernel[y][x] * norm + residual->kernel[y][x];
-                double value = PS_SQR(resid);
-                if (weight) {
-                    float wtVal = weight->kernel[y][x];
-                    value *= wtVal;
-                }
-                if (window) {
-                    float  winVal = window->kernel[y][x];
-                    value *= winVal;
-                }
-                sum += value;
-            }
-        }
-    }
-    psFree (polyValues);
-    psFree (residual);
-
-    return sum;
-}
-
-bool p_pmSubSolve_SetWeights (psVector *wApply, psVector *w, psVector *wMask) {
-
-    for (int i = 0; i < w->n; i++) {
-        wApply->data.F64[i] = wMask->data.U8[i] ? 0.0 : w->data.F64[i];
-    }
-    return true;
-}
-
-// we are supplied V and w; w represents a diagonal matrix (also, we apply 1/w instead of w)
-psImage *p_pmSubSolve_Xvar (psImage *V, psVector *w) {
-
-    psAssert (w->n == V->numCols, "w and U dimensions do not match");
-
-    psImage *Vn = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64);
-    psImageInit (Vn, 0.0);
-
-    // generate Vn = V * w^{-1}
-    for (int j = 0; j < Vn->numRows; j++) {
-        for (int i = 0; i < Vn->numCols; i++) {
-            if (!isfinite(w->data.F64[i])) continue;
-            if (w->data.F64[i] == 0.0) continue;
-            Vn->data.F64[j][i] = V->data.F64[j][i] / w->data.F64[i];
-        }
-    }
-
-    psImage *Xvar = psImageAlloc (V->numCols, V->numRows, PS_TYPE_F64);
-    psImageInit (Xvar, 0.0);
-
-    // generate Xvar = Vn * Vn^T
-    for (int j = 0; j < Vn->numRows; j++) {
-        for (int i = 0; i < Vn->numCols; i++) {
-            double sum = 0.0;
-            for (int k = 0; k < Vn->numCols; k++) {
-                sum += Vn->data.F64[k][i]*Vn->data.F64[k][j];
-            }
-            Xvar->data.F64[j][i] = sum;
-        }
-    }
-    return Xvar;
 }
 
@@ -2166 +1880 @@
 // of the elements of an image A(x,y) = A->data.F64[y][x] = A_x,y, a matrix
 // multiplication is: A_k,j * B_i,k = C_i,j
-
 
 bool psFitsWriteImageSimple (char *filename, psImage *image, psMetadata *header) {

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.h

@@ -4 +4 @@
 #include "pmSubtractionStamps.h"
 #include "pmSubtractionKernels.h"
+#include "pmSubtraction.h"
 
-typedef enum {
-    PM_SUBTRACTION_EQUATION_NONE    = 0x00,
-    PM_SUBTRACTION_EQUATION_NORM    = 0x01,
-    PM_SUBTRACTION_EQUATION_BG      = 0x02,
-    PM_SUBTRACTION_EQUATION_KERNELS = 0x04,
-    PM_SUBTRACTION_EQUATION_ALL     = 0x07, // value should be NORM | BG | KERNELS
-} pmSubtractionEquationCalculationMode;
+// typedef enum {
+//     PM_SUBTRACTION_EQUATION_NONE    = 0x00,
+//     PM_SUBTRACTION_EQUATION_NORM    = 0x01,
+//     PM_SUBTRACTION_EQUATION_BG      = 0x02,
+//     PM_SUBTRACTION_EQUATION_KERNELS = 0x04,
+//     PM_SUBTRACTION_EQUATION_ALL     = 0x07, // value should be NORM | BG | KERNELS
+// } pmSubtractionEquationCalculationMode;
 
 /// Execute a thread job to calculate the least-squares equation for a stamp
@@ -20 +21 @@
 bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, ///< Stamps
                                          pmSubtractionKernels *kernels, ///< Kernel parameters
-                                         int index, ///< Index of stamp
-                                         const pmSubtractionEquationCalculationMode mode
+                                         int index ///< Index of stamp
     );
 
 /// Calculate the least-squares equation to match the image quality
 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
-                                    pmSubtractionKernels *kernels, ///< Kernel parameters
-                                    const pmSubtractionEquationCalculationMode mode
+                                    pmSubtractionKernels *kernels ///< Kernel parameters
     );
 
 /// Solve the least-squares equation to match the image quality
 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters
-                                const pmSubtractionStampList *stamps, ///< Stamps
-                                const pmSubtractionEquationCalculationMode mode
+                                const pmSubtractionStampList *stamps ///< Stamps
     );
 
@@ -92 +90 @@
 bool pmSubtractionCalculateMomentsKernel(double *Mxx, double *Myy, psKernel *image, int footprint, int window);
 
-bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window);
+bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window);
 
-bool pmSubtractionCalculateChisqAndMoments(pmSubtractionStampList *stamps,
-                                           pmSubtractionKernels *kernels);
-
+bool pmSubtractionCalculateChisqAndMoments(pmSubtractionQuality **bestMatch, pmSubtractionStampList *stamps, pmSubtractionKernels *kernels);
 #endif

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c

@@ -122 +122 @@
             if ((image1 && image1->data.F32[y][x] < thresh1) ||
                 (image2 && image2->data.F32[y][x] < thresh2)) {
+                // fprintf (stderr, "%f,%f : thresh\n", xRaw, yRaw);
                 continue;
             }
@@ -429 +430 @@
         stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing,
                                              normFrac, sysErr, skyErr);
+    }
+
+    // XXX TEST : dump all stars in the stamps here
+    if (0) {
+        FILE *f = fopen ("stamp.dat", "w");
+        for (int i = 0; i < stamps->num; i++) {
+            psVector *xList = stamps->x->data[i];
+            psVector *yList = stamps->y->data[i]; // Coordinate lists
+            psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes
+
+            for (int j = 0; j < xList->n; j++) {
+                fprintf (f, "%d %d  %f %f  %f\n", i, j, xList->data.F32[j], yList->data.F32[j], fluxList->data.F32[j]);
+            }
+        }
+        fclose (f);
     }
 
@@ -636 +652 @@
     }
 
+    int nTotal = 0;
+
     // Sort the list by flux, with the brightest last
     for (int i = 0; i < numStamps; i++) {
@@ -662 +680 @@
         stamps->y->data[i] = ySorted;
         stamps->flux->data[i] = fluxSorted;
-    }
-
+        nTotal += num;
+    }
+    // fprintf (stderr, "nTotal %d\n", nTotal);
+    
     return stamps;
 }
@@ -809 +829 @@
     float R2 = Sr2 / Sf2;
 
-    stamps->normWindow1 = 2.0*R1;
-    stamps->normWindow2 = 2.0*R2;
+    stamps->normWindow1 = 2.75*R1;
+    stamps->normWindow2 = 2.75*R2;
     psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2);
 
     // Generate a weighting window based on the kron radii
     { 
-        float radius = 1.5 * PS_MAX(R1, R2);
+        float radius = 2.0 * PS_MAX(R1, R2);
 
         psImageInit(stamps->window->image, 0.0);
@@ -826 +846 @@
         }
     }
+
+    // complain if normWindow1 or normWindow2 are larger than size
+    psAssert(stamps->normWindow1 < size, "norm Window 1 too large");
+    psAssert(stamps->normWindow2 < size, "norm Window 2 too large");
 
 # else
@@ -1142 +1166 @@
             continue;
         }
+       
+        // XXX TEST
+        if (source->errMag > 0.1) continue;
+
         if (source->modelPSF) {
             x->data.F32[numOut] = source->modelPSF->params->data.F32[PM_PAR_XPOS];

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.c

@@ -20 +20 @@
 #include "pmSubtractionEquation.h"
 #include "pmSubtractionKernels.h"
+#include "pmSubtractionVisual.h"
 
 #include "pmVisual.h"
@@ -423 +424 @@
     }
 
-    // XXX clear the overlay(s) (red at least!)
+    // clear the overlay (red at least!)
+    KiiEraseOverlay (kapa2, "red");
 
     // get the kernel sizes
@@ -456 +458 @@
     int nKernels = 0;
 
+    // paste in the kernel images, scaled by sum2
     if (maxStamp->convolutions1) {
         // output image is a grid of NXsub by NYsub sub-images
@@ -478 +481 @@
             int yPix = ySub * (2*footprint + 1 + 3) + footprint;
             
-            double sum = 0.0;
             double sum2 = 0.0;
             for (int y = -footprint; y <= footprint; y++) {
                 for (int x = -footprint; x <= footprint; x++) {
-                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
-                    sum += kernel->kernel[y][x];
                     sum2 += PS_SQR(kernel->kernel[y][x]);
                 }
             }
-            // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
+            float scale = sqrt(sum2) / PS_SQR(2*footprint + 1);
+            for (int y = -footprint; y <= footprint; y++) {
+                for (int x = -footprint; x <= footprint; x++) {
+                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x] / scale;
+                }
+            }
         }               
         pmVisualScaleImage(kapa2, output, "Image", 0, true);
@@ -520 +525 @@
             int yPix = ySub * (2*footprint + 1 + 3) + footprint;
             
-            double sum = 0.0;
             double sum2 = 0.0;
             for (int y = -footprint; y <= footprint; y++) {
                 for (int x = -footprint; x <= footprint; x++) {
-                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
-                    sum += kernel->kernel[y][x];
                     sum2 += PS_SQR(kernel->kernel[y][x]);
                 }
             }
-            // fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
+            float scale = sqrt(sum2) / PS_SQR(2*footprint + 1);
+            for (int y = -footprint; y <= footprint; y++) {
+                for (int x = -footprint; x <= footprint; x++) {
+                    output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x] / scale;
+                }
+            }
         }               
         pmVisualScaleImage(kapa2, output, "Image", 1, true);
@@ -587 +594 @@
     KiiLoadOverlay (kapa2, overlay, Noverlay, "red");
     FREE (overlay);
-    return true;
-}
-
-static int footprint = 0;
-static int NX = 0;
-static int NY = 0;
-static psImage *sourceImage      = NULL;
-static psImage *targetImage      = NULL;
-static psImage *residualImage    = NULL;
-static psImage *fresidualImage   = NULL;
-static psImage *differenceImage  = NULL;
-static psImage *convolutionImage = NULL;
-
-bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
-
-    if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
-
-    // generate 4 storage images large enough to hold the N stamps:
-
-    footprint = stamps->footprint;
-
-    float NXf = sqrt(stamps->num);
-    NX = (int) NXf == NXf ? NXf : NXf + 1.0;
-    
-    float NYf = stamps->num / NX;
-    NY = (int) NYf == NY ? NYf : NYf + 1.0;
-
-    int NXpix = (2*footprint + 1) * NX;
-    NXpix += (NX > 1) ? 3 * NX : 0;
-
-    int NYpix = (2*footprint + 1) * NY;
-    NYpix += (NY > 1) ? 3 * NY : 0;
-
-    sourceImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    targetImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    residualImage    = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    fresidualImage   = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
-    
-    psImageInit (sourceImage,      0.0);
-    psImageInit (targetImage,      0.0);
-    psImageInit (residualImage,    0.0);
-    psImageInit (fresidualImage,   0.0);
-    psImageInit (differenceImage,  0.0);
-    psImageInit (convolutionImage, 0.0);
-
-    return true;
-}
-
-bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
-
-    if (!pmVisualTestLevel("ppsub.stamp", 1)) return true;
-
-    double sum;
-
-    int NXoff = index % NX;
-    int NYoff = index / NX;
-
-    int NXpix = NXoff * (2*footprint + 1 + 3) + footprint;
-    int NYpix = NYoff * (2*footprint + 1 + 3) + footprint;
-
-    // insert the (target) kernel into the (target) image:
-    sum = 0.0;
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
-            sum += targetImage->data.F32[y + NYpix][x + NXpix];
-        }
-    }
-    targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
-
-    // insert the (source) kernel into the (source) image:
-    sum = 0.0;
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
-            sum += sourceImage->data.F32[y + NYpix][x + NXpix];
-        }
-    }
-    sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
-
-    // insert the (convolution) kernel into the (convolution) image:
-    sum = 0.0;
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            convolutionImage->data.F32[y + NYpix][x + NXpix] = convolution->kernel[y][x];
-            sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
-        }
-    }
-    convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
-    
-    // insert the (difference) kernel into the (difference) image:
-    sum = 0.0;
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
-            sum += differenceImage->data.F32[y + NYpix][x + NXpix];
-        }
-    }
-    differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
-
-    // insert the (residual) kernel into the (residual) image:
-    sum = 0.0;
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x];
-            sum += residualImage->data.F32[y + NYpix][x + NXpix];
-        }
-    }
-    residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
-
-    // insert the (fresidual) kernel into the (fresidual) image:
-    for (int y = -footprint; y <= footprint; y++) {
-        for (int x = -footprint; x <= footprint; x++) {
-            fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
-        }
-    }
     return true;
 }
@@ -730 +619 @@
 }
 
-bool pmSubtractionVisualShowFit(double norm) {
-
-    // for testing, dump the residual image and exit
-    if (0) {
-        psMetadata *header = psMetadataAlloc();
-        psMetadataAddF32 (header, PS_LIST_TAIL, "NORM", 0, "Normalization", norm);
-        psFits *fits = psFitsOpen("resid.fits", "w");
-        psFitsWriteImage(fits, header, residualImage, 0, NULL);
-        psFitsClose(fits);
-        // exit (0);
-    }
+static int footprint = 0;
+static int NX = 0;
+static int NY = 0;
+static psImage *sourceImage      = NULL;
+static psImage *targetImage      = NULL;
+static psImage *residualImage    = NULL;
+static psImage *fresidualImage   = NULL;
+static psImage *differenceImage  = NULL;
+static psImage *convolutionImage = NULL;
+
+bool pmSubtractionVisualShowFit(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels) {
 
     if (!pmVisualTestLevel("ppsub.fit", 1)) return true;
@@ -746 +635 @@
     if (!pmVisualInitWindow(&kapa1, "ppSub:Images")) return false;
     if (!pmVisualInitWindow(&kapa2, "ppSub:Misc")) return false;
+
+    // set up holding images for the visualization
+    pmSubtractionVisualShowFitInit (stamps);
+
+    int numKernels = kernels->num;      // Number of kernels
+
+    psImage *polyValues = NULL;         // Polynomial values
+    psKernel *residual = psKernelAlloc(-stamps->footprint, stamps->footprint, -stamps->footprint, stamps->footprint); // Residual image
+
+    double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) { continue; }
+
+        // Calculate coefficients of the kernel basis functions
+        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
+        double background = p_pmSubtractionSolutionBackground(kernels, polyValues); // Difference in background
+
+        psImageInit(residual->image, 0.0);
+
+        if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
+            psKernel *target;           // Target postage stamp
+            psKernel *source;           // Source postage stamp
+            psArray *convolutions;      // Convolution postage stamps for each kernel basis function
+            switch (kernels->mode) {
+              case PM_SUBTRACTION_MODE_1:
+                target = stamp->image2;
+                source = stamp->image1;
+                convolutions = stamp->convolutions1;
+                break;
+              case PM_SUBTRACTION_MODE_2:
+                target = stamp->image1;
+                source = stamp->image2;
+                convolutions = stamp->convolutions2;
+                break;
+              default:
+                psAbort("Unsupported subtraction mode: %x", kernels->mode);
+            }
+
+            for (int j = 0; j < numKernels; j++) {
+                psKernel *convolution = convolutions->data[j]; // Convolution
+                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient
+                for (int y = - footprint; y <= footprint; y++) {
+                    for (int x = - footprint; x <= footprint; x++) {
+                        residual->kernel[y][x] += convolution->kernel[y][x] * coefficient;
+                    }
+                }
+            }
+            // visualize the target, source, convolution and residual
+            pmSubtractionVisualShowFitAddStamp (target, source, residual, background, norm, i);
+        } else {
+            // Dual convolution
+            psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image
+            psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
+            psKernel *image1 = stamp->image1; // The first image
+            psKernel *image2 = stamp->image2; // The second image
+
+            for (int j = 0; j < numKernels; j++) {
+                psKernel *conv1 = convolutions1->data[j]; // Convolution of first image
+                psKernel *conv2 = convolutions2->data[j]; // Convolution of second image
+                double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1
+                double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2
+
+                for (int y = - footprint; y <= footprint; y++) {
+                    for (int x = - footprint; x <= footprint; x++) {
+                        residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 + conv1->kernel[y][x] * coeff1;
+                    }
+                }
+            }
+            // visualize the target, source, convolution and residual
+            pmSubtractionVisualShowFitAddStamp (image2, image1, residual, background, norm, i);
+        }
+    }
+    pmSubtractionVisualShowFitImage(norm);
+
+    return true;
+}
+
+// generate 4 storage images large enough to hold the stamps:
+bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps) {
+
+    footprint = stamps->footprint;
+
+    float NXf = sqrt(stamps->num);
+    NX = (int) NXf == NXf ? NXf : NXf + 1.0;
+    
+    float NYf = stamps->num / NX;
+    NY = (int) NYf == NY ? NYf : NYf + 1.0;
+
+    int NXpix = (2*footprint + 1) * NX;
+    NXpix += (NX > 1) ? 3 * NX : 0;
+
+    int NYpix = (2*footprint + 1) * NY;
+    NYpix += (NY > 1) ? 3 * NY : 0;
+
+    sourceImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    targetImage      = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    residualImage    = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    fresidualImage   = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    differenceImage  = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    convolutionImage = psImageAlloc (NXpix, NYpix, PS_TYPE_F32);
+    
+    psImageInit (sourceImage,      0.0);
+    psImageInit (targetImage,      0.0);
+    psImageInit (residualImage,    0.0);
+    psImageInit (fresidualImage,   0.0);
+    psImageInit (differenceImage,  0.0);
+    psImageInit (convolutionImage, 0.0);
+
+    return true;
+}
+
+bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index) {
+
+    double sum;
+
+    int NXoff = index % NX;
+    int NYoff = index / NX;
+
+    int NXpix = NXoff * (2*footprint + 1 + 3) + footprint;
+    int NYpix = NYoff * (2*footprint + 1 + 3) + footprint;
+
+    // insert the (target) kernel into the (target) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            targetImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x];
+            sum += targetImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    targetImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (source) kernel into the (source) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            sourceImage->data.F32[y + NYpix][x + NXpix] = source->kernel[y][x];
+            sum += sourceImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    sourceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (convolution) kernel into the (convolution) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            convolutionImage->data.F32[y + NYpix][x + NXpix] = convolution->kernel[y][x];
+            sum += convolutionImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    convolutionImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+    
+    // insert the (difference) kernel into the (difference) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            differenceImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm;
+            sum += differenceImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    differenceImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (residual) kernel into the (residual) image:
+    sum = 0.0;
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            residualImage->data.F32[y + NYpix][x + NXpix] = target->kernel[y][x] - background - source->kernel[y][x] * norm - convolution->kernel[y][x];
+            sum += residualImage->data.F32[y + NYpix][x + NXpix];
+        }
+    }
+    residualImage->data.F32[footprint + 1 + NYpix][NXpix] = sum;
+
+    // insert the (fresidual) kernel into the (fresidual) image:
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            fresidualImage->data.F32[y + NYpix][x + NXpix] = residualImage->data.F32[y + NYpix][x + NXpix] / sqrt(PS_MAX(target->kernel[y][x], 100.0));
+        }
+    }
+    return true;
+}
+
+bool pmSubtractionVisualShowFitImage(double norm) {
 
     KiiEraseOverlay (kapa1, "red");
@@ -792 +864 @@
     psVector *x = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
     psVector *y = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
+    psVector *dy = psVectorAllocEmpty (kernels->num, PS_TYPE_F32);
 
     graphdata.xmin = -1.0;
@@ -804 +877 @@
         x->data.F32[i] = i;
         y->data.F32[i] = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, false);
+        dy->data.F32[i] = kernels->solution1err->data.F64[i];
         graphdata.ymin = PS_MIN(graphdata.ymin, y->data.F32[i]);
         graphdata.ymax = PS_MAX(graphdata.ymax, y->data.F32[i]);
     }
-    x->n = y->n = kernels->num;
+    x->n = y->n = dy->n = kernels->num;
 
     float range;
@@ -828 +902 @@
     graphdata.size = 0.5;
     graphdata.style = 2;
+    graphdata.etype |= 0x01;
 
     KapaPrepPlot   (kapa3, x->n, &graphdata);
     KapaPlotVector (kapa3, x->n, x->data.F32, "x");
     KapaPlotVector (kapa3, x->n, y->data.F32, "y");
+    KapaPlotVector (kapa3, x->n, dy->data.F32, "dym");
+    KapaPlotVector (kapa3, x->n, dy->data.F32, "dyp");
 
     psFree (x);
     psFree (y);
+    psFree (dy);
     psFree (polyValues);
 

branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionVisual.h

@@ -8 +8 @@
 bool pmSubtractionVisualPlotLeastSquaresResid (const pmSubtractionStampList *stamps, psImage *matrixIn, int nUsed);
 bool pmSubtractionVisualShowSubtraction(psImage *image, psImage *ref, psImage *sub);
+
+bool pmSubtractionVisualShowFit(pmSubtractionStampList *stamps, pmSubtractionKernels *kernels);
 bool pmSubtractionVisualShowFitInit(pmSubtractionStampList *stamps);
 bool pmSubtractionVisualShowFitAddStamp(psKernel *target, psKernel *source, psKernel *convolution, double background, double norm, int index);
-bool pmSubtractionVisualShowFit(double norm);
+bool pmSubtractionVisualShowFitImage(double norm);
+
 bool pmSubtractionVisualPlotFit(const pmSubtractionKernels *kernels);
 bool pmSubtractionVisualShowKernels(pmSubtractionKernels *kernels);