Changeset 29148 for branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionEquation.c

Timestamp:

Sep 12, 2010, 12:38:23 PM (16 years ago)

Author:

eugene

Message:

re-scale kernels by fwhm^{2 to avoid excessive dynamic range; modify residual stats to be more sensible; add threshold to min significant element of SVD inversion (1e-6); adjust penalties again (still not clear what is best...)}

File:

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

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

@@ -972 +972 @@
             pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
             if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
+                
                 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
                 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
+
                 psVectorAppend(norms, stamp->norm);
+
                 pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
                 numStamps++;
@@ -1107 +1110 @@
         }
 
-#ifdef TESTING
-        psFitsWriteImageSimple ("sumMatrix.fits", sumMatrix, NULL);
+#if 0
+        psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
+        psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
         psVectorWriteFile("sumVector.dat", sumVector);
+        psFree (save);
 #endif
 
@@ -1191 +1196 @@
             sumVector->data.F64[normIndex] = 0.0;
 
-            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
-
+// save the matrix and vector after the NULLs have been set
+#if 0
+            psImage *save = psImageCopy(NULL, sumMatrix, PS_TYPE_F32);
+            psFitsWriteImageSimple ("sumMatrix.fits", save, NULL);
+            psVectorWriteFile("sumVector.dat", sumVector);
+            psFree (save);
+#endif
+
+            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-6);
+            // solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 3e-4);
+            // psVectorCopy (solution, sumVector, PS_TYPE_F64);
+            // psMatrixGJSolve(sumMatrix, solution);
             solution->data.F64[normIndex] = normValue;
         }
@@ -1265 +1280 @@
 }
 
-bool pmSubtractionResidualStats(psVector *fSigRes, psVector *fMaxRes, psVector *fMinRes, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) {
-
-    // XXX measure some useful stats on the residuals
+// measure some useful stats on the stamp residuals:
+// fResSigma : the residual stdev / total flux
+// fResOuter : the residual fabs / total flux for R > 2 pix
+// fResTotal : the residual fabs / total flux for R > 0 pix
+bool pmSubtractionResidualStats(psVector *fResSigma, psVector *fResOuter, psVector *fResTotal, psKernel *target, psKernel *source, psKernel *residual, double norm, int footprint) {
+
     float sum = 0.0;
     float peak = 0.0;
@@ -1277 +1295 @@
     }
 
-    // only count pixels with more than X% of the source flux
-    // calculate stdev(dflux)
+    // init counters
+    int npix = 0;
     float dflux1 = 0.0;
     float dflux2 = 0.0;
-    int npix = 0;
-
-    float dmax = 0.0;
-    float dmin = 0.0;
-
-    // XXX update these with a bit more rigour
+    float dOuter = 0.0;
+    float dTotal = 0.0;
+
     for (int y = - footprint; y <= footprint; y++) {
         for (int x = - footprint; x <= footprint; x++) {
-          // float dflux = 0.5*(target->kernel[y][x] + source->kernel[y][x] * norm);
-            // if (dflux < 0.02*sum) continue;
             dflux1 += residual->kernel[y][x];
             dflux2 += PS_SQR(residual->kernel[y][x]);
-            // dmax = PS_MAX(residual->kernel[y][x], dmax);
-            // dmin = PS_MIN(residual->kernel[y][x], dmin);
-            dmax += fabs(residual->kernel[y][x]);
-
+            dTotal += fabs(residual->kernel[y][x]);
             if (hypot(x,y) > 2.0) {
-              dmin += fabs(residual->kernel[y][x]);
+              dOuter += fabs(residual->kernel[y][x]);
             }
             npix ++;
@@ -1305 +1315 @@
     float sigma = sqrt(dflux2 / npix - PS_SQR(dflux1/npix));
     if (!isfinite(sum))  return false;
-    if (!isfinite(dmax)) return false;
-    if (!isfinite(dmin)) return false;
     if (!isfinite(peak)) return false;
-
-    fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dmax/sum, dmin/sum);
-    psVectorAppend(fSigRes, sigma/sum);
-    psVectorAppend(fMaxRes, dmax/sum);
-    psVectorAppend(fMinRes, dmin/sum);
+    if (!isfinite(dOuter)) return false;
+    if (!isfinite(dTotal)) return false;
+
+    fprintf (stderr, "sum: %f, peak: %f, sigma: %f, fsigma: %f, fmax: %f, fmin: %f\n", sum, peak, sigma, sigma/sum, dOuter/sum, dTotal/sum);
+    psVectorAppend(fResSigma, sigma/sum);
+    psVectorAppend(fResOuter, dOuter/sum);
+    psVectorAppend(fResTotal, dTotal/sum);
     return true;
 }
@@ -1335 +1345 @@
     pmSubtractionVisualShowFitInit (stamps);
 
-    psVector *fSigRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
-    psVector *fMinRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
-    psVector *fMaxRes = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *fResSigma = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *fResOuter = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
+    psVector *fResTotal = psVectorAllocEmpty(stamps->num, PS_TYPE_F32);
 
     // we want to save the residual images for the 9 brightest stamps.
@@ -1449 +1459 @@
             }
 
-            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, target, source, residual, norm, footprint);
+            pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, target, source, residual, norm, footprint);
 
         } else {
@@ -1485 +1495 @@
             }
 
-            pmSubtractionResidualStats(fSigRes, fMaxRes, fMinRes, image1, image2, residual, norm, footprint);
+            pmSubtractionResidualStats(fResSigma, fResOuter, fResTotal, image1, image2, residual, norm, footprint);
         }
 
@@ -1564 +1574 @@
 
         psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-        psVectorStats (stats, fSigRes, NULL, NULL, 0);
-        kernels->fSigResMean = stats->robustMedian;
-        kernels->fSigResStdev = stats->robustStdev;
+        psVectorStats (stats, fResSigma, NULL, NULL, 0);
+        kernels->fResSigmaMean  = stats->robustMedian;
+        kernels->fResSigmaStdev = stats->robustStdev;
 
         psStatsInit (stats);
-        psVectorStats (stats, fMaxRes, NULL, NULL, 0);
-        kernels->fMaxResMean = stats->robustMedian;
-        kernels->fMaxResStdev = stats->robustStdev;
+        psVectorStats (stats, fResOuter, NULL, NULL, 0);
+        kernels->fResOuterMean  = stats->robustMedian;
+        kernels->fResOuterStdev = stats->robustStdev;
 
         psStatsInit (stats);
-        psVectorStats (stats, fMinRes, NULL, NULL, 0);
-        kernels->fMinResMean = stats->robustMedian;
-        kernels->fMinResStdev = stats->robustStdev;
+        psVectorStats (stats, fResTotal, NULL, NULL, 0);
+        kernels->fResTotalMean  = stats->robustMedian;
+        kernels->fResTotalStdev = stats->robustStdev;
 
         // XXX save these values somewhere
-        psLogMsg("psModules.imcombine", PS_LOG_INFO, "fSigma: %f +/- %f, fMaxRes: %f +/- %f, fMinRes: %f +/- %f",
-                 kernels->fSigResMean, kernels->fSigResStdev,
-                 kernels->fMaxResMean, kernels->fMaxResStdev,
-                 kernels->fMinResMean, kernels->fMinResStdev);
-
-        psFree (fSigRes);
-        psFree (fMaxRes);
-        psFree (fMinRes);
+        psLogMsg("psModules.imcombine", PS_LOG_INFO, "fResSigma: %f +/- %f, fResOuter: %f +/- %f, fResTotal: %f +/- %f",
+                 kernels->fResSigmaMean, kernels->fResSigmaStdev,
+                 kernels->fResOuterMean, kernels->fResOuterStdev,
+                 kernels->fResTotalMean, kernels->fResTotalStdev);
+
+        psFree (fResSigma);
+        psFree (fResOuter);
+        psFree (fResTotal);
         psFree (stats);
     }
@@ -1592 +1602 @@
     psFree(residual);
     psFree(polyValues);
-
 
     return deviations;