Changeset 23839

Timestamp:

Apr 13, 2009, 4:05:58 PM (17 years ago)

Author:

Paul Price

Message:

Generate ds9 region files for stamps at interesting places.

Location:

trunk/psModules/src/imcombine

Files:

• pmSubtraction.c (modified) (3 diffs)
• pmSubtractionEquation.c (modified) (4 diffs)
• pmSubtractionStamps.c (modified) (9 diffs)
• pmSubtractionStamps.h (modified) (1 diff)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -850 +850 @@
     psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit);
 
+
+    psString ds9name = NULL;            // Filename for ds9 region file
+    static int ds9num = 0;              // File number for ds9 region file
+    psStringAppend(&ds9name, "stamps_reject_%d.ds9", ds9num);
+    FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "rejected stamps");
+    psFree(ds9name);
+    ds9num++;
+
     int numRejected = 0;                // Number of stamps rejected
     int numGood = 0;                    // Number of good stamps
@@ -872 +880 @@
                     }
                 }
+                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
 
                 // Set stamp for replacement
@@ -897 +906 @@
                 numGood++;
                 newMean += deviations->data.F32[i];
+                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
             }
         }
     }
     newMean /= numGood;
+
+    if (ds9) {
+        fclose(ds9);
+    }
 
     if (numRejected > 0) {

trunk/psModules/src/imcombine/pmSubtractionEquation.c

@@ -749 +749 @@
     }
 
+    psString ds9name = NULL;            // Filename for ds9 region file
+    static int ds9num = 0;              // File number for ds9 region file
+    psStringAppend(&ds9name, "stamps_solution_%d.ds9", ds9num);
+    FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "solution stamps");
+    psFree(ds9name);
+    ds9num++;
+
     if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) {
         // Accumulate the least-squares matricies and vectors
@@ -761 +768 @@
                 (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1);
                 (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1);
+                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
                 numStamps++;
+            } else if (stamp->status == PM_SUBTRACTION_STAMP_REJECTED) {
+                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "red");
             }
         }
@@ -830 +840 @@
                 (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1);
                 (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2);
+                pmSubtractionStampPrint(ds9, stamp->x, stamp->y, stamps->footprint, "green");
                 numStamps++;
             }
@@ -990 +1001 @@
         // XXXXX Free temporary matrices and vectors
 
+    }
+
+    if (ds9) {
+        fclose(ds9);
     }
 

trunk/psModules/src/imcombine/pmSubtractionStamps.c

@@ -149 +149 @@
     return clean;
 }
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Functions for generating ds9 region files
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+static bool ds9regions = false;         // Save ds9 region files?
+
+void pmSubtractionRegions(bool state)
+{
+    ds9regions = state;
+}
+
+FILE *pmSubtractionStampsFile(const pmSubtractionStampList *stamps, const char *filename,
+                              const char *description)
+{
+    if (!ds9regions || !stamps || !filename) {
+        return NULL;
+    }
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Writing %s to ds9 region file: %s",
+             description, filename);
+
+    FILE *file = fopen(filename, "w");
+
+    // Outline the stamps
+    for (int i = 0; i < stamps->num; i++) {
+        psRegion *region = stamps->regions->data[i]; // Region of interest
+        float xCentre = 0.5 * (region->x0 + region->x1), yCentre = 0.5 * (region->y0 + region->y1);
+        fprintf(file, "image;box(%f,%f,%f,%f,0) # color=blue\nimage;text(%f,%f,{%d}) # color=blue\n",
+                xCentre, yCentre, region->x1 - region->x0, region->y1 - region->y0,
+                xCentre, yCentre, i);
+    }
+
+    return file;
+}
+
+void pmSubtractionStampPrint(FILE *ds9, float x, float y, float size, const char *color)
+{
+    if (!ds9regions || !ds9) {
+        return;
+    }
+    fprintf(ds9, "image;circle(%f,%f,%f)", x, y, size);
+    if (color && strlen(color) > 0) {
+        fprintf(ds9, " # color=%s", color);
+    }
+    fprintf(ds9, "\n");
+    return;
+}
+
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -366 +415 @@
 }
 
-bool pmSubtractionStampsRegions(pmSubtractionStampList *stamps, const char *filename)
-{
-    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
-    PS_ASSERT_STRING_NON_EMPTY(filename, false);
-
-    FILE *outFile = fopen(filename, "w"); // File to write
-    if (!outFile) {
-        psError(PS_ERR_IO, false, "Unable to open %s for writing", filename);
-        return false;
-    }
-
-    for (int i = 0; i < stamps->num; i++) {
-        // Outline the stamp
-        psRegion *region = stamps->regions->data[i]; // Region of interest
-        fprintf(outFile, "image;box(%f,%f,%f,%f,0)\n",
-                0.5 * (region->x0 + region->x1), 0.5 * (region->y0 + region->y1),
-                region->x1 - region->x0, region->y1 - region->y0);
-
-        psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Coordinate lists
-
-        for (int j = 0; j < xList->n; j++) {
-            fprintf(outFile, "image;circle(%f,%f,%d)\n",
-                    xList->data.F32[j], yList->data.F32[j], stamps->footprint);
-        }
-    }
-
-    fclose(outFile);
-
-    return true;
-}
+
 
 pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y,
@@ -424 +444 @@
     int numStamps = stamps->num;        // Number of stamps
 
+    psString ds9name = NULL;            // Filename for ds9 region file
+    static int ds9num = 0;              // File number for ds9 region file
+    psStringAppend(&ds9name, "stamps_all_%d.ds9", ds9num);
+    FILE *ds9 = pmSubtractionStampsFile(stamps, ds9name, "all stamps"); // ds9 region file
+    ds9num++;
+
     // Initialise the lists
     stamps->x = psArrayAlloc(numStamps);
@@ -442 +468 @@
             psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because outside region",
                     xPix, yPix);
+            pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "magenta");
             continue;
         }
@@ -448 +475 @@
             psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because bad mask",
                     xPix, yPix);
+            pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "red");
             continue;
         }
@@ -471 +499 @@
                 psTrace("psModules.imcombine", 9, "Putting input stamp (%d,%d) into subregion %d",
                         xPix, yPix, j);
+                pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "green");
             }
         }
@@ -477 +506 @@
             psTrace("psModules.imcombine", 9, "Unable to find subregion for stamp (%d,%d)",
                     xPix, yPix);
-        }
+            pmSubtractionStampPrint(ds9, xPix, yPix, footprint, "yellow");
+        }
+    }
+
+    if (ds9) {
+        fclose(ds9);
     }
 
@@ -508 +542 @@
     }
 
-    if (psTraceGetLevel("psModules.imcombine") >= 9) {
-        pmSubtractionStampsRegions(stamps, "stamps_all.ds9");
-    }
-
     return stamps;
 }
@@ -581 +611 @@
 }
 
-#if 0
-bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int kernelSize)
-{
-    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
-    PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, false);
-    PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
-
-    int size = kernelSize + stamps->footprint; // Size of postage stamps
-    int num = stamps->num;              // Number of stamps
-    float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma
-
-    for (int i = 0; i < num; i++) {
-        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (!(stamp->status & PM_SUBTRACTION_STAMP_CALCULATE)) {
-            continue;
-        }
-
-        float x = stamp->x, y = stamp->y; // Coordinates of stamp
-        float flux = stamp->flux; // Flux of star
-        if (!isfinite(flux)) {
-            psWarning("Unable to generate PSF for stamp %d --- bad flux.", i);
-            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-            continue;
-        }
-
-        float xStamp = x - (int)(x + 0.5); // x coordinate of star in stamp frame
-        float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame
-
-        psFree(stamp->image2);
-        stamp->image2 = psKernelAlloc(-size, size, -size, size);
-        psKernel *target = stamp->image2; // Target stamp
-
-        // Put in a Waussian, just for fun!
-        for (int v = -size; v <= size; v++) {
-            for (int u = -size; u <= size; u++) {
-                float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma));
-                target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z));
-            }
-        }
-
-    }
-
-    return true;
-}
-#endif
-
 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *image,
                                                           const psImage *subMask, const psRegion *region,

trunk/psModules/src/imcombine/pmSubtractionStamps.h

@@ -125 +125 @@
 
 
-/// Write a ds9 region file with stamp positions
+/// Turn on/off generation of ds9 region files
 ///
 /// Intended for debugging
-bool pmSubtractionStampsRegions(pmSubtractionStampList *stamps, ///< Stamps
-                                const char *filename ///< Filename to which to write regions
+void pmSubtractionRegions(bool state    ///< Generate ds9 region files?
     );
 
+/// Open a file for ds9 regions
+///
+/// Intended for debugging
+FILE *pmSubtractionStampsFile(const pmSubtractionStampList *stamps, ///< List of stamps, for outlines
+                              const char *filename, ///< Filename to write
+                              const char *description ///< Description of file
+    );
+
+/// Print a stamp position to ds9 region file
+///
+/// Intended for debugging
+void pmSubtractionStampPrint(FILE *ds9, ///< ds9 region file
+                             float x, float y, ///< Position of stamp
+                             float size,///< Size of circle
+                             const char *color ///< Colour
+    );
 
 #endif