Changeset 27838 for branches/tap_branches/psphot/src/psphotPetrosianRadialBins.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotPetrosianRadialBins.c (modified) (7 diffs)

branches/tap_branches/psphot/src

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

branches/tap_branches/psphot/src/psphotPetrosianRadialBins.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
 
 // convert the flux vs elliptical radius to annular bins
@@ -5 +6 @@
 // we are guaranteed to be limited by either the seeing (1 - few pixels) or by the pixels
 // themselves.  this function does not attempt to measure the radial profiles accurately
-// for radii that are smaller than a minimum (currently 1.0 pixels).  
+// for radii that are smaller than ~2 pixels
 
 // for small radii, we are measuring the mean surface brightness in non-overlapping radial
@@ -13 +14 @@
 // track the non-overlapping radius values.
 
+// Photo interpolates the image of interest to place the peak on the center of the central
+// pixel, and then uses the exact fractions of the pixels in each of the first few annuli.
+// Seems like a reasonable thing, but is there any significance to the difference?
+
 // XXX move the resulting elements from profile to extpars->petrosian?
 bool psphotPetrosianRadialBins (pmSource *source, float radiusMax, float skynoise) {
 
-    pmSourceRadialProfile *profile = source->extpars->profile;
-
-    float skyModelErrorSQ = PS_SQR(skynoise);
-
-    psVector *radius = profile->radiusElliptical;
-    psVector *flux = profile->fluxElliptical;
+    psAssert (source, "missing source");
+    psAssert (source->extpars, "missing extpars");
+    psAssert (source->extpars->ellipticalFlux, "missing ellipticalFlux");
+
+    psVector *radius = source->extpars->ellipticalFlux->radiusElliptical;
+    psVector *flux = source->extpars->ellipticalFlux->fluxElliptical;
 
     // sort incoming vectors by radius
     pmSourceRadialProfileSortPair (radius, flux);
+
+    if (!source->extpars->petProfile) {
+        source->extpars->petProfile = pmSourceRadialProfileAlloc();
+    }
+    pmSourceRadialProfile *profile = source->extpars->petProfile;
+
+    float skyModelErrorSQ = PS_SQR(skynoise);
 
     int nMax = radiusMax;
@@ -107 +119 @@
     psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    // psStats *stats = psStatsAlloc(PS_STAT_FITTED_MEAN_V4 | PS_STAT_FITTED_STDEV_V4);
-
-    // integrate flux, radius for each of these bins.  since flux is sorted by radius, 
-    // we can do this fairly quickly
 
     bool done = false;
@@ -133 +141 @@
                 dvalue = NAN;
             }
-            // binSB->data.F32[nOut] = stats->sampleMedian;
+
             binSB->data.F32[nOut] = value;
             binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ);
-            // binSB->data.F32[nOut] = stats->fittedMean;
-            // binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->fittedStdev) / values->n + skyModelErrorSQ);
 
             // error in the SB is the stdev per bin / sqrt (number of pixels) 
@@ -143 +149 @@
             // residual flux, but the sky from the sky model)
 
-            psTrace ("psphot", 5, "%3d  %5.1f %5.1f : %5.1f  %5.2f\n", 
-                     nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]);
+            psTrace ("psphot", 5, "%3d  %5.1f %5.1f : %5.1f  %5.2f\n", nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]);
 
             nOut ++;
@@ -163 +168 @@
     // XXX I think this misses the last radial bin -- do we care?
 
+    // interpolate any bins that were empty (extrapolate to center if needed)
+    if (!isfinite(binSB->data.F32[0]) && !isfinite(binSB->data.F32[1])) {
+        psWarning ("center 2 bins of source at %f, %f are NAN, skipping this source", source->peak->xf, source->peak->yf);
+        // XXX raise a flag
+        psFree(binSB);
+        psFree(binSBstdev);
+        psFree(binRad);
+        psFree(binArea);
+        psFree(radMin);
+        psFree(radMax);
+        psFree(radAlp);
+        psFree(radBet);
+        psFree(values);
+        psFree(stats);
+        return false;
+    }
+
+    // if center bin is empty assume same SB as next radius (probably true due to PSF)
+    if (!isfinite(binSB->data.F32[0])) {
+        binSB->data.F32[0] = binSB->data.F32[1];
+        binSBstdev->data.F32[0] = binSBstdev->data.F32[1];
+    }
+
+    // interpolate any bins that were empty (if center if needed)
+    for (int i = 1; i < binSB->n - 1; i++) {
+        if (isfinite(binSB->data.F32[i])) continue;
+        binSB->data.F32[i] = InterpolateValues (binRad->data.F32[i-1], binSB->data.F32[i-1], binRad->data.F32[i+1], binSB->data.F32[i+1], binRad->data.F32[i]);
+        binSBstdev->data.F32[i] = InterpolateValues (binRad->data.F32[i-1], binSBstdev->data.F32[i-1], binRad->data.F32[i+1], binSBstdev->data.F32[i+1], binRad->data.F32[i]);
+    }
+
+    psFree(profile->binSB);
+    psFree(profile->binSBstdev);
+    psFree(profile->radialBins);
+    psFree(profile->area);
+
     // save the vectors
     profile->radialBins = binRad;

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack