Changeset 30289 for branches/eam_branches/ipp-20101205/psModules/src/imcombine/pmSubtractionStamps.c

Timestamp:

Jan 17, 2011, 5:18:09 PM (16 years ago)

Author:

eugene

Message:

apply window only to the moments portions of the calculation; apply weights only to the chisq portion of the calculation; do not renormalize the chisq portion by the flux (keeps SINGLE and DUAL on same scale); move stamp convolution out of CalculateEquation; explicit tests if the convolutions are needed; remove the old unused matrix solution functions; calculate the coefficient errors; add function to calculate chisq and moments statistics, and to generate a subtraction quality statistic; pull model visualization out of CalculateDeviations; adjust the size of the normalization window (if too small, we likely miss some of the flux); exclude excessively faint sources from model fit

File:

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

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

@@ -122 +122 @@
             if ((image1 && image1->data.F32[y][x] < thresh1) ||
                 (image2 && image2->data.F32[y][x] < thresh2)) {
+                // fprintf (stderr, "%f,%f : thresh\n", xRaw, yRaw);
                 continue;
             }
@@ -429 +430 @@
         stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing,
                                              normFrac, sysErr, skyErr);
+    }
+
+    // XXX TEST : dump all stars in the stamps here
+    if (0) {
+        FILE *f = fopen ("stamp.dat", "w");
+        for (int i = 0; i < stamps->num; i++) {
+            psVector *xList = stamps->x->data[i];
+            psVector *yList = stamps->y->data[i]; // Coordinate lists
+            psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes
+
+            for (int j = 0; j < xList->n; j++) {
+                fprintf (f, "%d %d  %f %f  %f\n", i, j, xList->data.F32[j], yList->data.F32[j], fluxList->data.F32[j]);
+            }
+        }
+        fclose (f);
     }
 
@@ -636 +652 @@
     }
 
+    int nTotal = 0;
+
     // Sort the list by flux, with the brightest last
     for (int i = 0; i < numStamps; i++) {
@@ -662 +680 @@
         stamps->y->data[i] = ySorted;
         stamps->flux->data[i] = fluxSorted;
-    }
-
+        nTotal += num;
+    }
+    // fprintf (stderr, "nTotal %d\n", nTotal);
+    
     return stamps;
 }
@@ -809 +829 @@
     float R2 = Sr2 / Sf2;
 
-    stamps->normWindow1 = 2.0*R1;
-    stamps->normWindow2 = 2.0*R2;
+    stamps->normWindow1 = 2.75*R1;
+    stamps->normWindow2 = 2.75*R2;
     psLogMsg ("psModules.imcombine", PS_LOG_DETAIL, "Kron Radii: %f for 1, %f for 2\n", stamps->normWindow1, stamps->normWindow2);
 
     // Generate a weighting window based on the kron radii
     { 
-        float radius = 1.5 * PS_MAX(R1, R2);
+        float radius = 2.0 * PS_MAX(R1, R2);
 
         psImageInit(stamps->window->image, 0.0);
@@ -826 +846 @@
         }
     }
+
+    // complain if normWindow1 or normWindow2 are larger than size
+    psAssert(stamps->normWindow1 < size, "norm Window 1 too large");
+    psAssert(stamps->normWindow2 < size, "norm Window 2 too large");
 
 # else
@@ -1142 +1166 @@
             continue;
         }
+       
+        // XXX TEST
+        if (source->errMag > 0.1) continue;
+
         if (source->modelPSF) {
             x->data.F32[numOut] = source->modelPSF->params->data.F32[PM_PAR_XPOS];