Changeset 21106 for branches/pap_branch_20090108/ppStack/src/ppStackMatch.c

Timestamp:

Jan 12, 2009, 1:03:29 PM (17 years ago)

Author:

Paul Price

Message:

Adding more debugging output. Think I've solved the problem of the poor photometry: the background being added in to the fake images produces the bad photometry. Without it, the photometry s.d. goes from 0.08 mag to 0.02 mag. Want to replace it with the proper psphot background model (not done yet). These files are currently in an intermediate state.

File:

• branches/pap_branch_20090108/ppStack/src/ppStackMatch.c (modified) (11 diffs)

branches/pap_branch_20090108/ppStack/src/ppStackMatch.c

@@ -47 +47 @@
 #endif
 
+// Get coordinates from a source
+static void coordsFromSource(float *x, float *y, const pmSource *source)
+{
+    assert(x && y);
+    assert(source);
+
+    if (source->modelPSF) {
+        *x = source->modelPSF->params->data.F32[PM_PAR_XPOS];
+        *y = source->modelPSF->params->data.F32[PM_PAR_YPOS];
+    } else {
+        *x = source->peak->xf;
+        *y = source->peak->yf;
+    }
+    return;
+}
+
 static psArray *stackSourcesFilter(psArray *sources, // Source list to filter
                                    int exclusion // Exclusion zone, pixels
@@ -64 +80 @@
             continue;
         }
-        x->data.F32[numGood] = source->peak->xf;
-        y->data.F32[numGood] = source->peak->yf;
+        coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source);
         numGood++;
     }
@@ -80 +95 @@
             continue;
         }
-        coords->data.F64[0] = source->peak->xf;
-        coords->data.F64[1] = source->peak->yf;
+        float xSource, ySource;         // Coordinates of source
+        coordsFromSource(&xSource, &ySource, source);
+
+        coords->data.F64[0] = xSource;
+        coords->data.F64[1] = ySource;
 
         long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone
@@ -120 +138 @@
         }
 
-        int x = source->peak->xf + 0.5, y = source->peak->yf + 0.5; // Coordinates of source
+        float xSource, ySource;         // Coordinates of source
+        coordsFromSource(&xSource, &ySource, source);
+
+        int x = xSource + 0.5, y = ySource + 0.5; // Coordinates of source
 
         int xMin = PS_MAX(x - BG_SIZE, 0), xMax = PS_MIN(x + BG_SIZE, numCols);
@@ -319 +340 @@
                 pmSource *source = pmSourceAlloc();     // Source
                 sources->data[0] = source;
-                source->peak = pmPeakAlloc(50, 50, 10000, PM_PEAK_LONE);
+                source->peak = pmPeakAlloc(500, 500, 10000, PM_PEAK_LONE);
                 source->psfMag = -13.0;
-                pmReadoutFakeFromSources(test, 100, 100, sources, NULL, NULL, psf, 0.1, 0, false, true);
+                pmReadoutFakeFromSources(test, 1000, 1000, sources, NULL, NULL, psf, 0.1, 0, false, true);
                 float sum = 0.0;
                 for (int y = 0; y < test->image->numRows; y++) {
@@ -351 +372 @@
             }
 
+
+#if 0
             // Add the background into the target image
             psImage *bgImage = stackBackgroundModel(readout, maskVal, stampSources, footprint); // Image of background
             psBinaryOp(fake->image, fake->image, "+", bgImage);
             psFree(bgImage);
+#endif
+
 
 #ifdef TESTING
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "fake_%03d.fits", numInput);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
-                psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
                 psFitsClose(fits);
             }
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "real_%03d.fits", numInput);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
-                psFitsWriteImage(fits, NULL, readout->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL);
                 psFitsClose(fits);
             }
@@ -410 +437 @@
 #ifdef TESTING
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "conv_%03d.fits", numInput);
                 psFits *fits = psFitsOpen(name, "w");
                 psFree(name);
-                psFitsWriteImage(fits, NULL, output->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
                 psFitsClose(fits);
             }
             {
+                pmHDU *hdu = pmHDUFromCell(readout->parent);
                 psString name = NULL;
                 psStringAppend(&name, "diff_%03d.fits", numInput);
@@ -423 +452 @@
                 psFree(name);
                 psBinaryOp(fake->image, output->image, "-", fake->image);
-                psFitsWriteImage(fits, NULL, fake->image, 0, NULL);
+                psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL);
                 psFitsClose(fits);
             }
@@ -582 +611 @@
 #define RADIUS 10                       // Radius of photometry
 #define MIN_ERR 0.05                    // Minimum photometric error, mag
+#define MAX_MAG -13                     // Maximum magnitude for source
 
     // Ensure the normalisation is correct
@@ -595 +625 @@
         pmSource *source = sources->data[i]; // Source of interest
         if (!source || source->mode & SOURCE_MASK || !isfinite(source->psfMag) || !isfinite(source->errMag) ||
-            source->errMag > MIN_ERR) {
+            source->errMag > MIN_ERR || source->psfMag > MAX_MAG) {
             continue;
         }
 
-        float xSrc = source->peak->xf, ySrc = source->peak->yf; // Source coordinates
+        float xSrc, ySrc;              // Source coordinates
+        coordsFromSource(&xSrc, &ySrc, source);
         int xMin = PS_MAX(0, xSrc - RADIUS), xMax = PS_MIN(numCols - 1, xSrc + RADIUS); // Bounds in x
         int yMin = PS_MAX(0, ySrc - RADIUS), yMax = PS_MIN(numRows - 1, ySrc + RADIUS); // Bounds in y
@@ -635 +666 @@
 #ifdef TESTING
     {
+        pmHDU *hdu = pmHDUFromCell(readout->parent);
         psString name = NULL;
         psStringAppend(&name, "convolved_%03d.fits", numInput);
         psFits *fits = psFitsOpen(name, "w");
         psFree(name);
-        psFitsWriteImage(fits, NULL, output->image, 0, NULL);
+        psFitsWriteImage(fits, hdu->header, output->image, 0, NULL);
         psFitsClose(fits);
     }