Changeset 31154 for trunk/psphot/src/psphotRadialApertures.c

Timestamp:

Apr 4, 2011, 1:12:39 PM (15 years ago)

Author:

eugene

Message:

use the smoothed image for footprint culling; use a more stringent culling for saturated stars; more tweaks to get sat stars to be fitted; various updates to psphotStack; unify psphotImageLoop varients; be a bit careful about number of stars in a PSF clump

Location:

trunk/psphot

Files:

• . (modified) (2 props)
• src/psphotRadialApertures.c (modified) (12 diffs)

trunk/psphot

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM

trunk/psphot/src/psphotRadialApertures.c

@@ -37 +37 @@
     int Nradial = 0;
 
+    // perform full non-linear fits / extended source analysis?
+    if (!psMetadataLookupBool (&status, recipe, "RADIAL_APERTURES")) {
+        psLogMsg ("psphot", PS_LOG_INFO, "skipping radial apertures\n");
+        return true;
+    }
+
     psTimerStart ("psphot.radial");
 
@@ -62 +68 @@
     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)
@@ -67 +74 @@
     assert (maskVal);
 
-    // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
-    float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
-
     // S/N limit to perform full non-linear fits
     float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
@@ -75 +79 @@
     // source analysis is done in S/N order (brightest first)
     // XXX are we getting the objects out of order? does it matter?
-    sources = psArraySort (sources, pmSourceSortBySN);
+    sources = psArraySort (sources, pmSourceSortByFlux);
 
     // option to limit analysis to a specific region
@@ -95 +99 @@
         // limit selection to some SN limit
         assert (source->peak); // how can a source not have a peak?
-        if (source->peak->SN < SN_LIM) continue;
+        if (sqrt(source->peak->detValue) < SN_LIM) continue;
 
         // limit selection by analysis region
@@ -112 +116 @@
             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 ++;
 
         // 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);
-
-        if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, 0)) {
+        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
+
+        if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, 0)) {
             psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
         } else {
@@ -127 +135 @@
         }
 
+        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);
@@ -135 +148 @@
 }
 
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
+bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
 
     // if we are a child source, save the results to the parent source radial aperture array
@@ -163 +176 @@
     }
 
-    // center of the apertures
-    float xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
-    float yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+    float xCM = NAN, yCM = NAN;
+    if (pmSourcePositionUseMoments(source)) {
+        xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
+        yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+    } else {
+        xCM = source->peak->xf - 0.5 - source->pixels->col0; // coord of peak in subimage
+        yCM = source->peak->yf - 0.5 - source->pixels->row0; // coord of peak in subimage
+    }
 
     // one pass through the pixels to select the valid pixels and calculate R^2
@@ -203 +221 @@
     psVector *flux    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
     psVector *fluxErr = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
+    psVector *fluxStd = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
     psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
 
     psVectorInit (flux,    0.0);
+    psVectorInit (fluxStd, 0.0);
     psVectorInit (fluxErr, 0.0);
     psVectorInit (fill,    0.0);
@@ -212 +232 @@
     for (int i = 0; i < pixRadius2->n; i++, rPix2++) {
 
+        int j = 0;
         float *aRad2 = aperRadii2->data.F32;
-        for (int j = 0; (*rPix2 < *aRad2) && (j < aperRadii2->n); j++, aRad2++) {
+        for (; (*aRad2 < *rPix2) && (j < aperRadii2->n); j++, aRad2++);
+        for (; j < aperRadii2->n; j++, aRad2++) {
             flux->data.F32[j]    += pixFlux->data.F32[i];
+            fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]);
             fluxErr->data.F32[j] += pixVar->data.F32[i];
             fill->data.F32[j]    += 1.0;
         }
     }
+
+    /* for each radial bin, R(i), we measure:
+       1) the flux within that aperture: F(i) = \sum_{r_j<R_i}(F_j)
+       2) the fractional fill factor (count of valid pixels / effective area of the aperture
+       3) the error on the flux within that aperture
+     */
 
     for (int i = 0; i < flux->n; i++) {
         // calculate the total flux for bin 'nOut'
         float Area = M_PI*aperRadii2->data.F32[i];
-        fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
-        fill->data.F32[i] /= Area;
-        psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %4.2f %6.1f\n", 
-                 i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], fill->data.F32[i], Area);
+
+        int nPix = fill->data.F32[i];
+        float SBmean = flux->data.F32[i] / nPix;
+        float SBstdv = sqrt((fluxStd->data.F32[i] / nPix) - PS_SQR(SBmean));
+
+        // flux->data.F32[i]    = SBmean * Area;
+        fluxStd->data.F32[i] = SBstdv * Area;
+        fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
+
+        // fill->data.F32[i] /= Area;
+        psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %8.1f +/- %8.1f : %4.2f %6.1f\n", 
+                 i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], SBmean, SBstdv, fill->data.F32[i], Area);
     }
     
     radialAper->flux = flux;
+    radialAper->fluxStdev = fluxStd;
     radialAper->fluxErr = fluxErr;
     radialAper->fill = fill;
@@ -245 +283 @@
 static int nPix = 0;
 
-bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry) {
+bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *radMax, int entry) {
 
     psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?");

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM