Changeset 26505 for branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c

Timestamp:

Jan 4, 2010, 10:06:09 AM (17 years ago)

Author:

eugene

Message:

add skyErr, a way of increasing the constant portion of the variance in the kernel analysis: var = (readNoise^{2 + skyErr + gain*flux + sysErr*flux}2); clean up testing verbosity; temporary dump residuals and quit in visualization code

File:

• branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c (modified) (15 diffs)

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

@@ -181 +181 @@
 
 pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region,
-                                                    int footprint, float spacing, float sysErr)
+                                                    int footprint, float spacing, float sysErr, float skyErr)
 {
     pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return
@@ -224 +224 @@
     list->footprint = footprint;
     list->sysErr = sysErr;
+    list->skyErr = skyErr;
 
     return list;
@@ -326 +327 @@
                                                 const psImage *image2, const psImage *subMask,
                                                 const psRegion *region, float thresh1, float thresh2,
-                                                int size, int footprint, float spacing, float sysErr,
+                                                int size, int footprint, float spacing, float sysErr, float skyErr,
                                                 pmSubtractionMode mode)
 {
@@ -383 +384 @@
 
     if (!stamps) {
-        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr);
+        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing, sysErr, skyErr);
     }
 
@@ -441 +442 @@
                     xStamp = xList->data.F32[j];
                     yStamp = yList->data.F32[j];
-                    fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp);
+                    // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f\n", xCentre, yCentre, xStamp, yStamp);
                 }
             } else {
@@ -496 +497 @@
                                                const psImage *image, const psImage *subMask,
                                                const psRegion *region, int size, int footprint,
-                                               float spacing, float sysErr, pmSubtractionMode mode)
+                                               float spacing, float sysErr, float skyErr, pmSubtractionMode mode)
 
 {
@@ -517 +518 @@
     pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
                                                                  region, footprint, spacing,
-                                                                 sysErr); // Stamp list
+                                                                 sysErr, skyErr); // Stamp list
     int numStamps = stamps->num;        // Number of stamps
 
@@ -548 +549 @@
         }
 
-        fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix);
+        // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix);
 
         bool found = false;
@@ -742 +743 @@
             return false;
         }
-        fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y);
+        // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", stamp->x, stamp->y, x, y);
 
         // Catch memory leaks --- these should have been freed and NULLed before
@@ -765 +766 @@
             psImage *varSub = psImageSubset(variance, region); // Subimage with stamp
             psKernel *var = psKernelAllocFromImage(varSub, size, size); // Variance postage stamp
-            if (isfinite(stamps->sysErr) && stamps->sysErr > 0) {
+
+            if ((isfinite(stamps->skyErr) && (stamps->skyErr > 0)) ||
+                (isfinite(stamps->sysErr) && (stamps->sysErr > 0))) {
                 float sysErr = 0.25 * PS_SQR(stamps->sysErr); // Systematic error
+                float skyErr = stamps->skyErr;
                 psKernel *image1 = stamp->image1, *image2 = stamp->image2; // Input images
                 for (int y = -size; y <= size; y++) {
                     for (int x = -size; x <= size; x++) {
                         float additional = image1->kernel[y][x] + image2->kernel[y][x];
-                        weight->kernel[y][x] = 1.0 / (var->kernel[y][x] + sysErr * PS_SQR(additional));
+                        weight->kernel[y][x] = 1.0 / (skyErr + var->kernel[y][x] + sysErr * PS_SQR(additional));
                     }
                 }
@@ -795 +799 @@
 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image,
                                                           const psImage *subMask, const psRegion *region,
-                                                          int size, int footprint, float spacing, float sysErr,
+                                                          int size, int footprint, float spacing, float sysErr, float skyErr,
                                                           pmSubtractionMode mode)
 {
@@ -819 +823 @@
             y->data.F32[numOut] = source->peak->yf;
         }
-        fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
+        // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
         numOut++;
     }
@@ -826 +830 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size,
-                                                            footprint, spacing, sysErr, mode); // Stamps
+                                                            footprint, spacing, sysErr, skyErr, mode); // Stamps
     psFree(x);
     psFree(y);
@@ -840 +844 @@
 pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image,
                                                        const psImage *subMask, const psRegion *region,
-                                                       int size, int footprint, float spacing, float sysErr,
+                                                       int size, int footprint, float spacing, float sysErr, float skyErr,
                                                        pmSubtractionMode mode)
 {
@@ -858 +862 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint,
-                                                            spacing, sysErr, mode);
+                                                            spacing, sysErr, skyErr, mode);
     psFree(data);