Changeset 15809

Timestamp:

Dec 13, 2007, 11:32:44 AM (18 years ago)

Author:

Paul Price

Message:

Updates trying to get dual-convolution to work. It doesn't yet.

Location:

trunk/psModules/src/imcombine

Files:

• pmSubtraction.c (modified) (7 diffs)
• pmSubtractionEquation.c (modified) (29 diffs)
• pmSubtractionEquation.h (modified) (1 diff)
• pmSubtractionKernels.c (modified) (4 diffs)
• pmSubtractionMatch.c (modified) (11 diffs)
• pmSubtractionStamps.h (modified) (2 diffs)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-12-07 01:57:15 $
+ *  @version $Revision: 1.75 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-12-13 21:32:44 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -165 +165 @@
     }
 
-    if (!wantDual) {
-        // Put in the normalisation component
-        kernel->kernel[0][0] += p_pmSubtractionSolutionNorm(kernels);
-    }
-
+    // Put in the normalisation component
+    kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels));
 
     return kernel;
@@ -279 +276 @@
                                   float background, // Background value to apply
                                   psRegion region, // Region to convolve
-                                  bool useFFT  // Use FFT to convolve?
+                                  bool useFFT,  // Use FFT to convolve?
+                                  bool wantDual // Want the dual convolution?
     )
 {
-    *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, false);
+    *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual);
     if (weight) {
         *kernelWeight = varianceKernel(*kernelWeight, *kernelImage);
@@ -559 +557 @@
                 stamp->x = 0;
                 stamp->y = 0;
+                stamp->xNorm = NAN;
+                stamp->yNorm = NAN;
                 stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
                 // Recalculate convolutions
@@ -568 +568 @@
                 psFree(stamp->weight);
                 stamp->image1 = stamp->image2 = stamp->weight = NULL;
+                psFree(stamp->matrix1);
+                psFree(stamp->matrix2);
+                psFree(stamp->matrixX);
+                stamp->matrix1 = stamp->matrix2 = stamp->matrixX = NULL;
+                psFree(stamp->vector1);
+                psFree(stamp->vector2);
+                stamp->vector1 = stamp->vector2 = NULL;
             } else {
                 numGood++;
@@ -651 +658 @@
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false);
-    if (kernels->mode != PM_SUBTRACTION_MODE_1) {
+    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         PS_ASSERT_PTR_NON_NULL(out2, false);
         PS_ASSERT_PTR_NON_NULL(ro2, false);
@@ -812 +819 @@
 
             convolveRegion(convImage1, convWeight1, &kernelImage, &kernelWeight, image, weight,
-                           subMask, maskSource, kernels, polyValues, background, subRegion, useFFT);
+                           subMask, maskSource, kernels, polyValues, background, subRegion, useFFT, false);
 
             if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
                 convolveRegion(convImage2, convWeight2, &kernelImage, &kernelWeight, ro2->image, ro2->weight,
-                               subMask, maskSource, kernels, polyValues, background, subRegion, useFFT);
+                               subMask, maskSource, kernels, polyValues, background, subRegion, useFFT, true);
             }
 

trunk/psModules/src/imcombine/pmSubtractionEquation.c

@@ -58 +58 @@
         for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex += numKernels) {
             double convPoly = sum * polyValues->data.F64[iyOrder][ixOrder];
-            // XXX Not sure I've got this the right way around
+
+            assert(iIndex < matrix->numRows && j < matrix->numCols);
+
             matrix->data.F64[iIndex][j] = convPoly;
             if (symmetric) {
+
+                assert(iIndex < matrix->numCols && j < matrix->numRows);
+
                 matrix->data.F64[j][iIndex] = convPoly;
             }
@@ -93 +98 @@
                 for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex += numKernels) {
                     double convPoly = sum * iPoly * polyValues->data.F64[jyOrder][jxOrder];
+
+                    assert(iIndex < matrix->numRows && jIndex < matrix->numCols);
+
                     matrix->data.F64[iIndex][jIndex] = convPoly;
                     if (symmetric) {
+
+                        assert(iIndex < matrix->numCols && jIndex < matrix->numRows);
+
                         matrix->data.F64[jIndex][iIndex] = convPoly;
                     }
@@ -139 +150 @@
                             const psArray *convolutions, // Convolutions of source with kernels
                             const psKernel *input, // Input stamp, or NULL
-                            const psKernel *target, // Target stamp, or NULL
                             const psKernel *weight, // Weight stamp
                             const psImage *polyValues, // Spatial polynomial values
@@ -156 +166 @@
     assert(matrix->numCols == numTerms);
     assert(convolutions && convolutions->n == numKernels);
-    assert(target);
     assert(polyValues);
     assert(!normAndBG || input);        // If we want the normalisation and BG, then we need the input image
@@ -168 +177 @@
     // XXX To support higher-order background model than simply constant, the below code needs to be updated.
     if (normAndBG) {
-        // Normalisation and background terms
-        int normIndex, bgIndex;         // Indices in matrix for normalisation and background
-        PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels);
+        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background
 
         for (int i = 0; i < numKernels; i++) {
             psKernel *conv = convolutions->data[i]; // Convolution for i-th element
-            // Normalisation
+
+            // Normalisation-convolution terms
             if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, weight, polyValues, numKernels,
                                          footprint, spatialOrder, true)) {
@@ -181 +190 @@
             }
 
-            // Background
+            // Background-convolution terms
             double sumC = 0.0;          // Sum of the convolution
             for (int y = - footprint; y <= footprint; y++) {
@@ -192 +201 @@
                 return false;
             }
+
             for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
                 for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
-                    matrix->data.F64[index][bgIndex] = matrix->data.F64[bgIndex][index] =
-                        sumC * polyValues->data.F64[yOrder][xOrder];
+                    double value = sumC * polyValues->data.F64[yOrder][xOrder];
+                    matrix->data.F64[index][bgIndex] = value;
+                    matrix->data.F64[bgIndex][index] = value;
                 }
             }
         }
 
-        // Normalisation only and background only terms in the matrix
+        // Background only, normalisation only, and background-normalisation terms
         double sum1 = 0.0;              // Sum of the weighting
         double sumI = 0.0;              // Sum of the input
@@ -206 +217 @@
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
-                float invNoise2 = 1.0 / weight->kernel[y][x];
-                float value = input->kernel[y][x] * invNoise2;
+                double invNoise2 = 1.0 / weight->kernel[y][x];
+                double value = input->kernel[y][x] * invNoise2;
                 sumI += value;
-                sumII += input->kernel[y][x] * value;
+                sumII += value * input->kernel[y][x];
                 sum1 += invNoise2;
             }
+        }
+        if (!isfinite(sumI)) {
+            psTrace("psModules.imcombine", 2, "Bad sumI detected");
+            return false;
         }
         if (!isfinite(sumII)) {
             psTrace("psModules.imcombine", 2, "Bad sumII detected");
-            return false;
-        }
-        if (!isfinite(sumI)) {
-            psTrace("psModules.imcombine", 2, "Bad sumI detected");
             return false;
         }
@@ -239 +250 @@
                             const pmSubtractionKernels *kernels, // Kernel components
                             const psArray *convolutions, // Convolutions of source with kernels
-                            const psKernel *input, // Input stamp, or NULL
-                            const psKernel *target, // Target stamp, or NULL
+                            const psKernel *input, // Input stamp, or NULL if !normAndBG
+                            const psKernel *target, // Target stamp
                             const psKernel *weight, // Weight stamp
                             const psImage *polyValues, // Spatial polynomial values
-                            int footprint // (Half-)Size of stamp
+                            int footprint, // (Half-)Size of stamp
+                            bool normAndBG // Calculate normalisation and background terms?
     )
 {
@@ -250 +262 @@
     int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variation terms
     int bgOrder = kernels->bgOrder;     // Maximum order of background fit
-    int numBackground = PM_SUBTRACTION_POLYTERMS(bgOrder); // Number of background terms
-    int numTerms = numKernels * numSpatial + 1 + numBackground; // Total number of terms
+    int numBackground = normAndBG ? PM_SUBTRACTION_POLYTERMS(bgOrder) : 0; // Number of background terms
+    int numTerms = numKernels * numSpatial + (normAndBG ? 1 + numBackground : 0); // Total number of terms
     assert(vector && vector->n == numTerms);
     assert(convolutions && convolutions->n == numKernels);
-    assert(input);
     assert(target);
     assert(polyValues);
+    assert(!normAndBG || input);       // If we want the normalisation and BG, then we need the input image
 
     // Convolution terms
@@ -278 +290 @@
     }
 
-    // Normalisation and background terms
-    double sumIT = 0.0;         // Sum of the target and input
-    double sumT = 0.0;          // Sum of the target
-    for (int y = - footprint; y <= footprint; y++) {
-        for (int x = - footprint; x <= footprint; x++) {
-            float value = target->kernel[y][x] / weight->kernel[y][x];
-            sumIT += value * input->kernel[y][x];
-            sumT += value;
-        }
-    }
-    if (!isfinite(sumIT)) {
-        psTrace("psModules.imcombine", 2, "Bad sumIT detected");
-        return false;
-    }
-    if (!isfinite(sumT)) {
-        psTrace("psModules.imcombine", 2, "Bad sumI detected");
-        return false;
-    }
-
-    int normIndex, bgIndex;             // Indices in vector for normalisation and background terms
-    PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels);
-    vector->data.F64[normIndex] = sumIT;
-    vector->data.F64[bgIndex] = sumT;
+    if (normAndBG) {
+        // Background terms
+        double sumT = 0.0;              // Sum of the target
+        double sumIT = 0.0;             // Sum of the input-target product
+        for (int y = - footprint; y <= footprint; y++) {
+            for (int x = - footprint; x <= footprint; x++) {
+                float value = target->kernel[y][x] / weight->kernel[y][x];
+                sumIT += value * input->kernel[y][x];
+                sumT += value;
+            }
+        }
+        if (!isfinite(sumT)) {
+            psTrace("psModules.imcombine", 2, "Bad sumI detected");
+            return false;
+        }
+        if (!isfinite(sumIT)) {
+            psTrace("psModules.imcombine", 2, "Bad sumIT detected");
+            return false;
+        }
+
+        int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation term
+        vector->data.F64[normIndex] = sumIT;
+        int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background term
+        vector->data.F64[bgIndex] = sumT;
+    }
 
     return true;
@@ -383 +397 @@
     for (int yOrder = 0; yOrder <= order; yOrder++) {
         for (int xOrder = 0; xOrder <= order - yOrder; xOrder++, index += step) {
+
+            assert(index < coeff->n);
+
             sum += coeff->data.F64[index] * polyValues->data.F64[yOrder][xOrder];
         }
@@ -468 +485 @@
         }
 
+#ifdef TESTING
+        for (int j = 0; j < numKernels; j++) {
+            if (stamp->convolutions1) {
+                psString convName = NULL;
+                psStringAppend(&convName, "conv1_%03d_%03d.fits", i, 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", i, j);
+                psFits *fits = psFitsOpen(convName, "w");
+                psFree(convName);
+                psKernel *conv = stamp->convolutions2->data[j];
+                psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
+                psFitsClose(fits);
+            }
+        }
+#endif
+
         polyValues = p_pmSubtractionPolynomial(polyValues, spatialOrder, stamp->xNorm, stamp->yNorm);
 
-        stamp->vector = psVectorAlloc(numParams, PS_TYPE_F64);
+
+            stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+            stamp->vector1 = psVectorAlloc(numParams, PS_TYPE_F64);
 #ifdef TESTING
-        psVectorInit(stamp->vector, NAN);
+            psImageInit(stamp->matrix1, NAN);
+            psVectorInit(stamp->vector1, NAN);
 #endif
 
@@ -479 +523 @@
           case PM_SUBTRACTION_MODE_TARGET:
           case PM_SUBTRACTION_MODE_1:
-            stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
-#ifdef TESTING
-            psImageInit(stamp->matrix1, NAN);
-#endif
             status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
-                                     stamp->image2, stamp->weight, polyValues, footprint, true);
-            status &= calculateVector(stamp->vector, kernels, stamp->convolutions1, stamp->image1,
-                                      stamp->image2, stamp->weight, polyValues, footprint);
+                                     stamp->weight, polyValues, footprint, true);
+            status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1,
+                                      stamp->image2, stamp->weight, polyValues, footprint, true);
             break;
           case PM_SUBTRACTION_MODE_2:
-            stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
-#ifdef TESTING
-            psImageInit(stamp->matrix1, NAN);
-#endif
             status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2,
-                                     stamp->image1, stamp->weight, polyValues, footprint, true);
-            status &= calculateVector(stamp->vector, kernels, stamp->convolutions2, stamp->image2,
-                                      stamp->image1, stamp->weight, polyValues, footprint);
+                                     stamp->weight, polyValues, footprint, true);
+            status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2,
+                                      stamp->image1, stamp->weight, polyValues, footprint, true);
             break;
           case PM_SUBTRACTION_MODE_DUAL:
-            stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64);
             stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64);
             stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64);
+            stamp->vector2 = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64);
 #ifdef TESTING
-            psImageInit(stamp->matrix1, NAN);
             psImageInit(stamp->matrix2, NAN);
             psImageInit(stamp->matrixX, NAN);
-#endif
-            status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
-                                     stamp->image2, stamp->weight, polyValues, footprint, true);
-            status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, stamp->image2,
-                                      stamp->image1, stamp->weight, polyValues, footprint, false);
-            // XXX Not sure I'm passing the correct stamp->image to calculateMatrixCross:
+            psVectorInit(stamp->vector2, NAN);
+#endif
+            status  = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1,
+                                      stamp->weight, polyValues, footprint, true);
+            status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL,
+                                      stamp->weight, polyValues, footprint, false);
             status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1,
-                                           stamp->convolutions2, stamp->image2, stamp->weight, polyValues,
+                                           stamp->convolutions2, stamp->image1, stamp->weight, polyValues,
                                            footprint);
-            status &= calculateVector(stamp->vector, kernels, stamp->convolutions1, stamp->image1,
-                                      stamp->image2, stamp->weight, polyValues, footprint);
+            status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1,
+                                      stamp->image2, stamp->weight, polyValues, footprint, true);
+            status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL,
+                                      stamp->image2, stamp->weight, polyValues, footprint, false);
             break;
           default:
@@ -540 +577 @@
 
             matrixName = NULL;
-            psStringAppend(&matrixName, "vector_%d.fits", i);
-            psImage *dummy = psImageAlloc(stamp->vector->n, 1, PS_TYPE_F64);
-            memcpy(dummy->data.F64[0], stamp->vector->data.F64,
-                   PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector->n);
+            psStringAppend(&matrixName, "vector1_%d.fits", i);
+            psImage *dummy = psImageAlloc(stamp->vector1->n, 1, PS_TYPE_F64);
+            memcpy(dummy->data.F64[0], stamp->vector1->data.F64,
+                   PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector1->n);
             matrixFile = psFitsOpen(matrixName, "w");
             psFree(matrixName);
@@ -550 +587 @@
             psFitsClose(matrixFile);
 
+            if (stamp->vector2) {
+                matrixName = NULL;
+                psStringAppend(&matrixName, "vector2_%d.fits", i);
+                dummy = psImageAlloc(stamp->vector2->n, 1, PS_TYPE_F64);
+                memcpy(dummy->data.F64[0], stamp->vector2->data.F64,
+                       PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector2->n);
+                matrixFile = psFitsOpen(matrixName, "w");
+                psFree(matrixName);
+                psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL);
+                psFree(dummy);
+                psFitsClose(matrixFile);
+            }
 
             if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
@@ -575 +624 @@
 }
 
-bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps,
-                                float tolerance)
+bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps)
 {
     PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
-    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        PS_ASSERT_FLOAT_LARGER_THAN(tolerance, 0.0, false);
-    }
 
     // Check inputs
@@ -594 +639 @@
         }
 
+        PS_ASSERT_VECTOR_NON_NULL(stamp->vector1, false);
+        if (numParams == -1) {
+            numParams = stamp->vector1->n;
+        }
+        PS_ASSERT_VECTOR_SIZE(stamp->vector1, (long)numParams, false);
+        PS_ASSERT_VECTOR_TYPE(stamp->vector1, PS_TYPE_F64, false);
+        PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false);
+        PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false);
+        PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false);
+
         if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-            PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false);
             PS_ASSERT_IMAGE_NON_NULL(stamp->matrix2, false);
             PS_ASSERT_IMAGE_NON_NULL(stamp->matrixX, false);
-            if (numParams == -1) {
-                numParams = stamp->matrix1->numCols;
+            if (numParams2 == 0) {
                 numParams2 = stamp->matrix2->numCols;
             }
-            PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false);
             PS_ASSERT_IMAGE_SIZE(stamp->matrix2, numParams2, numParams2, false);
             PS_ASSERT_IMAGE_SIZE(stamp->matrixX, numParams, numParams2, false);
-            PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false);
             PS_ASSERT_IMAGE_TYPE(stamp->matrix2, PS_TYPE_F64, false);
             PS_ASSERT_IMAGE_TYPE(stamp->matrixX, PS_TYPE_F64, false);
-            PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false);
-            PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false);
-            PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false);
-        } else {
-            PS_ASSERT_VECTOR_NON_NULL(stamp->vector, false);
-            if (numParams == -1) {
-                numParams = stamp->vector->n;
-            }
-            PS_ASSERT_VECTOR_SIZE(stamp->vector, (long)numParams, false);
-            PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, false);
-            PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false);
-            PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false);
-            PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false);
+            PS_ASSERT_VECTOR_NON_NULL(stamp->vector2, false);
+            PS_ASSERT_VECTOR_SIZE(stamp->vector2, (long)numParams2, false);
+            PS_ASSERT_VECTOR_TYPE(stamp->vector2, PS_TYPE_F64, false);
         }
     }
@@ -638 +679 @@
             if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
                 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1);
-                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
+                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1);
             }
         }
@@ -667 +708 @@
         psImage *sumMatrix2 = psImageAlloc(numParams2, numParams2, PS_TYPE_F64);
         psImage *sumMatrixX = psImageAlloc(numParams, numParams2, PS_TYPE_F64);
-        psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64);
+        psVector *sumVector1 = psVectorAlloc(numParams, PS_TYPE_F64);
+        psVector *sumVector2 = psVectorAlloc(numParams, PS_TYPE_F64);
         psImageInit(sumMatrix1, 0.0);
         psImageInit(sumMatrix2, 0.0);
         psImageInit(sumMatrixX, 0.0);
-        psVectorInit(sumVector, 0.0);
+        psVectorInit(sumVector1, 0.0);
+        psVectorInit(sumVector2, 0.0);
 
         for (int i = 0; i < stamps->num; i++) {
@@ -679 +722 @@
                 (void)psBinaryOp(sumMatrix2, sumMatrix2, "+", stamp->matrix2);
                 (void)psBinaryOp(sumMatrixX, sumMatrixX, "+", stamp->matrixX);
-                (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
-            }
-        }
-
-        psVector *permutation1 = NULL, *permutation2 = NULL; // Permutation vectors for LUD
-        psImage *lu1 = psMatrixLUD(NULL, &permutation1, sumMatrix1); // LUD of matrix 1
-        psImage *lu2 = psMatrixLUD(NULL, &permutation2, sumMatrix2); // LUD of matrix 2
-        psFree(sumMatrix1);
-        psFree(sumMatrix2);
-        if (!lu1 || !lu2) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to generate LU decomposed matrices.");
-            psFree(permutation1);
-            psFree(permutation2);
-            psFree(lu1);
-            psFree(lu2);
-            psFree(sumMatrixX);
-            psFree(sumVector);
-            return NULL;
-        }
-
-        psVector *lastSol2 = psVectorAlloc(numParams2, PS_TYPE_F64); // Value of sol2 in last iteration
-        psVectorInit(lastSol2, 0.0);
-
-        double diff = NAN;              // Difference between iterations
-        psVector *vector = NULL;        // RHS for equations
-        int iter = 0;                   // Iteration number
-        do {
-            if (!vector) {
-                // Start with traditional Alard-Lupton solution
-                vector = psVectorCopy(NULL, sumVector, PS_TYPE_F64);
-            } else {
-                // Generate vector for RHS of equation 1
-                vector->n = numParams;
-                psVector *sol2 = kernels->solution2;
-                assert(sol2);
-                for (int i = 0; i < numParams; i++) {
-                    double value = 0.0;             // Value of vector element
-                    for (int j = 0; j < numParams2; j++) {
-                        // XXX Not sure I've got this the right way around
-                        value += sol2->data.F64[j] * sumMatrixX->data.F64[j][i];
-                    }
-                    vector->data.F64[i] = value;
-                }
-            }
-
-            kernels->solution1 = psMatrixLUSolve(kernels->solution1, lu1, vector, permutation1);
-            if (!kernels->solution1) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to solve initial matrix equation.");
-                psFree(permutation1);
-                psFree(permutation2);
-                psFree(lu1);
-                psFree(lu2);
-                psFree(sumMatrixX);
-                psFree(sumVector);
-                psFree(vector);
-                psFree(lastSol2);
-                return NULL;
-            }
-
-            // Generate vector for RHS of equation 2
-            vector->n = numParams2;
-            psVector *sol1 = kernels->solution1; // Solution of first equation
-            for (int i = 0; i < numParams2; i++) {
-                double value = 0.0;             // Value of vector element
-                for (int j = 0; j < numParams; j++) {
-                    // XXX Not sure I've got this the right way around
-                    value += sol1->data.F64[j] * sumMatrixX->data.F64[i][j];
-                }
-                vector->data.F64[i] = value;
-            }
-
-            kernels->solution2 = psMatrixLUSolve(kernels->solution2, lu2, vector, permutation2);
-            if (!kernels->solution2) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to solve matrix equation 2.");
-                psFree(permutation1);
-                psFree(permutation2);
-                psFree(lu1);
-                psFree(lu2);
-                psFree(sumMatrixX);
-                psFree(sumVector);
-                psFree(vector);
-                psFree(lastSol2);
-                return NULL;
-            }
-
-            // Get difference produced by iteration
-            diff = 0.0;
+                (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1);
+                (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2);
+            }
+        }
+
 #if 0
-            psVector *sol2 = kernels->solution2; // Solution of second equation
-            for (int i = 0; i < numParams2; i++) {
-                diff += fabs(sol2->data.F64[i] - lastSol2->data.F64[i]);
-            }
-            lastSol2 = psVectorCopy(lastSol2, sol2, PS_TYPE_F64);
-            psTrace("psModules.imcombine", 4, "Solution iteration %d: %lf\n", iter, diff);
-#endif
-
-            psVectorInit(kernels->solution2, 0.0);
-
-            iter++;
-        } while (diff > tolerance);
-
-        psFree(permutation1);
-        psFree(permutation2);
-        psFree(lu1);
-        psFree(lu2);
-        psFree(sumMatrixX);
-        psFree(sumVector);
-        psFree(vector);
-        psFree(lastSol2);
+        // Apply weighting to maximise the difference between solution coefficients for the same functions
+        double fudge = PS_SQR(2 * stamps->footprint + 1); // Fudge factor
+        for (int i = 0; i < kernels->num; i++) {
+            sumMatrix1->data.F64[i][i] -= fudge;
+            sumMatrix2->data.F64[i][i] += fudge;
+        }
+#endif
+
+        // Pure matrix operations
+
+        // A * a = Ct * b + d
+        // C * a = B  * b + e
+        //
+        // a = (Ct * Bi * C - A)i (Ct * Bi * e - d)
+        // b = Bi * (C * a - e)
+        psVector *a = psVectorRecycle(kernels->solution1, numParams, PS_TYPE_F64);
+        psVector *b = psVectorRecycle(kernels->solution2, numParams2, PS_TYPE_F64);
+#ifdef TESTING
+        psVectorInit(a, NAN);
+        psVectorInit(b, NAN);
+#endif
+        psImage *A = sumMatrix1;
+        psImage *B = sumMatrix2;
+        psImage *C = sumMatrixX;
+        psVector *d = sumVector1;
+        psVector *e = sumVector2;
+
+        assert(a->n == numParams);
+        assert(b->n == numParams2);
+        assert(A->numRows == numParams && A->numCols == numParams);
+        assert(B->numRows == numParams2 && B->numCols == numParams2);
+        assert(C->numRows == numParams2 && C->numCols == numParams);
+        assert(d->n == numParams);
+        assert(e->n == numParams2);
+
+        psImage *Bi = psMatrixInvert(NULL, B, NULL);
+        assert(Bi->numRows == numParams2 && Bi->numCols == numParams2);
+        psImage *Ct = psMatrixTranspose(NULL, C);
+        assert(Ct->numRows == numParams && Ct->numCols == numParams2);
+
+        psImage *BiC = psMatrixMultiply(NULL, Bi, C);
+        assert(BiC->numRows == numParams2 && BiC->numCols == numParams);
+        psImage *CtBi = psMatrixMultiply(NULL, Ct, Bi);
+        assert(CtBi->numRows == numParams && CtBi->numCols == numParams2);
+
+        psImage *CtBiC = psMatrixMultiply(NULL, Ct, BiC);
+        assert(CtBiC->numRows == numParams && CtBiC->numCols == numParams);
+
+        psImage *F = (psImage*)psBinaryOp(NULL, CtBiC, "-", A);
+        assert(F->numRows == numParams && F->numCols == numParams);
+        float det = NAN;
+        psImage *Fi = psMatrixInvert(NULL, F, &det);
+        assert(Fi->numRows == numParams && Fi->numCols == numParams);
+        psTrace("psModules.imcombine", 4, "Determinant of F: %f\n", det);
+
+        psVector *g = psVectorAlloc(numParams, PS_TYPE_F64);
+#ifdef TESTING
+        psVectorInit(g, NAN);
+#endif
+        assert(CtBi->numRows == numParams && CtBi->numCols == numParams2);
+        assert(e->n == numParams2);
+        assert(d->n == numParams);
+        for (int i = 0; i < numParams; i++) {
+            double value = 0.0;
+            for (int j = 0; j < numParams2; j++) {
+                value += CtBi->data.F64[i][j] * e->data.F64[j];
+            }
+            g->data.F64[i] = value - d->data.F64[i];
+        }
+
+        assert(Fi->numRows == numParams && Fi->numCols == numParams);
+        assert(g->n == numParams);
+        for (int i = 0; i < numParams; i++) {
+            double value = 0.0;
+            for (int j = 0; j < numParams; j++) {
+                value += Fi->data.F64[i][j] * g->data.F64[j];
+            }
+            a->data.F64[i] = value;
+        }
+
+        psVector *h = psVectorAlloc(numParams2, PS_TYPE_F64);
+#ifdef TESTING
+        psVectorInit(h, NAN);
+#endif
+        assert(C->numRows == numParams2 && C->numCols == numParams);
+        assert(a->n == numParams);
+        assert(e->n == numParams2);
+        for (int i = 0; i < numParams2; i++) {
+            double value = 0.0;
+            for (int j = 0; j < numParams; j++) {
+                value += C->data.F64[i][j] * a->data.F64[j];
+            }
+            h->data.F64[i] = value - e->data.F64[i];
+        }
+
+        assert(Bi->numRows == numParams2 && Bi->numCols == numParams2);
+        assert(h->n == numParams2);
+        for (int i = 0; i < numParams2; i++) {
+            double value = 0.0;
+            for (int j = 0; j < numParams2; j++) {
+                value += Bi->data.F64[i][j] * h->data.F64[j];
+            }
+            b->data.F64[i] = value;
+        }
+
+
+#if 0
+        for (int i = 0; i < numParams; i++) {
+            double aVal1 = 0.0, bVal1 = 0.0;
+            for (int j = 0; j < numParams2; j++) {
+                aVal1 += A->data.F64[i][j] * a->data.F64[j];
+                bVal1 += Ct->data.F64[i][j] * b->data.F64[j];
+            }
+            bVal1 += d->data.F64[i];
+            for (int j = numParams2; j < numParams; j++) {
+                aVal1 += A->data.F64[i][j] * a->data.F64[j];
+            }
+            printf("%d: %lf\n", i, aVal1 - bVal1);
+        }
+
+        for (int i = 0; i < numParams2; i++) {
+            double aVal2 = 0.0, bVal2 = 0.0;
+            for (int j = 0; j < numParams2; j++) {
+                aVal2 += C->data.F64[i][j] * a->data.F64[j];
+                bVal2 += B->data.F64[i][j] * b->data.F64[j];
+            }
+            bVal2 += e->data.F64[i];
+            for (int j = numParams2; j < numParams; j++) {
+                aVal2 += C->data.F64[i][j] * a->data.F64[j];
+            }
+            printf("%d: %lf\n", i, aVal2 - bVal2);
+        }
+#endif
+
+        {
+            psFits *fits = psFitsOpen("sumMatrix1.fits", "w");
+            psFitsWriteImage(fits, NULL, sumMatrix1, 0, NULL);
+            psFitsClose(fits);
+        }
+        {
+            psFits *fits = psFitsOpen("sumMatrix2.fits", "w");
+            psFitsWriteImage(fits, NULL, sumMatrix2, 0, NULL);
+            psFitsClose(fits);
+        }
+        {
+            psFits *fits = psFitsOpen("sumMatrixX.fits", "w");
+            psFitsWriteImage(fits, NULL, sumMatrixX, 0, NULL);
+            psFitsClose(fits);
+        }
+        {
+            psFits *fits = psFitsOpen("sumFinverse.fits", "w");
+            psFitsWriteImage(fits, NULL, Fi, 0, NULL);
+            psFitsClose(fits);
+        }
+
+
+        kernels->solution1 = a;
+        kernels->solution2 = b;
+
+        // XXXXX Free temporary matrices and vectors
+
     }
 
@@ -817 +918 @@
     int numKernels = kernels->num;      // Number of kernels
 
-    psVector *coeff1 = psVectorAlloc(numKernels, PS_TYPE_F64); // Coefficients for first image
-    psVector *coeff2 = kernels->mode == PM_SUBTRACTION_MODE_DUAL ?
-        psVectorAlloc(numKernels, PS_TYPE_F64) : NULL;            // Coefficients for second image
     psImage *polyValues = NULL;         // Polynomial values
     psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
@@ -832 +930 @@
         // Calculate coefficients of the kernel basis functions
         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
-        for (int j = 0; j < numKernels; j++) {
-            coeff1->data.F64[j] = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false);
-        }
         double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
@@ -863 +958 @@
             for (int j = 0; j < numKernels; j++) {
                 psKernel *convolution = convolutions->data[j]; // Convolution
-                double coefficient = coeff1->data.F64[j]; // Coefficient
+                double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j,
+                                                                  false); // Coefficient
                 for (int y = - footprint; y <= footprint; y++) {
                     for (int x = - footprint; x <= footprint; x++) {
@@ -880 +976 @@
             psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image
             psKernel *image1 = stamp->image1; // The first image
-
-            for (int j = 0; j < numKernels; j++) {
-                coeff2->data.F64[j] = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true);
-            }
+            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 coefficient1 = coeff1->data.F64[j]; // Coefficient for convolution 1
-                double coefficient2 = coeff2->data.F64[j]; // Coefficient for convolution 2
+                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] += conv1->kernel[y][x] * coefficient1 +
-                            conv2->kernel[y][x] * coefficient2;
+                        residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 - conv1->kernel[y][x] * coeff1;
                     }
                 }
@@ -901 +992 @@
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
-                    residual->kernel[y][x] += image1->kernel[y][x] * norm - background;
+                    residual->kernel[y][x] += image2->kernel[y][x] - background - image1->kernel[y][x] * norm;
                 }
             }
@@ -909 +1000 @@
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
-                deviation += PS_SQR(residual->kernel[y][x]) / weight->kernel[y][x];
+                double dev = PS_SQR(residual->kernel[y][x]) / weight->kernel[y][x];
+                deviation += dev;
+#ifdef TESTING
+                residual->kernel[y][x] = dev;
+#endif
             }
         }
@@ -932 +1027 @@
             psFitsClose(fits);
         }
+        {
+            psString filename = NULL;
+            psStringAppend(&filename, "stamp_image_%03d.fits", i);
+            psFits *fits = psFitsOpen(filename, "w");
+            psFree(filename);
+            psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL);
+            psFitsClose(fits);
+        }
+        {
+            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
 
@@ -937 +1048 @@
     psFree(residual);
     psFree(polyValues);
-    psFree(coeff1);
 
     return deviations;

trunk/psModules/src/imcombine/pmSubtractionEquation.h

@@ -12 +12 @@
 /// Solve the least-squares equation to match the image quality
 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters
-                                const pmSubtractionStampList *stamps, ///< Stamps
-                                float tolerance ///< Tolerance for iteration (dual mode only)
+                                const pmSubtractionStampList *stamps ///< Stamps
     );
 

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -136 +136 @@
                         }
                     }
+                    preCalc->kernel[0][0] -= 1.0;
                 }
 
@@ -216 +217 @@
     }
 
-    // Subtract normalisation from all the others to maintain flux scaling
-    if (spatialOrder > 0) {
-        int numGaussians = fwhms->n;       // Number of Gaussians
-        for (int i = 0, index = 0; i < numGaussians; i++) {
-            for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
-                for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
-                    if (uOrder % 2 == 0 && vOrder % 2 == 0) {
-                        psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract
-                        kernel->kernel[0][0] -= 1.0;
-                    }
-                }
-            }
-        }
-    }
-
     if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) {
         for (int i = 0; i < kernels->num; i++) {
@@ -246 +232 @@
                 }
             }
-            psTrace("psModules.imcombine.kernel", 10, "Kernel %d sum: %lf\n", i, sum);
+            psTrace("psModules.imcombine.kernel", 10, "Kernel %d sum: %le\n", i, sum);
         }
     }
@@ -447 +433 @@
     psStringPrepend(&kernels->description, "GUNK=");
     psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder);
-
-    // Subtract unity from the kernels to maintain photometric flux scaling
-    int numGaussians = fwhms->n;       // Number of Gaussians
-    if (spatialOrder > 0) {
-        for (int i = 0, index = 0; i < numGaussians; i++) {
-            for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
-                for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
-                    if (uOrder % 2 == 0 && vOrder % 2 == 0) {
-                        psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract
-                        kernel->kernel[0][0] -= 1.0;
-                    }
-                }
-            }
-        }
-    }
 
     int numISIS = kernels->num;         // Number of ISIS kernels

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -18 +18 @@
 #include "pmSubtractionMatch.h"
 
-#define TOL 1.0e-3                      // Tolerance for subtraction solution
 #define KERNEL_MOSAIC 2                 // Half-number of kernel instances in the mosaic image
 
@@ -24 +23 @@
 
 
-//#define TESTING
+#define TESTING
+//#define TESTING_MEMORY
 
 // Output memory usage information
 static void memCheck(const char *where)
 {
-#ifdef TESTING
+#ifdef TESTING_MEMORY
     psMemBlock **leaks = NULL;
     int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true);
@@ -216 +216 @@
     psString regionString = NULL;       // String for region
     pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF
-    psVector *solution = NULL;          // Solution to match PSF
     pmSubtractionKernels *kernels = NULL; // Kernel basis functions
     int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
@@ -261 +260 @@
                                                            stampSpacing, mode);
             } else if (stampsName && strlen(stampsName) > 0) {
-                stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, footprint,
+                stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, footprint,
                                                         stampSpacing, mode);
             }
@@ -331 +330 @@
                 psTrace("psModules.imcombine", 3, "Solving equation...\n");
 
-                if (!pmSubtractionSolveEquation(kernels, stamps, TOL)) {
+                if (!pmSubtractionSolveEquation(kernels, stamps)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
@@ -360 +359 @@
             if (numRejected > 0) {
                 psTrace("psModules.imcombine", 3, "Solving equation...\n");
-                if (!pmSubtractionSolveEquation(kernels, stamps, TOL)) {
+                if (!pmSubtractionSolveEquation(kernels, stamps)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
@@ -413 +412 @@
             }
 
-#ifdef TESTING
+#if 0
             {
                 // Generate images of the kernel components
@@ -442 +441 @@
             kernels = NULL;
 
+#if 0
             // Put the solution on the metadata
             {
@@ -460 +460 @@
                 }
             }
-
-            psFree(solution);
-            solution = NULL;
+#endif
 
             // There is data in the readout now
@@ -496 +494 @@
 
 
+#ifdef TESTING
+    {
+        psFits *fits = psFitsOpen("convolved1.fits", "w");
+        psFitsWriteImage(fits, NULL, conv1->image, 0, NULL);
+        psFitsClose(fits);
+
+        if (mode == PM_SUBTRACTION_MODE_DUAL) {
+            psFits *fits = psFitsOpen("convolved2.fits", "w");
+            psFitsWriteImage(fits, NULL, conv2->image, 0, NULL);
+            psFitsClose(fits);
+        }
+    }
+#endif
+
     return true;
 
@@ -504 +516 @@
     psFree(kernels);
     psFree(stamps);
-    psFree(solution);
     psFree(weight);
     return false;

trunk/psModules/src/imcombine/pmSubtractionStamps.h

@@ -63 +63 @@
     psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component, or NULL
     psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component, or NULL
-    psImage *matrix1, *matrix2;         ///< Matrices for each image, or NULL
+    psImage *matrix1, *matrix2;         ///< Least-squares matrices for each image, or NULL
     psImage *matrixX;                   ///< Cross-matrix (for mode DUAL), or NULL
-    psVector *vector;                   ///< Associated vector (when mode not DUAL), or NULL
+    psVector *vector1, *vector2;        ///< Least-squares vectors for each image, or NULL
     pmSubtractionStampStatus status;    ///< Status of stamp
 } pmSubtractionStamp;
@@ -108 +108 @@
 pmSubtractionStampList *pmSubtractionStampsSetFromFile(
     const char *filename,               ///< Filename of file containing x,y (or x,y,flux) on each line
+    const psImage *image,               ///< Image for flux of stamp
     const psImage *subMask,             ///< Mask, or NULL
     const psRegion *region,             ///< Region to search, or NULL