Changeset 27838 for branches/tap_branches/psModules/src/imcombine/pmSubtractionStamps.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/imcombine/pmSubtractionStamps.c (modified) (38 diffs)

branches/tap_branches/psModules/src/imcombine/pmSubtractionStamps.c

@@ -46 +46 @@
     psFree(list->y);
     psFree(list->flux);
+    psFree(list->window);
 }
 
@@ -54 +55 @@
     psFree(stamp->image1);
     psFree(stamp->image2);
-    psFree(stamp->variance);
+    psFree(stamp->weight);
     psFree(stamp->convolutions1);
     psFree(stamp->convolutions2);
-
-    psFree(stamp->matrix1);
-    psFree(stamp->matrix2);
-    psFree(stamp->matrixX);
-    psFree(stamp->vector1);
-    psFree(stamp->vector2);
-
+    psFree(stamp->matrix);
+    psFree(stamp->vector);
 }
 
@@ -80 +76 @@
 
 // Search a region for a suitable stamp
-bool stampSearch(int *xStamp, int *yStamp, // Coordinates of stamp, to return
+bool stampSearch(float *xStamp, float *yStamp, // Coordinates of stamp, to return
                  float *fluxStamp, // Flux of stamp, to return
                  const psImage *image1, const psImage *image2, // Images to search
@@ -93 +89 @@
     *fluxStamp = -INFINITY;             // Flux of best stamp
 
+    // fprintf (stderr, "xMin, xMax: %d, %d -> ", xMin, xMax);
+
     // Ensure we're not going to go outside the bounds of the image
     xMin = PS_MAX(border, xMin);
@@ -99 +97 @@
     yMax = PS_MIN(numRows - border - 1, yMax);
 
+    if (xMax < xMin) return false;
+    if (yMax < yMin) return false;
+
+    psAssert (xMin <= xMax, "x mismatch?");
+    psAssert (yMin <= yMax, "y mismatch?");
+
     for (int y = yMin; y <= yMax; y++) {
         for (int x = xMin; x <= xMax; x++) {
@@ -115 +119 @@
                 ((image1) ? image1->data.F32[y][x] : image2->data.F32[y][x]); // Flux at pixel
             if (flux > *fluxStamp) {
-                *xStamp = x;
-                *yStamp = y;
+                *xStamp = x + 0.5;
+                *yStamp = y + 0.5;
                 *fluxStamp = flux;
                 found = true;
@@ -180 +184 @@
 
 pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region,
-                                                    int footprint, float spacing)
+                                                    int footprint, float spacing, float normFrac,
+                                                    float sysErr, float skyErr)
 {
     pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return
@@ -220 +225 @@
     list->y = NULL;
     list->flux = NULL;
+    list->window = NULL;
+    list->normFrac = normFrac;
+    list->normWindow = 0;
     list->footprint = footprint;
+    list->sysErr = sysErr;
+    list->skyErr = skyErr;
 
     return list;
@@ -238 +248 @@
     out->y = NULL;
     out->flux = NULL;
+    out->window = psMemIncrRefCounter(in->window);
     out->footprint = in->footprint;
+    out->normWindow = in->normWindow;
 
     for (int i = 0; i < num; i++) {
@@ -256 +268 @@
         outStamp->image1 = inStamp->image1 ? psKernelCopy(inStamp->image1) : NULL;
         outStamp->image2 = inStamp->image2 ? psKernelCopy(inStamp->image2) : NULL;
-        outStamp->variance = inStamp->variance ? psKernelCopy(inStamp->variance) : NULL;
+        outStamp->weight = inStamp->weight ? psKernelCopy(inStamp->weight) : NULL;
 
         if (inStamp->convolutions1) {
@@ -279 +291 @@
         }
 
-        outStamp->matrix1 = inStamp->matrix1 ? psImageCopy(NULL, inStamp->matrix1, PS_TYPE_F64) : NULL;
-        outStamp->matrix2 = inStamp->matrix2 ? psImageCopy(NULL, inStamp->matrix2, PS_TYPE_F64) : NULL;
-        outStamp->matrixX = inStamp->matrixX ? psImageCopy(NULL, inStamp->matrixX, PS_TYPE_F64) : NULL;
-        outStamp->vector1 = inStamp->vector1 ? psVectorCopy(NULL, inStamp->vector1, PS_TYPE_F64) : NULL;
-        outStamp->vector2 = inStamp->vector2 ? psVectorCopy(NULL, inStamp->vector2, PS_TYPE_F64) : NULL;
+        outStamp->matrix = inStamp->matrix ? psImageCopy(NULL, inStamp->matrix, PS_TYPE_F64) : NULL;
+        outStamp->vector = inStamp->vector ? psVectorCopy(NULL, inStamp->vector, PS_TYPE_F64) : NULL;
 
         out->stamps->data[i] = outStamp;
@@ -305 +314 @@
     stamp->image1 = NULL;
     stamp->image2 = NULL;
-    stamp->variance = NULL;
+    stamp->weight = NULL;
     stamp->convolutions1 = NULL;
     stamp->convolutions2 = NULL;
 
-    stamp->matrix1 = NULL;
-    stamp->matrix2 = NULL;
-    stamp->matrixX = NULL;
-    stamp->vector1 = NULL;
-    stamp->vector2 = NULL;
+    stamp->matrix = NULL;
+    stamp->vector = NULL;
+    stamp->norm = NAN;
 
     return stamp;
@@ -322 +329 @@
                                                 const psImage *image2, const psImage *subMask,
                                                 const psRegion *region, float thresh1, float thresh2,
-                                                int size, int footprint, float spacing,
-                                                pmSubtractionMode mode)
+                                                int size, int footprint, float spacing, float normFrac,
+                                                float sysErr, float skyErr, pmSubtractionMode mode)
 {
     if (!image1 && !image2) {
@@ -362 +369 @@
         if (psRegionIsNaN(*region)) {
             psString string = psRegionToString(*region);
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string);
+            psError(PM_ERR_PROG, true, "Input region (%s) contains NAN values", string);
             psFree(string);
             return false;
@@ -369 +376 @@
             region->y0 < 0 || region->y1 > numRows) {
             psString string = psRegionToString(*region);
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)",
+            psError(PM_ERR_PROG, true, "Input region (%s) does not fit in image (%dx%d)",
                     string, numCols, numRows);
             psFree(string);
@@ -379 +386 @@
 
     if (!stamps) {
-        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing);
+        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing,
+                                             normFrac, sysErr, skyErr);
     }
 
@@ -401 +409 @@
             numSearch++;
 
-            int xStamp = 0, yStamp = 0; // Coordinates of stamp
+            float xStamp = 0.0, yStamp = 0.0; // Coordinates of stamp
             float fluxStamp = -INFINITY; // Flux of stamp
             bool goodStamp = false;     // Found a good stamp?
@@ -413 +421 @@
                 // Take stamps off the top of the (sorted) list
                 for (int j = xList->n - 1; j >= 0 && !goodStamp; j--) {
-                    int xCentre = xList->data.F32[j] + 0.5, yCentre = yList->data.F32[j] + 0.5;// Stamp centre
-
                     // Chop off the top of the list
                     xList->n = j;
@@ -420 +426 @@
                     fluxList->n = j;
 
+#if 0
                     // Fish around a bit to see if we can find a pixel that isn't masked
+                    // This is not a good idea if we're using the window feature
                     psTrace("psModules.imcombine", 7, "Searching for stamp %d around %d,%d\n",
                             i, xCentre, yCentre);
 
                     // Search bounds
+                    int xCentre = xList->data.F32[j] - 0.5, yCentre = yList->data.F32[j] - 0.5;// Stamp centre
                     int search = footprint - size; // Search radius
                     int xMin = PS_MAX(border, xCentre - search);
@@ -433 +442 @@
                     goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
                                             subMask, xMin, xMax, yMin, yMax, numCols, numRows, border);
+                    // fprintf (stderr, "find: %d %d ==> %5.1f %5.1f (\n", xCentre, yCentre, xStamp, yStamp);
+#else
+                    // Only search the exact centre pixel
+                    goodStamp = stampSearch(&xStamp, &yStamp, &fluxStamp, image1, image2, thresh1, thresh2,
+                                            subMask, xList->data.F32[j], xList->data.F32[j],
+                                            yList->data.F32[j], yList->data.F32[j], numCols, numRows, border);
+#endif
+                    if (0 && goodStamp) {
+                        fprintf (stderr, "find: %6.1f < %6.1f < %6.1f, %6.1f < %6.1f %6.1f\n",
+                                 region->x0 + size, xStamp, region->x1 - size,
+                                 region->y0 + size, yStamp, region->y1 - size);
+                    }
                 }
             } else {
@@ -452 +473 @@
                 psFree(stamp->image1);
                 psFree(stamp->image2);
-                psFree(stamp->variance);
+                psFree(stamp->weight);
                 psFree(stamp->convolutions1);
                 psFree(stamp->convolutions2);
-                stamp->image1 = stamp->image2 = stamp->variance = NULL;
+                stamp->image1 = stamp->image2 = stamp->weight = NULL;
                 stamp->convolutions1 = stamp->convolutions2 = NULL;
 
@@ -487 +508 @@
                                                const psImage *image, const psImage *subMask,
                                                const psRegion *region, int size, int footprint,
-                                               float spacing, pmSubtractionMode mode)
+                                               float spacing, float normFrac, float sysErr, float skyErr,
+                                               pmSubtractionMode mode)
 
 {
@@ -507 +529 @@
     int numStars = x->n;                // Number of stars
     pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
-                                                                 region, footprint, spacing); // Stamp list
+                                                                 region, footprint, spacing,
+                                                                 normFrac, sysErr, skyErr); // Stamp list
     int numStamps = stamps->num;        // Number of stamps
 
@@ -514 +537 @@
     psStringAppend(&ds9name, "stamps_all_%d.ds9", ds9num);
     FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "all stamps"); // ds9 region file
+    psFree(ds9name);
     ds9num++;
 
@@ -529 +553 @@
     for (int i = 0; i < numStars; i++) {
         float xStamp = x->data.F32[i], yStamp = y->data.F32[i]; // Coordinates of stamp
-        int xPix = xStamp + 0.5, yPix = yStamp + 0.5; // Pixel coordinate of stamp
+        int xPix = xStamp - 0.5, yPix = yStamp - 0.5; // Pixel coordinate of stamp
         if (!checkStampRegion(xPix, yPix, region)) {
             // It's not in the big region
@@ -537 +561 @@
             continue;
         }
+
+        // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d\n", xStamp, yStamp, xPix, yPix);
 
         bool found = false;
@@ -605 +631 @@
 
 
+bool pmSubtractionStampsGetWindow(pmSubtractionStampList *stamps, int kernelSize)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
+    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
+
+    int size = stamps->footprint; // Size of postage stamps
+
+    psFree (stamps->window);
+    stamps->window = psKernelAlloc(-size, size, -size, size);
+    psImageInit(stamps->window->image, 0.0);
+
+    // generate normalizations for each stamp
+    psVector *norm1 = psVectorAlloc(stamps->num, PS_TYPE_F32);
+    psVector *norm2 = psVectorAlloc(stamps->num, PS_TYPE_F32);
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (!stamp) continue;
+        if (!stamp->image1) continue;
+        if (!stamp->image2) continue;
+
+        float sum1 = 0.0;
+        float sum2 = 0.0;
+        for (int y = -size; y <= size; y++) {
+            for (int x = -size; x <= size; x++) {
+                sum1 += stamp->image1->kernel[y][x];
+                sum2 += stamp->image2->kernel[y][x];
+            }
+        }
+        norm1->data.F32[i] = sum1;
+        norm2->data.F32[i] = sum2;
+    }
+
+    // storage vector for flux data
+    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psVector *flux1 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
+    psVector *flux2 = psVectorAllocEmpty(2*stamps->num, PS_TYPE_F32);
+
+    // generate the window pixels
+    double sum = 0.0;                   // Sum inside the window
+    float maxValue = 0.0;               // Maximum value, for normalisation
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+
+            flux1->n = 0;
+            flux2->n = 0;
+            for (int i = 0; i < stamps->num; i++) {
+                pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+                if (!stamp) continue;
+                if (!stamp->image1) continue;
+                if (!stamp->image2) continue;
+
+                psVectorAppend(flux1, stamp->image1->kernel[y][x] / norm1->data.F32[i]);
+                psVectorAppend(flux2, stamp->image2->kernel[y][x] / norm2->data.F32[i]);
+            }
+
+            psStatsInit (stats);
+            if (!psVectorStats (stats, flux1, NULL, NULL, 0)) {
+                psAbort ("failed to generate stats");
+            }
+            float f1 = stats->robustMedian;
+            psStatsInit (stats);
+            if (!psVectorStats (stats, flux2, NULL, NULL, 0)) {
+                psAbort ("failed to generate stats");
+            }
+            float f2 = stats->robustMedian;
+
+            stamps->window->kernel[y][x] = f1 + f2;
+            if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(size)) {
+                sum += stamps->window->kernel[y][x];
+            }
+            maxValue = PS_MAX(maxValue, stamps->window->kernel[y][x]);
+        }
+    }
+
+    psTrace("psModules.imcombine", 3, "Window total: %f, threshold: %f\n",
+            sum, (1.0 - stamps->normFrac) * sum);
+    bool done = false;
+    for (int radius = 1; radius <= size && !done; radius++) {
+        double within = 0.0;
+        for (int y = -radius; y <= radius; y++) {
+            for (int x = -radius; x <= radius; x++) {
+                if (PS_SQR(x) + PS_SQR(y) <= PS_SQR(radius)) {
+                    within += stamps->window->kernel[y][x];
+                }
+            }
+        }
+        psTrace("psModules.imcombine", 5, "Radius %d: %f\n", radius, within);
+        if (within > (1.0 - stamps->normFrac) * sum) {
+            stamps->normWindow = radius;
+            done = true;
+        }
+    }
+
+    psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow);
+    if (stamps->normWindow == 0 || stamps->normWindow >= size) {
+        psError(PM_ERR_STAMPS, true, "Unable to determine normalisation window size.");
+        return false;
+    }
+
+    // re-normalize so chisquare values are sensible
+    for (int y = -size; y <= size; y++) {
+        for (int x = -size; x <= size; x++) {
+            stamps->window->kernel[y][x] /= maxValue;
+        }
+    }
+
+#if 0
+    {
+        psFits *fits = psFitsOpen ("window.fits", "w");
+        psFitsWriteImage (fits, NULL, stamps->window->image, 0, NULL);
+        psFitsClose (fits);
+    }
+#endif
+
+    psFree (stats);
+    psFree (flux1);
+    psFree(flux2);
+    psFree (norm1);
+    psFree (norm2);
+    return true;
+}
+
 bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
-                                psImage *variance, int kernelSize)
+                                psImage *variance, int kernelSize, psRegion bounds)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -616 +764 @@
         PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false);
     }
-    PS_ASSERT_IMAGE_NON_NULL(variance, false);
-    PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false);
-    PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false);
-    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
-
-    int numCols = image1->numCols, numRows = image1->numRows; // Size of images
+    if (variance) {
+        PS_ASSERT_IMAGE_NON_NULL(variance, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(variance, image1, false);
+        PS_ASSERT_IMAGE_TYPE(variance, PS_TYPE_F32, false);
+        PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
+    }
+
     int size = kernelSize + stamps->footprint; // Size of postage stamps
 
@@ -630 +779 @@
         }
 
-        if (isnan(stamp->xNorm)) {
-            stamp->xNorm = 2.0 * (stamp->x - (float)numCols/2.0) / (float)numCols;
-        }
-        if (isnan(stamp->yNorm)) {
-            stamp->yNorm = 2.0 * (stamp->y - (float)numRows/2.0) / (float)numRows;
-        }
-
-        int x = stamp->x + 0.5, y = stamp->y + 0.5; // Stamp coordinates
-        if (x < size || x > numCols - size || y < size || y > numRows - size) {
-            psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the image border.\n", i, x, y);
+        p_pmSubtractionPolynomialNormCoords(&stamp->xNorm, &stamp->yNorm, stamp->x, stamp->y,
+                                            bounds.x0, bounds.x1, bounds.y0, bounds.y1);
+
+        int x = stamp->x - 0.5, y = stamp->y - 0.5; // Stamp coordinates
+        // fprintf (stderr, "stamp: %5.1f %5.1f == %d %d (size: %d)\n", stamp->x, stamp->y, x, y, size);
+
+        if (x < bounds.x0 + size || x > bounds.x1 - size || y < bounds.y0 + size || y > bounds.y1 - size) {
+            psError(PM_ERR_PROG, false, "Stamp %d (%d,%d) is within the region border.\n", i, x, y);
             return false;
         }
@@ -646 +793 @@
         assert(stamp->image1 == NULL);
         assert(stamp->image2 == NULL);
-        assert(stamp->variance == NULL);
+        assert(stamp->weight == NULL);
 
         psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest
@@ -660 +807 @@
         }
 
-        psImage *wtSub = psImageSubset(variance, region); // Subimage with stamp
-        stamp->variance = psKernelAllocFromImage(wtSub, size, size);
-        psFree(wtSub);                  // Drop reference
+        psKernel *weight = stamp->weight = psKernelAlloc(-size, size, -size, size); // Weight = 1/variance
+        if (variance) {
+            psImage *varSub = psImageSubset(variance, region); // Subimage with stamp
+            psKernel *var = psKernelAllocFromImage(varSub, size, size); // Variance postage stamp
+
+            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 / (skyErr + var->kernel[y][x] + sysErr * PS_SQR(additional));
+                    }
+                }
+            } else {
+                for (int y = -size; y <= size; y++) {
+                    for (int x = -size; x <= size; x++) {
+                        weight->kernel[y][x] = 1.0 / var->kernel[y][x];
+                    }
+                }
+            }
+            psFree(var);
+            psFree(varSub);
+        } else {
+            psImageInit(weight->image, 1.0);
+        }
 
         stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
@@ -673 +845 @@
                                                           const psImage *subMask, const psRegion *region,
                                                           int size, int footprint, float spacing,
+                                                          float normFrac, float sysErr, float skyErr,
                                                           pmSubtractionMode mode)
 {
@@ -696 +869 @@
             y->data.F32[numOut] = source->peak->yf;
         }
+        // fprintf (stderr, "stamp: %5.1f %5.1f\n", x->data.F32[numOut], y->data.F32[numOut]);
         numOut++;
     }
@@ -702 +876 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size,
-                                                            footprint, spacing, mode); // Stamps to return
+                                                            footprint, spacing, normFrac,
+                                                            sysErr, skyErr, mode); // Stamps
     psFree(x);
     psFree(y);
 
     if (!stamps) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to set stamps from sources.");
+        psError(psErrorCodeLast(), false, "Unable to set stamps from sources.");
     }
 
@@ -716 +891 @@
 pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image,
                                                        const psImage *subMask, const psRegion *region,
-                                                       int size, int footprint, float spacing,
-                                                       pmSubtractionMode mode)
+                                                       int size, int footprint, float spacing, float normFrac,
+                                                       float sysErr, float skyErr, pmSubtractionMode mode)
 {
     PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
@@ -724 +899 @@
     psArray *data = psVectorsReadFromFile(filename, "%f %f");
     if (!data) {
-        psError(PS_ERR_IO, false, "Unable to read stamps file %s", filename);
+        psError(psErrorCodeLast(), false, "Unable to read stamps file %s", filename);
         return NULL;
     }
@@ -734 +909 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, image, subMask, region, size, footprint,
-                                                            spacing, mode);
+                                                            spacing, normFrac, sysErr, skyErr, mode);
     psFree(data);
 
@@ -740 +915 @@
 
 }
+
+
+bool pmSubtractionStampsResetStatus (pmSubtractionStampList *stamps) {
+
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (!stamp) continue;
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) continue;
+        stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
+    }
+    return true;
+}
+