Changeset 26332

Timestamp:

Dec 3, 2009, 12:20:59 PM (16 years ago)

Author:

eugene

Message:

add code to calculate norm separately from kernel coeffs

Location:

branches/eam_branches/20091201/psModules/src/imcombine

Files:

• pmSubtraction.c (modified) (1 diff)
• pmSubtractionEquation.c (modified) (21 diffs)
• pmSubtractionEquation.h (modified) (2 diffs)
• pmSubtractionMatch.c (modified) (4 diffs)
• pmSubtractionVisual.c (modified) (7 diffs)

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

@@ -462 +462 @@
 {
     psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
-#if 1
+#if 0
     // Convolving using separable convolution
     psVector *xKernel = preCalc->data[0]; // Kernel in x

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

@@ -26 +26 @@
 
 // Calculate the least-squares matrix and vector
+// XXX why are we calculating these values on each iteration?  aren't they identical (no change in pixel values...)
 static bool calculateMatrixVector(psImage *matrix, // Least-squares matrix, updated
                                   psVector *vector, // Least-squares vector, updated
@@ -36 +37 @@
                                   const psImage *polyValues, // Spatial polynomial values
                                   int footprint, // (Half-)Size of stamp
-                                  const bool normOnly
+                                  const pmSubtractionEquationCalculationMode mode
                                   )
 {
@@ -54 +55 @@
 
     // Evaluate polynomial-polynomial terms
-    // XXX we can skip this for normOnly
+    // XXX we can skip this if we are not calculating kernel coeffs
     for (int iyOrder = 0, iIndex = 0; iyOrder <= spatialOrder; iyOrder++) {
         for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex++) {
@@ -68 +69 @@
     }
 
-
-    // if we only want to calculate the normalization (and background model?)
-    // then we should assign these values in the matrix and vector to 1 (diagonal) or 0 (off-diagonal)
-    // XXX need to disable normalization calculation if !normOnly
-    if (normOnly) {
-        psImageInit(matrix, 0.0);
-        psVectorInit(vector, 1.0);
-        for (int i = 0; i < matrix->numCols; i++) {
-            matrix->data.F64[i][i] = 1.0;
-        }
-    } else {
-        for (int i = 0; i < numKernels; i++) {
-            psKernel *iConv = convolutions->data[i]; // Convolution for index i
-            for (int j = i; j < numKernels; j++) {
-                psKernel *jConv = convolutions->data[j]; // Convolution for index j
-
-                double sumCC = 0.0;         // Sum of convolution products
-                for (int y = - footprint; y <= footprint; y++) {
-                    for (int x = - footprint; x <= footprint; x++) {
-                        double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
-#ifdef USE_WEIGHT
+    // initialize the matrix and vector for NOP on all coeffs.  we only fill in the coeffs we 
+    // choose to calculate
+    psImageInit(matrix, 0.0);
+    psVectorInit(vector, 1.0);
+    for (int i = 0; i < matrix->numCols; i++) {
+        matrix->data.F64[i][i] = 1.0;
+    }
+
+    // the order of the elements in the matrix and vector is:
+    // [kernel 0, x^0 y^0][kernel 1 x^0 y^0]...[kernel N, x^0 y^0]
+    // [kernel 0, x^1 y^0][kernel 1 x^1 y^0]...[kernel N, x^1 y^0]
+    // [kernel 0, x^n y^m][kernel 1 x^n y^m]...[kernel N, x^n y^m]
+    // normalization
+    // bg 0, bg 1, bg 2 (only 0 is currently used?)
+
+    for (int i = 0; i < numKernels; i++) {
+        psKernel *iConv = convolutions->data[i]; // Convolution for index i
+        for (int j = i; j < numKernels; j++) {
+            psKernel *jConv = convolutions->data[j]; // Convolution for index j
+
+            double sumCC = 0.0;         // Sum of convolution products
+            for (int y = - footprint; y <= footprint; y++) {
+                for (int x = - footprint; x <= footprint; x++) {
+                    double cc = iConv->kernel[y][x] * jConv->kernel[y][x];
+                    if (weight) {
                         cc *= weight->kernel[y][x];
-#endif
-#ifdef USE_WINDOW
-                        if (window) {
-                            cc *= window->kernel[y][x];
-                        }
-#endif
-                        sumCC += cc;
                     }
+                    if (window) {
+                        cc *= window->kernel[y][x];
+                    }
+                    sumCC += cc;
                 }
-
-                // Spatial variation
+            }
+
+            // Spatial variation of kernel coeffs
+            if (mode | PM_SUBTRACTION_EQUATION_KERNELS) {
                 for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
                     for (int jTerm = 0, jIndex = j; jTerm < numPoly; jTerm++, jIndex += numKernels) {
@@ -109 +113 @@
                 }
             }
-
-            double sumRC = 0.0;             // Sum of the reference-convolution products
-            double sumIC = 0.0;             // Sum of the input-convolution products
-            double sumC = 0.0;              // Sum of the convolution
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    float conv = iConv->kernel[y][x];
-                    float in = input->kernel[y][x];
-                    float ref = reference->kernel[y][x];
-                    double ic = in * conv;
-                    double rc = ref * conv;
-                    double c = conv;
-#ifdef USE_WEIGHT
+        }
+
+        double sumRC = 0.0;             // Sum of the reference-convolution products
+        double sumIC = 0.0;             // Sum of the input-convolution products
+        double sumC = 0.0;              // Sum of the convolution
+        for (int y = - footprint; y <= footprint; y++) {
+            for (int x = - footprint; x <= footprint; x++) {
+                float conv = iConv->kernel[y][x];
+                float in = input->kernel[y][x];
+                float ref = reference->kernel[y][x];
+                double ic = in * conv;
+                double rc = ref * conv;
+                double c = conv;
+                if (weight) {
                     float wtVal = weight->kernel[y][x];
                     ic *= wtVal;
                     rc *= wtVal;
                     c *= wtVal;
-#endif
-#ifdef USE_WINDOW
-                    if (window) {
-                        float winVal = window->kernel[y][x];
-                        ic *= winVal;
-                        rc *= winVal;
-                        c  *= winVal;
-                    }
-#endif
-                    sumIC += ic;
-                    sumRC += rc;
-                    sumC += c;
                 }
+                if (window) {
+                    float winVal = window->kernel[y][x];
+                    ic *= winVal;
+                    rc *= winVal;
+                    c  *= winVal;
+                }
+                sumIC += ic;
+                sumRC += rc;
+                sumC += c;
             }
-            // Spatial variation
-            for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
-                double normTerm = sumRC * poly[iTerm];
-                double bgTerm = sumC * poly[iTerm];
+        }
+        // Spatial variation
+        for (int iTerm = 0, iIndex = i; iTerm < numPoly; iTerm++, iIndex += numKernels) {
+            double normTerm = sumRC * poly[iTerm];
+            double bgTerm = sumC * poly[iTerm];
+            if ((mode | PM_SUBTRACTION_EQUATION_NORM) && (mode | PM_SUBTRACTION_EQUATION_KERNELS)) {
                 matrix->data.F64[iIndex][normIndex] = normTerm;
                 matrix->data.F64[normIndex][iIndex] = normTerm;
+            }
+            if ((mode | PM_SUBTRACTION_EQUATION_BG) && (mode | PM_SUBTRACTION_EQUATION_KERNELS)) {
                 matrix->data.F64[iIndex][bgIndex] = bgTerm;
                 matrix->data.F64[bgIndex][iIndex] = bgTerm;
+            }
+            if (mode | PM_SUBTRACTION_EQUATION_KERNELS) {
                 vector->data.F64[iIndex] = sumIC * poly[iTerm];
             }
         }
     }
-
-    // XXX need to disable normalization calculation if !normOnly
 
     double sumRR = 0.0;                 // Sum of the reference product
@@ -167 +173 @@
             double rr = PS_SQR(ref);
             double one = 1.0;
-#ifdef USE_WEIGHT
-            float wtVal = weight->kernel[y][x];
-            rr *= wtVal;
-            ir *= wtVal;
-            in *= wtVal;
-            ref *= wtVal;
-            one *= wtVal;
-#endif
-#ifdef USE_WINDOW
+            if (weight) {
+                float wtVal = weight->kernel[y][x];
+                rr *= wtVal;
+                ir *= wtVal;
+                in *= wtVal;
+                ref *= wtVal;
+                one *= wtVal;
+            }
             if (window) {
                 float  winVal = window->kernel[y][x];
@@ -184 +189 @@
                 one *= winVal;
             }
-#endif
             sumRR += rr;
             sumIR += ir;
@@ -192 +196 @@
         }
     }
-    matrix->data.F64[normIndex][normIndex] = sumRR;
-    matrix->data.F64[bgIndex][bgIndex] = sum1;
-    matrix->data.F64[normIndex][bgIndex] = matrix->data.F64[bgIndex][normIndex] = sumR;
-    vector->data.F64[normIndex] = sumIR;
-    vector->data.F64[bgIndex] = sumI;
-
+    if (mode | PM_SUBTRACTION_EQUATION_NORM) {
+        matrix->data.F64[normIndex][normIndex] = sumRR;
+        vector->data.F64[normIndex] = sumIR;
+    }
+    if (mode | PM_SUBTRACTION_EQUATION_BG) {
+        matrix->data.F64[bgIndex][bgIndex] = sum1;
+        vector->data.F64[bgIndex] = sumI;
+    }
+    if ((mode | PM_SUBTRACTION_EQUATION_NORM) && (mode | PM_SUBTRACTION_EQUATION_BG)) {
+        matrix->data.F64[normIndex][bgIndex] = sumR;
+        matrix->data.F64[bgIndex][normIndex] = sumR;
+    }
     return true;
 }
@@ -267 +277 @@
                     double bb = iConv2->kernel[y][x] * jConv2->kernel[y][x];
                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
-#ifdef USE_WEIGHT
-                    float wtVal = weight->kernel[y][x];
-                    aa *= wtVal;
-                    bb *= wtVal;
-                    ab *= wtVal;
-#endif
+                    if (weight) {
+                        float wtVal = weight->kernel[y][x];
+                        aa *= wtVal;
+                        bb *= wtVal;
+                        ab *= wtVal;
+                    }
                     sumAA += aa;
                     sumBB += bb;
@@ -299 +309 @@
                 for (int x = - footprint; x <= footprint; x++) {
                     double ab = iConv1->kernel[y][x] * jConv2->kernel[y][x];
-#ifdef USE_WEIGHT
-                    ab *= weight->kernel[y][x];
-#endif
+                    if (weight) {
+                        ab *= weight->kernel[y][x];
+                    }
                     sumAB += ab;
                 }
@@ -336 +346 @@
                 double i1i2 = i1 * i2;
 
-#ifdef USE_WEIGHT
-                float wtVal = weight->kernel[y][x];
-                ai2 *= wtVal;
-                bi2 *= wtVal;
-                ai1 *= wtVal;
-                bi1 *= wtVal;
-                i1i2 *= wtVal;
-                a *= wtVal;
-                b *= wtVal;
-                i2 *= wtVal;
-#endif
-
+                if (weight) {
+                    float wtVal = weight->kernel[y][x];
+                    ai2 *= wtVal;
+                    bi2 *= wtVal;
+                    ai1 *= wtVal;
+                    bi1 *= wtVal;
+                    i1i2 *= wtVal;
+                    a *= wtVal;
+                    b *= wtVal;
+                    i2 *= wtVal;
+                }
                 sumAI2 += ai2;
                 sumBI2 += bi2;
@@ -392 +401 @@
             double i1i2 = i1 * i2;
 
-#ifdef USE_WEIGHT
-            float wtVal = weight->kernel[y][x];
-            i1 *= wtVal;
-            i1i1 *= wtVal;
-            one *= wtVal;
-            i2 *= wtVal;
-            i1i2 *= wtVal;
-#endif
-
+            if (weight) {
+                float wtVal = weight->kernel[y][x];
+                i1 *= wtVal;
+                i1i1 *= wtVal;
+                one *= wtVal;
+                i2 *= wtVal;
+                i1i2 *= wtVal;
+            }
             sumI1 += i1;
             sumI1I1 += i1i1;
@@ -588 +596 @@
     const pmSubtractionKernels *kernels = job->args->data[1]; // Kernels
     int index = PS_SCALAR_VALUE(job->args->data[2], S32); // Stamp index
-
-    return pmSubtractionCalculateEquationStamp(stamps, kernels, index, normOnly);
+    pmSubtractionEquationCalculationMode mode  = PS_SCALAR_VALUE(job->args->data[3], S32); // calculation model
+
+    return pmSubtractionCalculateEquationStamp(stamps, kernels, index, mode);
 }
 
 bool pmSubtractionCalculateEquationStamp(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
-                                         int index, bool normOnly)
+                                         int index, const pmSubtractionEquationCalculationMode mode)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -662 +671 @@
 
     bool status;                    // Status of least-squares matrix/vector calculation
+
+    psKernel *weight = NULL;
+    psKernel *window = NULL;
+
+#ifdef USE_WEIGHT
+    weight = stamp->weight;
+#endif
+#ifdef USE_WINDOW
+    window = stamps->window;
+#endif
+
     switch (kernels->mode) {
       case PM_SUBTRACTION_MODE_1:
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image2, stamp->image1,
-                                       stamp->weight, stamps->window, stamp->convolutions1, kernels,
-                                       polyValues, footprint, normOnly);
+                                       weight, window, stamp->convolutions1, kernels,
+                                       polyValues, footprint, mode);
         break;
       case PM_SUBTRACTION_MODE_2:
         status = calculateMatrixVector(stamp->matrix1, stamp->vector1, stamp->image1, stamp->image2,
-                                       stamp->weight, stamps->window, stamp->convolutions2, kernels,
-                                       polyValues, footprint, normOnly);
+                                       weight, window, stamp->convolutions2, kernels,
+                                       polyValues, footprint, mode);
         break;
       case PM_SUBTRACTION_MODE_DUAL:
-        psAbort ("dual is disabled: need to add window and normOnly");
+        psAbort ("dual is disabled: need to add window and calculation mode");
         if (new) {
             stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64);
@@ -686 +706 @@
 #endif
         status = calculateDualMatrixVector(stamp->matrix1, stamp->vector1, stamp->matrix2, stamp->vector2,
-                                           stamp->matrixX, stamp->image1, stamp->image2, stamp->weight,
+                                           stamp->matrixX, stamp->image1, stamp->image2, weight,
                                            stamp->convolutions1, stamp->convolutions2, kernels, polyValues,
                                            footprint);
@@ -758 +778 @@
 }
 
-bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const bool normOnly)
+bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels, const pmSubtractionEquationCalculationMode mode)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -778 +798 @@
             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, normOnly, PS_TYPE_BOOL);
+            PS_ARRAY_ADD_SCALAR(job->args, mode, PS_TYPE_S32);
             if (!psThreadJobAddPending(job)) {
                 psFree(job);
@@ -785 +805 @@
             psFree(job);
         } else {
-            pmSubtractionCalculateEquationStamp(stamps, kernels, i, normOnly);
+            pmSubtractionCalculateEquationStamp(stamps, kernels, i, mode);
         }
     }
@@ -803 +823 @@
 }
 
-bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps)
+bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, 
+                                const pmSubtractionStampList *stamps, 
+                                const pmSubtractionEquationCalculationMode mode)
 {
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
@@ -941 +963 @@
             return NULL;
         }
-        kernels->solution1 = psMatrixLUSolution(kernels->solution1, luMatrix, sumVector, permutation);
+
+        psVector *solution = psMatrixLUSolution(NULL, luMatrix, sumVector, permutation);
+        psFree(sumVector);
+        psFree(luMatrix);
+        psFree(permutation);
+        if (!solution) {
+            psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
+            return NULL;
+        }
+
+        if (!kernels->solution1) {
+            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
+            psVectorInit (kernels->solution1, 0.0);
+        }
+
+        // only update the solutions that we chose to calculate:
+        if (mode | PM_SUBTRACTION_EQUATION_NORM) {
+            int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            kernels->solution1->data.F64[normIndex] = solution->data.F64[normIndex];
+        }
+        if (mode | PM_SUBTRACTION_EQUATION_BG) {
+            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
+            kernels->solution1->data.F64[bgIndex] = solution->data.F64[bgIndex];
+        }
+        if (mode | PM_SUBTRACTION_EQUATION_KERNELS) {
+            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];
+            }
+        }
+        psFree (solution);
 
 #ifdef TESTING
@@ -952 +1006 @@
 #endif
 
-        psFree(sumVector);
-        psFree(luMatrix);
-        psFree(permutation);
-        if (!kernels->solution1) {
-            psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
-            return NULL;
-        }
     } else {
         // Dual convolution solution

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h

@@ -4 +4 @@
 #include "pmSubtractionStamps.h"
 #include "pmSubtractionKernels.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     = 0x05, // value should be NORM | BG | KERNELS
+} pmSubtractionEquationCalculationMode;
 
 /// Execute a thread job to calculate the least-squares equation for a stamp
@@ -13 +21 @@
                                          const pmSubtractionKernels *kernels, ///< Kernel parameters
                                          int index, ///< Index of stamp
-                                         const bool normOnly
-                                    );
+                                         const pmSubtractionEquationCalculationMode mode
+    );
 
 /// Calculate the least-squares equation to match the image quality
 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
                                     const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                    const bool normOnly
-                                    );
+                                    const pmSubtractionEquationCalculationMode mode
+    );
 
 /// Solve the least-squares equation to match the image quality
 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters
-                                const pmSubtractionStampList *stamps ///< Stamps
+                                const pmSubtractionStampList *stamps, ///< Stamps
+                                const pmSubtractionEquationCalculationMode mode
     );
 

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

@@ -561 +561 @@
 
                 // XXX step 1: calculate normalization
-                psTrace("psModules.imcombine", 3, "Calculating equation...\n");
-                if (!pmSubtractionCalculateEquation(stamps, kernels, true)) {
+                psTrace("psModules.imcombine", 3, "Calculating equation for normalization...\n");
+                // if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_NORM)) {
+                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
                 }
 
-                // XXX step 2: calculate only other parameters
+                psTrace("psModules.imcombine", 3, "Solving equation for normalization...\n");
+                // if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_NORM)) {
+                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+                memCheck("  solve equation");
+
 # if (0)
-                psTrace("psModules.imcombine", 3, "Calculating equation...\n");
-                if (!pmSubtractionCalculateEquation(stamps, kernels)) {
+                // XXX step 2: calculate kernel parameters
+                psTrace("psModules.imcombine", 3, "Calculating equation for kernels...\n");
+                if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
                 }
+                memCheck("  calculate equation");
+
+                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
+                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+                    goto MATCH_ERROR;
+                }
+                memCheck("  solve equation");
 # endif
-
-                memCheck("  calculate equation");
-
-                psTrace("psModules.imcombine", 3, "Solving equation...\n");
-
-                if (!pmSubtractionSolveEquation(kernels, stamps)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto MATCH_ERROR;
-                }
-
-                memCheck("  solve equation");
 
                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
@@ -609 +615 @@
             // if we hit the max number of iterations and we have rejected stamps, re-solve
             if (numRejected > 0) {
-                psTrace("psModules.imcombine", 3, "Solving equation...\n");
-                if (!pmSubtractionSolveEquation(kernels, stamps)) {
+                // calculate kernel parameters
+                psTrace("psModules.imcombine", 3, "Solving equation for kernels...\n");
+                // if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
+                if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
                 }
+                memCheck("  solve equation");
+
                 psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
                 if (!deviations) {
@@ -918 +928 @@
     assert(kernels);
 
-    psTrace("psModules.imcombine", 3, "Calculating %s equation...\n", description);
-    if (!pmSubtractionCalculateEquation(stamps, kernels)) {
+    psTrace("psModules.imcombine", 3, "Calculating %s normalization equation...\n", description);
+    // if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_NORM)) {
+    if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_ALL)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
         return false;
     }
 
-    psTrace("psModules.imcombine", 3, "Solving %s equation...\n", description);
-    if (!pmSubtractionSolveEquation(kernels, stamps)) {
+    psTrace("psModules.imcombine", 3, "Solving %s normalization equation...\n", description);
+    // if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_NORM)) {
+    if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
         return false;
     }
+
+# if (0)
+    psTrace("psModules.imcombine", 3, "Calculating %s kernel coeffs equation...\n", description);
+    if (!pmSubtractionCalculateEquation(stamps, kernels, PM_SUBTRACTION_EQUATION_KERNELS)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+
+    psTrace("psModules.imcombine", 3, "Solving %s kernel coeffs equation...\n", description);
+    if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return false;
+    }
+# endif
 
     psTrace("psModules.imcombine", 3, "Calculate %s deviations...\n", description);
@@ -949 +975 @@
         // Allow re-fit with reduced stamps set
         psTrace("psModules.imcombine", 3, "Resolving %s equation...\n", description);
-        if (!pmSubtractionSolveEquation(kernels, stamps)) {
+        // if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_KERNELS)) {
+        if (!pmSubtractionSolveEquation(kernels, stamps, PM_SUBTRACTION_EQUATION_ALL)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
             return false;
         }
         psTrace("psModules.imcombine", 3, "Recalculate %s deviations...\n", description);
+
         psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
         if (!deviations) {

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

@@ -263 +263 @@
 static psImage *targetImage      = NULL;
 static psImage *residualImage    = NULL;
+static psImage *fresidualImage   = NULL;
 static psImage *differenceImage  = NULL;
 static psImage *convolutionImage = NULL;
@@ -289 +290 @@
     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);
@@ -295 +297 @@
     psImageInit (targetImage,      0.0);
     psImageInit (residualImage,    0.0);
+    psImageInit (fresidualImage,   0.0);
     psImageInit (differenceImage,  0.0);
     psImageInit (convolutionImage, 0.0);
@@ -343 +346 @@
         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];
+        }
+    }
+
+    // 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] / target->kernel[y][x];
         }
     }
@@ -358 +368 @@
     pmVisualScaleImage(kapa1, convolutionImage, "Convolution Stamps", 2, true);
 
+    // pmVisualScaleImage(kapa2, fresidualImage, "Frac Residual Stamps", 2, true);
+    pmVisualScaleImage(kapa2, residualImage, "Residual Stamps", 1, true);
     pmVisualScaleImage(kapa2, differenceImage, "Difference Stamps", 0, true);
-    pmVisualScaleImage(kapa2, residualImage, "Residual Stamps", 1, true);
     pmVisualAskUser(NULL);
 
@@ -367 +378 @@
     psFree(differenceImage);
     psFree(residualImage);
+    psFree(fresidualImage);
 
     targetImage = NULL;
@@ -373 +385 @@
     differenceImage = NULL;
     residualImage = NULL;
+    fresidualImage = NULL;
 
     return true;