Changeset 21199 for trunk/ppStack/src/ppStackMatch.c

Timestamp:

Jan 28, 2009, 10:49:50 AM (17 years ago)

Author:

Paul Price

Message:

Merging pap_branch_20090108: fixed photometric normalisation. Conflicts resolved, code compiles.

File:

• trunk/ppStack/src/ppStackMatch.c (modified) (16 diffs)

trunk/ppStack/src/ppStackMatch.c

@@ -6 +6 @@
 #include <pslib.h>
 #include <psmodules.h>
+#include <psphot.h>
 
 #include "ppStack.h"
@@ -47 +48 @@
 #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 +81 @@
             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 +96 @@
             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
@@ -101 +120 @@
 }
 
-#define BG_SIZE 64                      // Large box half-size in which to measure background
-
-// Generate a background model of the readout we're matching
-psImage *stackBackgroundModel(pmReadout *ro, psImageMaskType maskVal, const psArray *sources, int size)
+// Add background into the fake image
+// Based on ppSubBackground()
+static psImage *stackBackgroundModel(pmReadout *ro, // Readout for which to generate background model
+                                     const pmConfig *config // Configuration
+    )
 {
-    psImage *image = ro->image, *mask = ro->mask; // Image and mask of readout
+    psAssert(ro && ro->image, "Need readout image");
+    psAssert(config, "Need configuration");
+
+    psImage *image = ro->image;         // Image of interest
     int numCols = image->numCols, numRows = image->numRows; // Size of image
 
-    psImage *bgImage = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Background image
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for measuring background
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
-
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i]; // Source of interest
-        if (!source) {
-            continue;
-        }
-
-        int x = source->peak->xf + 0.5, y = source->peak->yf + 0.5; // Coordinates of source
-
-        int xMin = PS_MAX(x - BG_SIZE, 0), xMax = PS_MIN(x + BG_SIZE, numCols);
-        int yMin = PS_MAX(y - BG_SIZE, 0), yMax = PS_MIN(y + BG_SIZE, numCols);
-
-        psRegion region = psRegionSet(xMin, xMax, yMin, yMax); // Region of interest
-        psImage *subImage = psImageSubset(image, region); // Subimage
-        psImage *subMask = psImageSubset(mask, region); // Subimage mask
-        if (!psImageBackground(stats, NULL, subImage, subMask, maskVal, rng)) {
-            // Can't do anything about it
-            psErrorClear();
-            continue;
-        }
-        psFree(subImage);
-        psFree(subMask);
-
-        float value = stats->robustMedian; // Value to set background to
-
-        int uMin = PS_MAX(x - size, 0), uMax = PS_MIN(x + size, numCols);
-        int vMin = PS_MAX(y - size, 0), vMax = PS_MIN(y + size, numCols);
-
-        for (int v = vMin; v < vMax; v++) {
-            for (int u = uMin; u < uMax; u++) {
-                bgImage->data.F32[v][u] = value;
-            }
-        }
-    }
-
-    psFree(stats);
-    psFree(rng);
-
-    return bgImage;
+    psMetadata *ppStackRecipe = psMetadataLookupPtr(NULL, config->recipes, PPSTACK_RECIPE);
+    psAssert(ppStackRecipe, "Need PPSTACK recipe");
+    psMetadata *psphotRecipe = psMetadataLookupPtr(NULL, config->recipes, PSPHOT_RECIPE);
+    psAssert(psphotRecipe, "Need PSPHOT recipe");
+
+    psString maskBadStr = psMetadataLookupStr(NULL, ppStackRecipe, "MASK.BAD");// Name of bits to mask for bad
+    psMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad);
+
+    psImage *binned = psphotBackgroundModel(ro, config); // Binned background model
+    psImageBinning *binning = psMetadataLookupPtr(NULL, psphotRecipe,
+                                                  "PSPHOT.BACKGROUND.BINNING"); // Binning for model
+    psAssert(binning, "Need binning parameters");
+    psImage *unbinned = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Unbinned background model
+    if (!psImageUnbin(unbinned, binned, binning)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to unbin background model");
+        psFree(unbinned);
+        return NULL;
+    }
+
+    // XXX should these really be here?? (probably not...)
+    // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL");
+    // pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL.STDEV");
+
+     return unbinned;
 }
+
 
 bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting,
@@ -163 +172 @@
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
     psAssert(recipe, "We've thrown an error on this before.");
+
+    float deconvLimit = psMetadataLookupF32(NULL, recipe, "DECONV.LIMIT"); // Limit on deconvolution fraction
 
     // Look up appropriate values from the ppSub recipe
@@ -195 +206 @@
         assert(outName);
         // Read convolution kernel
-        {
-            psString filename = NULL;   // Output filename
-            psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
-            psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
-            psFree(filename);
-            psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
-            psFree(resolved);
-            if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
-                psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
-                psFitsClose(fits);
-                return false;
-            }
+        psString filename = NULL;   // Output filename
+        psStringAppend(&filename, "%s.%d.kernel", outName, numInput);
+        psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename
+        psFree(filename);
+        psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel
+        psFree(resolved);
+        if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) {
+            psError(PS_ERR_IO, false, "Unable to read previously produced kernel");
             psFitsClose(fits);
-
-            // Add in variance factor
-            pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
-                                                                PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
-            float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
-            psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
-            if (!isfinite(vf)) {
-                vf = 1.0;
-            }
-            if (isfinite(vfItem->data.F32)) {
-                vfItem->data.F32 *= vf;
-            } else {
-                vfItem->data.F32 = vf;
-            }
+            return false;
+        }
+        psFitsClose(fits);
+
+        // Add in variance factor
+        pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis,
+                                                            PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels
+        float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor
+        psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR");
+        if (!isfinite(vf)) {
+            vf = 1.0;
+        }
+        if (isfinite(vfItem->data.F32)) {
+            vfItem->data.F32 *= vf;
+        } else {
+            vfItem->data.F32 = vf;
         }
 
@@ -244 +253 @@
         psFree(maskName);
         psFree(weightName);
+
+        psRegion *region = psMetadataLookupPtr(NULL, output->analysis,
+                                               PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region
+
+        pmSubtractionAnalysis(readout->analysis, kernels, region,
+                              readout->image->numCols, readout->image->numRows);
     } else {
 #endif
@@ -285 +300 @@
             }
 
-            // Generate a fake image to match to
-            float minFlux = INFINITY;       // Minimum flux for fake image
-            {
-                psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for bg
-                if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal, rng)) {
-                    psWarning("Can't measure background for image.");
-                    psErrorClear();
-                } else {
-                    minFlux = bg->robustStdev;
-                }
-                psFree(bg);
-            }
-
-#if 0
-            float maxMag = -INFINITY;   // Maximum magnitude of sources
-            for (int i = 0; i < sources->n; i++) {
-                pmSource *source = sources->data[i]; // Source of interest
-                if (source->psfMag > maxMag && source->psfMag <= MAG_IGNORE &&
-                    !(source->mode & SOURCE_MASK)) {
-                    maxMag = source->psfMag;
-                }
-            }
-            minFlux = PS_MIN(FAINT_SOURCE_FRAC * powf(10.0, -0.4 * maxMag), minFlux);
-#endif
-
-
 #if 0
             // Testing the normalisation of the fake image
@@ -318 +307 @@
                 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++) {
@@ -342 +331 @@
 
             if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
-                                          stampSources, NULL, NULL, psf, minFlux, 0, false, true)) {
+                                          stampSources, NULL, NULL, psf, NAN, footprint + size,
+                                          false, true)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF.");
                 psFree(fake);
@@ -351 +341 @@
 
             // Add the background into the target image
-            psImage *bgImage = stackBackgroundModel(readout, maskVal, stampSources, footprint); // Image of background
+            psImage *bgImage = stackBackgroundModel(readout, config); // Image of background
             psBinaryOp(fake->image, fake->image, "+", bgImage);
             psFree(bgImage);
 
-
 #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 +401 @@
 #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 +416 @@
                 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);
             }
@@ -548 +541 @@
     }
 
+    // Reject image completely if the maximum deconvolution fraction exceeds the limit
+    float deconv = psMetadataLookupF32(NULL, output->analysis,
+                                       PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Maximum deconvolution fraction
+    if (deconv > deconvLimit) {
+        psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n",
+                  deconv, deconvLimit);
+        psFree(output);
+        return NULL;
+    }
+
     // Renormalise the variances if desired
     if (renorm) {
@@ -580 +583 @@
     psFree(bg);
 
+
+#if 0
+#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
+    // XXX Ideally, would like to do proper PSF-fit photometry, but will settle for aperture photometry
+    int numSources = sources->n;        // Number of sources
+    psVector *ratio = psVectorAlloc(numSources, PS_TYPE_F32); // Flux ratios for sources
+    psVector *ratioMask = psVectorAlloc(numSources, PS_TYPE_MASK); // Mask for flux ratios
+    psVectorInit(ratioMask, 0xFF);
+    psImage *image = readout->image;    // Image of interest
+    psImage *mask = readout->mask;      // Mask of interest
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    for (int i = 0; i < numSources; i++) {
+        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->psfMag > MAX_MAG) {
+            continue;
+        }
+
+        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
+        int numPix = 0;                 // Number of pixels
+        float sum = 0.0;                // Sum of pixels
+        for (int y = yMin; y <= yMax; y++) {
+            for (int x = xMin; x <= xMax; x++) {
+                if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) {
+                    continue;
+                }
+                sum += image->data.F32[y][x];
+                numPix++;
+            }
+        }
+        if (sum >= 0 && numPix > 0) {
+            float mag = -2.5 * log10(sum * M_PI * PS_SQR(RADIUS) / numPix); // Instrumental magnitude
+            ratio->data.F32[i] = mag - source->psfMag;
+            ratioMask->data.PS_TYPE_MASK_DATA[i] = 0;
+        }
+    }
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
+    if (!psVectorStats(stats, ratio, NULL, ratioMask, 0xFF)) {
+        psWarning("Unable to measure normalisation --- assuming correct.");
+    } else {
+        psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image by %f (+/- %f) mag\n",
+                 stats->robustMedian, stats->robustStdev);
+        float norm = powf(10.0, -0.4 * stats->robustMedian); // Normalisation to apply
+        psBinaryOp(image, image, "*", psScalarAlloc(norm, PS_TYPE_F32));
+    }
+    psFree(stats);
+    psFree(ratio);
+    psFree(ratioMask);
+#endif
+
+#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, hdu->header, output->image, 0, NULL);
+        psFitsClose(fits);
+    }
+#endif
+
     psFree(output);