Changeset 19765 for trunk/psModules/src/imcombine/pmSubtractionEquation.c

Timestamp:

Sep 26, 2008, 4:03:24 PM (18 years ago)

Author:

Paul Price

Message:

Two major changes. 1. Optimising ISIS kernels: instead of doing an FFT convolution, they're separable. 2. Trying to make the convolved weight map account for systematic errors in the convolution kernel by assuming a relative error in each kernel pixel (currently hardwired with a #define, but we can turn it into a recipe value later, or get an estimate from the data) and propagating that into the variance map. The effect is to increase the variance near the bright stars, which is hopefully sufficient to decrease the number of false sources detected in bad subtractions.

File:

• trunk/psModules/src/imcombine/pmSubtractionEquation.c (modified) (5 diffs)

trunk/psModules/src/imcombine/pmSubtractionEquation.c

@@ -758 +758 @@
         calculatePenalty(sumVector, kernels);
 
+#ifdef TESTING
+        {
+            psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
+            psFits *fits = psFitsOpen("matrixInv.fits", "w");
+            psFitsWriteImage(fits, NULL, inverse, 0, NULL);
+            psFitsClose(fits);
+            psFree(inverse);
+        }
+        {
+            psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
+            psImage *Xt = psMatrixTranspose(NULL, X);
+            psImage *XtX = psMatrixMultiply(NULL, Xt, X);
+            psFits *fits = psFitsOpen("matrixErr.fits", "w");
+            psFitsWriteImage(fits, NULL, XtX, 0, NULL);
+            psFitsClose(fits);
+            psFree(X);
+            psFree(Xt);
+            psFree(XtX);
+        }
+#endif
+
         psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
         psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix);
@@ -769 +790 @@
         }
         kernels->solution1 = psMatrixLUSolve(kernels->solution1, luMatrix, sumVector, permutation);
+
         psFree(sumVector);
         psFree(luMatrix);
@@ -806 +828 @@
         calculatePenalty(sumVector2, kernels);
 
-#if 0
-        // 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
 
@@ -941 +954 @@
 #endif
 
+#ifdef TESTING
         {
             psFits *fits = psFitsOpen("sumMatrix1.fits", "w");
@@ -961 +975 @@
             psFitsClose(fits);
         }
-
+#endif
 
         kernels->solution1 = a;