Changeset 29148 for branches/eam_branches/ipp-20100823/psModules/src/imcombine/pmSubtractionVisual.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/pmSubtractionVisual.c (modified) (5 diffs)

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

@@ -152 +152 @@
 
 /** Plot the least-squares matrix of each stamp */
-bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps, bool dual) {
+bool pmSubtractionVisualPlotLeastSquares (pmSubtractionStampList *stamps) {
 
     if (!pmVisualTestLevel("ppsub.chisq", 1)) return true;
@@ -209 +209 @@
     pmVisualScaleImage(kapa1, canvas32, "Least_Squares", 0, true);
 
+    if (0) {
+        static int count = 0;
+        char filename[64];
+        sprintf (filename, "chisq.%02d.fits", count);
+        count ++;
+        psFits *fits = psFitsOpen (filename, "w");
+        psFitsWriteImage (fits, NULL, canvas32, 0, NULL);
+        psFitsClose (fits);
+    }
+
     pmVisualAskUser(&plotLeastSquares);
     psFree(canvas);
@@ -299 +309 @@
         if (!isfinite(stamp->flux)) continue;
         if (!stamp->convolutions1 && !stamp->convolutions2) continue;
+        fprintf (stderr, "flux: %f, maxFlux: %f  ", stamp->flux, maxFlux);
         if (!maxStamp) {
             maxFlux = stamp->flux;
             maxStamp = stamp;
+            fprintf (stderr, "maxStamp %d\n", i);
             continue;
+        } else {
+            fprintf (stderr, "\n");
         }
         if (stamp->flux > maxFlux) {
@@ -339 +353 @@
             
             double sum = 0.0;
+            double sum2 = 0.0;
             for (int y = -footprint; y <= footprint; y++) {
                 for (int x = -footprint; x <= footprint; x++) {
                     output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
                     sum += kernel->kernel[y][x];
+                    sum2 += PS_SQR(kernel->kernel[y][x]);
                 }
             }
-            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
+            fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
         }               
         pmVisualScaleImage(kapa2, output, "Image", 0, true);
-    }                                    
+
+        if (0) {
+            psFits *fits = psFitsOpen("basis.1.fits", "w");
+            psFitsWriteImage(fits, NULL, output, 0, NULL);
+            psFitsClose(fits);
+        }
+    }
         
     if (maxStamp->convolutions2) {
@@ -373 +395 @@
             
             double sum = 0.0;
+            double sum2 = 0.0;
             for (int y = -footprint; y <= footprint; y++) {
                 for (int x = -footprint; x <= footprint; x++) {
                     output->data.F32[y + yPix][x + xPix] = kernel->kernel[y][x];
                     sum += kernel->kernel[y][x];
+                    sum2 += PS_SQR(kernel->kernel[y][x]);
                 }
             }
-            fprintf (stderr, "kernel %d, sum %f\n", i, sum);
+            fprintf (stderr, "kernel %d, sum %f, sum2: %e\n", i, sum, sum2);
         }               
         pmVisualScaleImage(kapa2, output, "Image", 1, true);
+
+        if (0) {
+            psFits *fits = psFitsOpen("basis.2.fits", "w");
+            psFitsWriteImage(fits, NULL, output, 0, NULL);
+            psFitsClose(fits);
+        }
     }