Changeset 30631 for branches/czw_branch/20101203/psphot/src/psphotRadialAperturesByObject.c

Timestamp:

Feb 14, 2011, 1:05:28 PM (15 years ago)

Author:

watersc1

Message:

Merge of trunk back into branch.

Location:

branches/czw_branch/20101203

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src/psphotRadialAperturesByObject.c (modified) (6 diffs)

branches/czw_branch/20101203/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);
@@ -52 +75 @@
         pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
         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);
@@ -77 +111 @@
             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->parent) {
+                if (!source->parent->radialAper) {
+                    source->parent->radialAper = psArrayAlloc(nPSFsizes);
+                }
+            } else {
+                if (!source->radialAper) {
+                    source->radialAper = psArrayAlloc(nPSFsizes);
+                }
+            }
+
             // replace object in image
             if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
                 pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
             }
+
+            // we need to change the view for the radial aperture analysis, but we want to recover exactly
+            // the original view; the following elements get destroyed by pmSourceRedefinePixels so save them:
+            psImage *oldMaskObj   = psMemIncrRefCounter(source->maskObj);
+            psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux);
+            psImage *oldPSFimage  = psMemIncrRefCounter(source->psfImage);
+            psRegion oldRegion    = source->region;
+
             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 {
@@ -99 +159 @@
             }
 
+            pmSourceRedefinePixelsByRegion (source, readout, oldRegion);
+            psFree(source->maskObj);   source->maskObj   = oldMaskObj;
+            psFree(source->modelFlux); source->modelFlux = oldModelFlux;
+            psFree(source->psfImage);  source->psfImage  = oldPSFimage;
+
             // 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);
         }
     }