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

Timestamp:

Mar 17, 2009, 10:00:02 AM (17 years ago)

Author:

Paul Price

Message:

Reorganising ppStackLoop: splitting into multiple separate functions which are joined using a common data carrier, ppStackOptions (which aren't really 'options', but that's the word I chose...).

File:

• trunk/ppStack/src/ppStackPhotometry.c (modified) (5 diffs)

trunk/ppStack/src/ppStackPhotometry.c

@@ -9 +9 @@
 
 #include "ppStack.h"
+#include "ppStackLoop.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)
+bool ppStackPhotometry(ppStackOptions *options, pmConfig *config)
 {
-    PM_ASSERT_READOUT_NON_NULL(ro, false);
-    PS_ASSERT_ARRAY_NON_NULL(sources, false);
+    psAssert(options, "Require options");
+    psAssert(config, "Require configuration");
 
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe
@@ -23 +20 @@
 
     bool mdok;                          // Status of MD lookup
-
-    float zpRadius = psMetadataLookupS32(&mdok, recipe, "PHOT.RADIUS"); // Radius for PHOT measurement
-    if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.RADIUS in recipe");
-        return false;
-    }
-    float zpSigma = psMetadataLookupF32(&mdok, recipe, "PHOT.SIGMA"); // Gaussian sigma for photometry
-    if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.SIGMA in recipe");
-        return false;
-    }
-    float zpFrac = psMetadataLookupF32(&mdok, recipe, "PHOT.FRAC"); // Fraction of good pixels for photometry
-    if (!mdok) {
-        psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.FRAC in recipe");
-        return false;
+    if (!psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) {
+        // Nothing to do
+        return true;
     }
 
-    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;
-    }
-    psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask
+    psTimerStart("PPSTACK_PHOT");
 
-    psImage *image = ro->image, *mask = ro->mask; // Image and mask from readout
+    ppStackFileActivation(config, PPSTACK_FILES_COMBINE, false);
+    ppStackFileActivation(config, PPSTACK_FILES_PHOT, true);
+    pmFPAview *photView = ppStackFilesIterateDown(config); // View to readout
+    ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true);
 
-    // Measure background
-    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // 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_IMAGE_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));
-            if (isfinite(source->psfMag)) {
-                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)
-{
-    pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT");
-    pmFPACopy(photFile->fpa, readout->parent->parent->parent);
+    pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File for photometry
+    pmFPACopy(photFile->fpa, options->outRO->parent->parent->parent);
 
     if (psMetadataLookupBool(NULL, config->arguments, "-visual")) {
@@ -122 +39 @@
     }
 
+#if 0
     // Need to ensure aperture residual is not calculated
     psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe
@@ -131 +49 @@
         item->data.B = false;
     }
+#endif
 
     // set maskValue and markValue in the psphot recipe
     psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask
     psImageMaskType markValue = pmConfigMaskGet("MARK.VALUE", config); // Bits to use for marking
-    psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "Bits to mask", maskValue);
-    psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE, "Bits to use for mark", markValue);
+    psMetadataAddImageMask(recipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE,
+                            "Bits to mask", maskValue);
+    psMetadataAddImageMask(recipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE,
+                           "Bits to use for mark", markValue);
 
-    // XXX replace with psphotReadoutKnownSources (config, view)
-    // XXX this function requires that we construct the input source list and place it on PSPHOT.INPUT.CMF
-    pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT");
-    psArray *inSources = psMetadataLookupPtr (NULL, sourcesCell->analysis, "PSPHOT.SOURCES");
+    pmCell *sourcesCell = pmFPAfileThisCell(config->files, photView, "PPSTACK.OUTPUT");
+    psArray *inSources = psMetadataLookupPtr(NULL, sourcesCell->analysis, "PSPHOT.SOURCES");
     if (!inSources) {
         psError(PS_ERR_UNKNOWN, false, "Unable to find input sources");
+        psFree(photView);
         return false;
     }
 
-    if (!psphotReadoutKnownSources(config, view, inSources)) {
+    if (!psphotReadoutKnownSources(config, photView, inSources)) {
         // Clear the error, so that the output files are written.
         psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on stacked image.");
+        psFree(photView);
         return false;
     }
+    psFree(photView);
 
     if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") ||
@@ -162 +84 @@
     pmFPAfileActivate(config->files, false, "PSPHOT.INPUT");
 
+    if (options->stats) {
+        pmReadout *photRO = pmFPAviewThisReadout(photView, photFile->fpa); // Readout with the sources
+        psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources
+        psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0,
+                         "Number of sources detected", sources->n);
+        psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0,
+                         "Time to do photometry", psTimerMark("PPSTACK_PHOT"));
+    }
+    psLogMsg("ppStack", PS_LOG_INFO, "Time to do photometry: %f sec", psTimerClear("PPSTACK_PHOT"));
+
     return true;
 }