Changeset 15516 for branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionStamps.c

Timestamp:

Nov 8, 2007, 12:59:04 PM (19 years ago)

Author:

eugene

Message:

accepting updates from HEAD

File:

• branches/eam_branch_20071023/psModules/src/imcombine/pmSubtractionStamps.c (modified) (21 diffs)

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

@@ -4 +4 @@
 
 #include <stdio.h>
+#include <string.h>
 #include <pslib.h>
 
@@ -45 +46 @@
     psFree(stamp->matrix);
     psFree(stamp->vector);
-    psFree(stamp->reference);
-    psFree(stamp->input);
+    psFree(stamp->image1);
+    psFree(stamp->image2);
     psFree(stamp->weight);
-    psFree(stamp->convolutions);
+    psFree(stamp->convolutions1);
+    psFree(stamp->convolutions2);
 }
 
@@ -65 +67 @@
 // Is this position unmasked?
 static bool checkStampMask(int x, int y, // Coordinates of stamp
-                           const psImage *mask
+                           const psImage *mask, // Mask
+                           pmSubtractionMode mode // Mode for subtraction
                            )
 {
@@ -74 +77 @@
         return false;
     }
-    return (mask->data.PS_TYPE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER |
-                                                  PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_REJ)) ?
-        false : true;
+
+    psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ; // Mask value
+    switch (mode) {
+      case PM_SUBTRACTION_MODE_1:
+      case PM_SUBTRACTION_MODE_TARGET:
+        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1;
+        break;
+      case PM_SUBTRACTION_MODE_2:
+        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_2;
+        break;
+      case PM_SUBTRACTION_MODE_UNSURE:
+      case PM_SUBTRACTION_MODE_DUAL:
+        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1 | PM_SUBTRACTION_MASK_FOOTPRINT_2;
+        break;
+      default:
+        psAbort("Unsupported subtraction mode: %x", mode);
+    }
+
+    return (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) ? false : true;
 }
 
@@ -140 +159 @@
     stamp->status = PM_SUBTRACTION_STAMP_INIT;
 
-    stamp->reference = NULL;
-    stamp->input = NULL;
+    stamp->image1 = NULL;
+    stamp->image2 = NULL;
     stamp->weight = NULL;
-    stamp->convolutions = NULL;
+    stamp->convolutions1 = NULL;
+    stamp->convolutions2 = NULL;
 
     return stamp;
@@ -151 +171 @@
 pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image,
                                                 const psImage *subMask, const psRegion *region,
-                                                float threshold, float spacing)
+                                                float threshold, float spacing, pmSubtractionMode mode)
 {
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
@@ -226 +246 @@
             for (int y = subRegion->y0; y <= subRegion->y1 ; y++) {
                 for (int x = subRegion->x0; x <= subRegion->y1 ; x++) {
-                    if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxStamp) {
+                    if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) {
                         fluxStamp = image->data.F32[y][x];
                         xStamp = x;
@@ -242 +262 @@
 
             // Reset the postage stamps since we're making a new stamp
-            psFree(stamp->reference);
-            psFree(stamp->input);
+            psFree(stamp->image1);
+            psFree(stamp->image2);
             psFree(stamp->weight);
-            stamp->reference = stamp->input = stamp->weight = NULL;
-
-            stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
+            psFree(stamp->convolutions1);
+            psFree(stamp->convolutions2);
+            stamp->image1 = stamp->image2 = stamp->weight = NULL;
+            stamp->convolutions1 = stamp->convolutions2 = NULL;
+
+            stamp->status = PM_SUBTRACTION_STAMP_FOUND;
             numFound++;
             psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n",
@@ -266 +289 @@
 pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux,
                                                const psImage *image, const psImage *subMask,
-                                               const psRegion *region, float spacing, int exclusionZone)
+                                               const psRegion *region, float spacing, pmSubtractionMode mode)
+
 {
     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
@@ -289 +313 @@
     }
     PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(exclusionZone, NULL);
 
     int numStars = x->n;                // Number of stars
-    psVector *exclude = psVectorAlloc(numStars, PS_TYPE_U8); // Exclude a star?
-    psVectorInit(exclude, 0);
-
-    // Apply exclusion zone around each star --- important when we're convolving to a specified PSF; less so
-    // when we're trying to get a good subtraction.
-    if (exclusionZone > 0) {
-        psTrace("psModules.imcombine", 2, "Applying exclusion zone of %d pixels for stamps\n", exclusionZone);
-        // We use something resembling a binary search --- coordinates are sorted in the x dimension, and then
-        // to exclude stars within a nominated distance, we need only examine (i.e., calculate the
-        // 2-dimensional distance to) other stars in the list that are within that distance of the x
-        // coordinate of the star of interest.  By marking both stars when we find a violation of the
-        // exclusion zone, we need only search upwards from the x coordinate of the star of interest.
-
-        int minDistance2 = PS_SQR(exclusionZone); // Minimum square distance for other stars
-        psVector *indexes = psVectorSortIndex(NULL, x); // Indices for sorting in x
-        for (int i = 0; i < numStars - 1; i++) {
-            int iIndex = indexes->data.S32[i]; // Index for star i
-            if (exclude->data.U8[iIndex]) {
-                continue;
-            }
-            float ix = x->data.F32[iIndex], iy = y->data.F32[iIndex]; // Coordinates for star i
-            float jx, jy;                   // Coordinates for star j
-            for (int j = i + 1, jIndex = indexes->data.S32[j];
-                 j < numStars && (jx = x->data.F32[jIndex]) < ix + exclusionZone;
-                 j++, jIndex = indexes->data.S32[j]) {
-                jy = y->data.F32[jIndex];
-                if (PS_SQR(ix - jx) + PS_SQR(iy - jy) < minDistance2) {
-                    exclude->data.U8[iIndex] = 0xff;
-                    exclude->data.U8[jIndex] = 0xff;
-                    psTrace("psModules.imcombine", 5, "Excluding stamps %d,%d and %d,%d\n",
-                            (int)x->data.F32[iIndex], (int)y->data.F32[iIndex],
-                            (int)x->data.F32[jIndex], (int)y->data.F32[jIndex]);
-                    // Would 'break' here, but there might be more stamps within the exclusion zone.
-                }
-            }
-        }
-        psFree(indexes);
-    }
-
     pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
                                                                  region, spacing); // Stamp list
@@ -347 +331 @@
     // Put the stars into their appropriate subregions
     for (int i = 0; i < numStars; i++) {
-        if (exclude->data.U8[i]) {
-            // Star has been excluded
-            continue;
-        }
         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
@@ -357 +337 @@
             continue;
         }
-        if (!checkStampMask(xPix, yPix, subMask)) {
+        if (!checkStampMask(xPix, yPix, subMask, mode)) {
             // Not a good stamp
             continue;
@@ -386 +366 @@
         }
     }
-    psFree(exclude);
 
     // Sort the list by flux, with the brightest last
@@ -421 +400 @@
 
 
-bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *reference, psImage *input,
+bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
                                 psImage *weight, int footprint, int kernelSize)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
-    PS_ASSERT_IMAGE_NON_NULL(reference, false);
-    PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false);
-    if (input) {
-        PS_ASSERT_IMAGE_NON_NULL(input, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false);
-        PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false);
-    }
-    if (weight) {
-        PS_ASSERT_IMAGE_NON_NULL(weight, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(weight, reference, false);
-        PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
-    }
+    PS_ASSERT_IMAGE_NON_NULL(image1, false);
+    PS_ASSERT_IMAGE_TYPE(image1, PS_TYPE_F32, false);
+    if (image2) {
+        PS_ASSERT_IMAGE_NON_NULL(image2, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(image2, image1, false);
+        PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false);
+    }
+    PS_ASSERT_IMAGE_NON_NULL(weight, false);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image1, false);
+    PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
     PS_ASSERT_INT_NONNEGATIVE(footprint, false);
     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
 
-    int numCols = reference->numCols, numRows = reference->numRows; // Size of images
+    int numCols = image1->numCols, numRows = image1->numRows; // Size of images
     int size = kernelSize + footprint; // Size of postage stamps
-
-    if (!weight) {
-        // Use the input (or if I must, the reference) as a rough approximation to the variance map, and HOPE
-        // that it's positive.
-        weight = input ? input : reference;
-    }
 
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
+        if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_FOUND) {
             continue;
         }
@@ -468 +439 @@
         }
 
+        // Catch memory leaks --- these should have been freed and NULLed before
+        assert(stamp->image1 == NULL);
+        assert(stamp->image2 == NULL);
+        assert(stamp->weight == NULL);
+
         psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest
 
-        psImage *refSub = psImageSubset(reference, region); // Subimage with stamp
-        stamp->reference = psKernelAllocFromImage(refSub, size, size);
-        psFree(refSub);                 // Drop reference
-
-        if (input) {
-            psImage *inSub = psImageSubset(input, region); // Subimage with stamp
-            stamp->input = psKernelAllocFromImage(inSub, size, size);
-            psFree(inSub);              // Drop reference
+        psImage *sub1 = psImageSubset(image1, region); // Subimage with stamp
+        stamp->image1 = psKernelAllocFromImage(sub1, size, size);
+        psFree(sub1);                   // Drop reference
+
+        if (image2) {
+            psImage *sub2 = psImageSubset(image2, region); // Subimage with stamp
+            stamp->image2 = psKernelAllocFromImage(sub2, size, size);
+            psFree(sub2);               // Drop reference
         }
 
         psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp
         stamp->weight = psKernelAllocFromImage(wtSub, size, size);
-        psFree(wtSub);                 // Drop reference
+        psFree(wtSub);                  // Drop reference
+
+        stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
     }
 
@@ -488 +466 @@
 }
 
-
+#if 0
 bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize)
 {
@@ -517 +495 @@
         float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame
 
-        psFree(stamp->input);
-        stamp->input = psKernelAlloc(-size, size, -size, size);
-        psKernel *input = stamp->input; // Target stamp
+        psFree(stamp->image2);
+        stamp->image2 = psKernelAlloc(-size, size, -size, size);
+        psKernel *target = stamp->image2; // Target stamp
 
         // Put in a Waussian, just for fun!
@@ -525 +503 @@
             for (int u = -size; u <= size; u++) {
                 float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma));
-                input->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z));
+                target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z));
             }
         }
@@ -533 +511 @@
     return true;
 }
-
+#endif
 
 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask,
                                                           const psRegion *region, float spacing,
-                                                          int exclusionZone)
+                                                          pmSubtractionMode mode)
 {
     PS_ASSERT_ARRAY_NON_NULL(sources, NULL);
@@ -561 +539 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region,
-                                                            spacing, exclusionZone); // Stamps to return
+                                                            spacing, mode); // Stamps to return
     psFree(x);
     psFree(y);
@@ -572 +550 @@
     return stamps;
 }
+
+
+pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *subMask,
+                                                       const psRegion *region, float spacing,
+                                                       pmSubtractionMode mode)
+{
+    PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
+    // Let pmSubtractionStampsSet take care of the rest of the assertions
+
+    const char *format = (mode == PM_SUBTRACTION_MODE_TARGET ? "%f %f %f" : "%f %f"); // Format of file
+    psArray *data = psVectorsReadFromFile(filename, format);
+    if (!data) {
+        psError(PS_ERR_IO, false, "Unable to read stamps file %s", filename);
+        return NULL;
+    }
+    psVector *x = data->data[0], *y = data->data[1]; // Stamp positions
+    psVector *flux = (mode == PM_SUBTRACTION_MODE_TARGET ? data->data[2] : NULL); // Stamp fluxes
+
+    // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
+    psBinaryOp(x, x, "-", psScalarAlloc(1.0, PS_TYPE_F32));
+    psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32));
+
+    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, spacing,
+                                                            mode);
+    psFree(data);
+
+    return stamps;
+
+}