Changeset 30101 for branches/eam_branches/ipp-20101205/psphot/src/psphotRadialAperturesByObject.c

Timestamp:

Dec 17, 2010, 9:57:34 AM (16 years ago)

Author:

eugene

Message:

various modifications to psphotStack: add needed metadata to header; enable multiple PSF-matched images; work to speed up radial aperture photometry; radialAper is now an array (one element per PSF size); better feedback on extended fit skips; clarify processing: main analysis is on the input images (raw or convolved) while radial apertures is on the psf-matched images

File:

• branches/eam_branches/ipp-20101205/psphot/src/psphotRadialAperturesByObject.c (modified) (7 diffs)

branches/eam_branches/ipp-20101205/psphot/src/psphotRadialAperturesByObject.c

@@ -3 +3 @@
 // aperture-like measurements for extended sources
 // flux in simple, circular apertures
-bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule) {
+
+// **** it looks like this function will re-point the source pixels at the specified FILERULE
+// **** I need to distinguish PSF-matched images from raw
+// **** save (somewhere) the PSF-matched PSF values
+
+// this function measures the radial aperture fluxes for the set of readouts.  this function
+// may be called multiple times (presumably with different matched PSF sizes).  we must have
+// already added an entry to the readout->analysis identifying the FWHM of this version.
+
+bool psphotRadialAperturesByObject (pmConfig *config, psArray *objects, const pmFPAview *view, const char *filerule, int nMatchedPSF) {
 
     bool status;
@@ -28 +37 @@
     psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
     psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
+    float outerRadius = radMax->data.F32[radMax->n - 1];
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -39 +49 @@
     float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
+    // how many target PSFs do we want?
+    int nPSFsizes = 0;
+    { 
+        psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM");
+        if (status) {
+            nPSFsizes = 1;
+        } else {
+            psVector *fwhmValues = psMetadataLookupVector(&status, recipe, "PSPHOT.STACK.TARGET.PSF.FWHM"); // Magnitude offsets
+            psAssert (status, "missing psphot recipe value PSPHOT.STACK.TARGET.PSF.FWHM");
+            nPSFsizes = fwhmValues->n;
+        }
+    }
+    
     // source analysis is done in S/N order (brightest first)
     objects = psArraySort (objects, pmPhotObjSortBySN);
@@ -53 +76 @@
         psAssert (readout, "missing readout?");
 
+        psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+        if (!fwhmValues) {
+            psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM is not defined for this readout");
+            return false;
+        }
+        if (fwhmValues->n != nMatchedPSF + 1) {
+            psError (PSPHOT_ERR_CONFIG, true, "convolved or measured FWHM sequence is inconsistent this readout");
+            return false;
+        }
+        psLogMsg ("psphot", PS_LOG_DETAIL, "PSF FWHM of %s : %f pixels\n", file->name, fwhmValues->data.F32[nMatchedPSF]);
+
         readouts->data[i] = psMemIncrRefCounter(readout);
     }
 
     // process the objects in order.  
-    for (int i = 0; i < objects->n; i++) {
+    // XXX TEST: fix this!
+    for (int i = 0; (i < 50) && (i < objects->n); i++) {
         pmPhotObj *object = objects->data[i];
         if (!object) continue;
@@ -77 +112 @@
             if (source->peak->SN < SN_LIM) continue;
 
+            int index = source->imageID;
+            if (index >= readouts->n) continue; // skip the sources generated by the chisq image
+            pmReadout *readout = readouts->data[index];
+
+            // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
+            // psphotVisualShowImage(readout);
+
+            // allocate pmSourceExtendedParameters, if not already defined
+            if (!source->radialAper) {
+                source->radialAper = psArrayAlloc(nPSFsizes);
+            }
+
             // replace object in image
             if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
@@ -83 +130 @@
             Nradial ++;
 
-            int index = source->imageID;
-            pmReadout *readout = readouts->data[index];
+            // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
+            // psphotVisualShowImage(readout);
 
             // force source image to be a bit larger...
-            float radius = source->peak->xf - source->pixels->col0;
-            radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
-            radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
-            radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
-            pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
+            // float radius = source->peak->xf - source->pixels->col0;
+            // radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
+            // radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
+            // radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
+            pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
 
-            if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) {
+            if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, nMatchedPSF)) {
                 psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
             } else {
@@ -101 +148 @@
             // re-subtract the object, leave local sky
             pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+
+            // psLogMsg("psphot", PS_LOG_INFO, "radial apertures for %d", index);
+            // psphotVisualShowImage(readout);
         }
     }