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

Timestamp:

Jan 26, 2011, 5:21:42 PM (16 years ago)

Author:

eugene

Message:

report the diff chisq as well as the residual chisq; truncate SVD at dynamic range of 1e10; fix the scaling process

File:

• branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionEquation.c (modified) (16 diffs)

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

@@ -1072 +1072 @@
         // SINGLE solution
 # if (1)
-        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10);
         invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
 # endif
@@ -1176 +1176 @@
         // DUAL solution
 # if (1)
-        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-10);
         invMatrix = psMatrixInvert(NULL, sumMatrix, NULL);
 # endif
@@ -1306 +1306 @@
 
 // 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, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *residual, psKernel *weight, psKernel *window) {
+bool pmSubtractionChisqStats(psVector *fluxesVector, psVector *chisqDVector, psVector *chisqRVector, psVector *momentVector, psVector *stampMask, psKernel *convolved1, psKernel *convolved2, psKernel *difference, psKernel *residual, psKernel *weight, psKernel *window) {
 
 # ifndef USE_WEIGHT
@@ -1316 +1316 @@
 
     int npix = 0;
-    float chisq = 0;
+    float chisqR = 0;
+    float chisqD = 0;
 
     // get the chisq
     for (int y = residual->yMin; y <= residual->yMax; y++) {
         for (int x = residual->xMin; x <= residual->xMax; x++) {
-            float value = PS_SQR(residual->kernel[y][x]);
+            float valueR = PS_SQR(residual->kernel[y][x]);
             if (weight) {
-                value *= weight->kernel[y][x];
+                valueR *= 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];
+            // XXX NOTE: do NOT apply the window to the chisq portions of the calculation (that would bias the chisq)
+            chisqR += valueR;
+
+            float valueD = PS_SQR(difference->kernel[y][x]);
+            if (weight) {
+                valueD *= weight->kernel[y][x];
             }
-            chisq += value;
+            chisqD += valueD;
             npix ++;
         }
     }
-    psVectorAppend(chisqVector, chisq / npix);
+    psVectorAppend(chisqRVector, chisqR / npix);
+    psVectorAppend(chisqDVector, chisqD / npix);
 
     float value1 = 0;
@@ -1430 +1435 @@
     int Nelem = fluxesVector->n - 1;
     bool valid = true;
-    valid &= isfinite(chisqVector->data.F32[Nelem]);
+    valid &= isfinite(chisqRVector->data.F32[Nelem]);
     valid &= isfinite(fluxesVector->data.F32[Nelem]);
     valid &= isfinite(momentVector->data.F32[Nelem]);
@@ -1453 +1458 @@
     // XXX need to save these somewhere
     psVector *fluxes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
-    psVector *chisq = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *chisqD = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *chisqR = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     psVector *moments = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
     psVector *stampMask = psVectorAllocEmpty(stamps->num, PS_TYPE_VECTOR_MASK);
@@ -1464 +1470 @@
     // storage for the image (convolved2 is not used in SINGLE mode)
     psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
+    psKernel *difference = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
     psKernel *convolved1 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
     psKernel *convolved2 = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image
@@ -1474 +1481 @@
             psVectorAppend(moments, NAN);
             psVectorAppend(fluxes, NAN);
-            psVectorAppend(chisq, NAN);
+            psVectorAppend(chisqD, NAN);
+            psVectorAppend(chisqR, NAN);
             psVectorAppend(stampMask, 0x01);
             continue;
@@ -1487 +1495 @@
         // Calculate residuals
         psImageInit(residual->image, 0.0);
+        psImageInit(difference->image, 0.0);
 
         psKernel *weight = NULL;
@@ -1537 +1546 @@
             }
 
-            // generate the residual image
+            // Generate the difference, residual, and convolved source images.  Note the we
+            // accumulate the convolution of (A-B), so we need to replace it to generate the
+            // images of the convolved source image.
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
+                    difference->kernel[y][x] = target->kernel[y][x] - source->kernel[y][x] * norm - background;
+                    residual->kernel[y][x] = difference->kernel[y][x] - convolved1->kernel[y][x];
                     convolved1->kernel[y][x] += source->kernel[y][x] * norm;
-                    residual->kernel[y][x] = target->kernel[y][x] - convolved1->kernel[y][x] - background;
                 }
             }
+
             // XXX if we want to have a weight and window, we'll need to pass through to here
-            pmSubtractionChisqStats(fluxes, chisq, moments, stampMask, convolved1, NULL, residual, weight, window);
+            pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, NULL, difference, residual, weight, window);
 
         } else {
@@ -1574 +1587 @@
             }
 
+            // Generate the difference, residual, and convolved source images.  Note the we
+            // accumulate the convolutions of (A-B), so we need to replace (A or B) to generate
+            // the images of the convolved source images.
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
+                    difference->kernel[y][x] = image2->kernel[y][x] - image1->kernel[y][x] * norm - background;
+                    residual->kernel[y][x] = difference->kernel[y][x] + convolved2->kernel[y][x] - convolved1->kernel[y][x];
                     convolved1->kernel[y][x] += image1->kernel[y][x] * norm;
                     convolved2->kernel[y][x] += image2->kernel[y][x];
-                    residual->kernel[y][x] = convolved2->kernel[y][x] - convolved1->kernel[y][x] - background;
                 }
             }
@@ -1589 +1606 @@
             } 
 
-            pmSubtractionChisqStats(fluxes, chisq, moments, stampMask, convolved1, convolved2, residual, weight, window);
+            pmSubtractionChisqStats(fluxes, chisqD, chisqR, moments, stampMask, convolved1, convolved2, difference, residual, weight, window);
         }
     }
@@ -1595 +1612 @@
     // find the mean chisq and mean moment
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
-    psVectorStats (stats, chisq, NULL, stampMask, 0xff);
-    float chisqValue = stats->sampleMean;
+    psVectorStats (stats, chisqD, NULL, stampMask, 0xff);
+    float chisqDValue = stats->sampleMean;
+
+    psStatsInit(stats);
+    psVectorStats (stats, chisqR, NULL, stampMask, 0xff);
+    float chisqRValue = stats->sampleMean;
 
     psStatsInit(stats);
@@ -1637 +1658 @@
     // penalized by increasing the score somewhat.  the 0.01 value is not well-chosen.
     float orderFactor = 0.01 * kernels->spatialOrder;
-    float score = 2.0 * chisqValue / (sumKernel1 + sumKernel2) + orderFactor;
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "chisq: %6.3f, moment: %6.3f, sumKernel_1: %6.3f, sumKernel_2, score: %6.3f: %6.3f\n", chisqValue, momentValue, sumKernel1, sumKernel2, score);
+    float score = 2.0 * chisqRValue / (sumKernel1 + sumKernel2) + orderFactor;
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "chisq: %6.3f, chisqD: %6.3f, moment: %6.3f, sumKernel_1: %6.3f, sumKernel_2, score: %6.3f: %6.3f\n", chisqRValue, chisqDValue, momentValue, sumKernel1, sumKernel2, score);
 
     // save this result if it is the first or the best (skip if bestMatch is NULL)
@@ -1663 +1684 @@
             match->nGood        = nGood;
             match->fluxes       = psMemIncrRefCounter(fluxes);
-            match->chisq        = psMemIncrRefCounter(chisq);
+            match->chisq        = psMemIncrRefCounter(chisqR);
             match->moments      = psMemIncrRefCounter(moments);
             match->stampMask    = psMemIncrRefCounter(stampMask);
@@ -1669 +1690 @@
     }
 
-    pmSubtractionVisualPlotChisqAndMoments(fluxes, chisq, moments);
+    pmSubtractionVisualPlotChisqAndMoments(fluxes, chisqR, moments);
 
     psFree(stats);
-    psFree(chisq);
+    psFree(chisqR);
+    psFree(chisqD);
     psFree(fluxes);
     psFree(moments);