Changeset 20427 for trunk/ppStack/src/ppStackPhotometry.c

Timestamp:

Oct 27, 2008, 4:27:35 PM (18 years ago)

Author:

Paul Price

Message:

Adding function to re-do photometry of input sources

File:

• trunk/ppStack/src/ppStackPhotometry.c (modified) (1 diff)

trunk/ppStack/src/ppStackPhotometry.c

@@ -9 +9 @@
 
 #include "ppStack.h"
+
+#define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \
+                          PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \
+                          PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources
+
+bool ppStackInputPhotometry(const pmReadout *ro, const psArray *sources, const pmConfig *config)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, false);
+    PS_ASSERT_ARRAY_NON_NULL(sources, false);
+
+    psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
+    psAssert(recipe, "We've thrown an error on this before.");
+
+    bool mdok;                          // Status of MD lookup
+
+    float zpRadius = psMetadataLookupS32(&mdok, recipe, "ZP.RADIUS"); // Radius for ZP measurement
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, true, "Unable to find ZP.RADIUS in recipe");
+        return false;
+    }
+    float zpSigma = psMetadataLookupF32(&mdok, recipe, "ZP.SIGMA"); // Gaussian sigma for photometry
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, true, "Unable to find ZP.SIGMA in recipe");
+        return false;
+    }
+    float zpFrac = psMetadataLookupF32(&mdok, recipe, "ZP.FRAC"); // Fraction of good pixels for photometry
+    if (!mdok) {
+        psError(PS_ERR_UNKNOWN, true, "Unable to find ZP.FRAC in recipe");
+        return false;
+    }
+
+    psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in
+    if (!mdok || !maskValStr) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe");
+        return false;
+    }
+    psMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask
+
+    psImage *image = ro->image, *mask = ro->mask; // Image and mask from readout
+
+    // Measure background
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
+    if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image");
+        psFree(stats);
+        psFree(rng);
+        return false;
+    }
+    float bg = stats->robustMedian; // Background level
+    psFree(stats);
+    psFree(rng);
+
+    // Photometer sources
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    int numSources = sources->n; // Number of sources
+    int numGood = 0;            // Number of good sources
+    float maxMag = -INFINITY;   // Maximum magnitude
+    for (int j = 0; j < numSources; j++) {
+        pmSource *source = sources->data[j]; // Source of interest
+        if (source->mode & PHOT_SOURCE_MASK || !isfinite(source->psfMag)) {
+            source->psfMag = NAN;
+            continue;
+        }
+        float x = source->peak->xf, y = source->peak->yf; // Coordinates of source
+        int xCentral = x + 0.5, yCentral = y + 0.5; // Central pixel
+        // Range of integration
+        int xMin = PS_MAX(0, xCentral - zpRadius), xMax = PS_MIN(numCols - 1, xCentral + zpRadius);
+        int yMin = PS_MAX(0, yCentral - zpRadius), yMax = PS_MIN(numRows - 1, yCentral + zpRadius);
+
+        // Integrate
+        double sumImage = 0.0, sumKernel = 0.0; // Sums from integration
+        int numBadPix = 0;         // Number of bad pixels
+        for (int v = yMin; v <= yMax; v++) {
+            float dy2 = PS_SQR(y - v); // Distance from centroid
+            for (int u = xMin; u <= xMax; u++) {
+                if (mask->data.PS_TYPE_MASK_DATA[v][u] & maskVal) {
+                    numBadPix++;
+                    continue;
+                }
+                float dx2 = PS_SQR(x - u); // Distance from centroid
+                double kernel = exp(-0.5 * (dx2 + dy2) / PS_SQR(zpSigma)); // Kernel value
+                sumImage += (image->data.F32[v][u] - bg) * kernel;
+                sumKernel += kernel;
+            }
+        }
+
+        if (numBadPix > zpFrac * M_PI * PS_SQR(zpRadius) || !isfinite(sumImage)) {
+            source->psfMag = NAN;
+        } else {
+            source->psfMag = - 2.5 * log10(sumImage / sumKernel * M_PI * PS_SQR(zpRadius));
+            numGood++;
+            maxMag = PS_MAX(maxMag, source->psfMag);
+        }
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Photometered %d good sources in image; max mag %f", numGood, maxMag);
+
+    return true;
+}
+
+
 
 bool ppStackPhotometry(pmConfig *config, const pmReadout *readout, const pmFPAview *view)