Changeset 26340

Timestamp:

Dec 5, 2009, 5:02:29 AM (16 years ago)

Author:

eugene

Message:

allow either separate or simultaneous coeff calculation (norm vs kernels) -- compile time option; report norm and background

File:

• branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (modified) (9 diffs)

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

@@ -104 +104 @@
 
             // Spatial variation of kernel coeffs
-            if (mode | PM_SUBTRACTION_EQUATION_KERNELS) {
+            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) {
@@ -147 +147 @@
             double normTerm = sumRC * poly[iTerm];
             double bgTerm = sumC * poly[iTerm];
-            if ((mode | PM_SUBTRACTION_EQUATION_NORM) && (mode | PM_SUBTRACTION_EQUATION_KERNELS)) {
+            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)) {
+            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) {
+            if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
                 vector->data.F64[iIndex] = sumIC * poly[iTerm];
+                if (!(mode & PM_SUBTRACTION_EQUATION_NORM)) {
+                    // subtract norm * sumRC * poly[iTerm]
+                    psAssert (kernels->solution1, "programming error: define solution first!");
+                    int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+                    double norm = fabs(kernels->solution1->data.F64[normIndex]);  // Normalisation
+                    vector->data.F64[iIndex] -= norm * normTerm;
+                }
             }
         }
@@ -196 +203 @@
         }
     }
-    if (mode | PM_SUBTRACTION_EQUATION_NORM) {
+    if (mode & PM_SUBTRACTION_EQUATION_NORM) {
         matrix->data.F64[normIndex][normIndex] = sumRR;
         vector->data.F64[normIndex] = sumIR;
-    }
-    if (mode | PM_SUBTRACTION_EQUATION_BG) {
+        // subtract sum over kernels * kernel solution
+    }
+    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)) {
+    if ((mode & PM_SUBTRACTION_EQUATION_NORM) && (mode & PM_SUBTRACTION_EQUATION_BG)) {
         matrix->data.F64[normIndex][bgIndex] = sumR;
         matrix->data.F64[bgIndex][normIndex] = sumR;
@@ -953 +961 @@
 #endif
 
+        // XXX test: save the matrix A and vector b:
+        {
+            psFits *fits = psFitsOpen("matrix.fits", "w");
+            psFitsWriteImage(fits, NULL, sumMatrix, 0, NULL);
+            psFitsClose(fits);
+
+            FILE *f = fopen ("vector.dat", "w");
+            int fd = fileno(f);
+            p_psVectorPrint (fd, sumVector, "B");
+            fclose (f);
+        }
+
         psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
         psImage *luMatrix = psMatrixLUDecomposition(NULL, &permutation, sumMatrix);
@@ -979 +999 @@
 
         // only update the solutions that we chose to calculate:
-        if (mode | PM_SUBTRACTION_EQUATION_NORM) {
+        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) {
+        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) {
+        if (mode & PM_SUBTRACTION_EQUATION_KERNELS) {
             int numKernels = kernels->num;
             int spatialOrder = kernels->spatialOrder;       // Order of spatial variation
@@ -1225 +1245 @@
      }
 
-    pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const
+    // pmSubtractionVisualPlotLeastSquares((pmSubtractionStampList *) stamps); //casting away const
     return true;
 }
@@ -1241 +1261 @@
     double devNorm = 1.0 / (double)numPixels; // Normalisation for deviations
     int numKernels = kernels->num;      // Number of kernels
-    double norm = NAN;
 
     psImage *polyValues = NULL;         // Polynomial values
@@ -1258 +1277 @@
         // Calculate coefficients of the kernel basis functions
         polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm);
-        norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
         double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
 
@@ -1404 +1423 @@
     
     }
+
+    // calculate and report the normalization and background for the image center
+    { 
+        polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, 0.0, 0.0);
+        double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation
+        double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background
+        psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f, background: %f", norm, background);
+    }
+
+    pmSubtractionVisualShowFit();
+    pmSubtractionVisualPlotFit(kernels);
+
     psFree(residual);
     psFree(polyValues);
 
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "normalization: %f", norm);
-    pmSubtractionVisualShowFit();
-    pmSubtractionVisualPlotFit(kernels);
-
     return deviations;
 }