Changeset 25128 for branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c

Timestamp:

Aug 19, 2009, 12:16:54 PM (17 years ago)

Author:

eugene

Message:

updates to make radial profiles a bit more robust

File:

• branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c (modified) (3 diffs)

branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c

@@ -30 +30 @@
 
     // examine data in the two ranges Fm - Fo and Fo - Fp to define the bin size
+    // XXX reconsider the fractional isophote value
     float Fm = fluxMin + 0.25*fluxRange;
     float Fp = fluxMin + 0.75*fluxRange;
     float Fo = fluxMin + 0.50*fluxRange;
+    int Rbin = 1;
       
-    // find the mean radius of the points in the flux range Fm - Fp:
-    float Rsum = 0;
-    float Rnpt = 0;
-    for (int i = 0; i < flux->n; i++) {
-        if (flux->data.F32[i] < Fm) continue;
-        if (flux->data.F32[i] > Fp) continue;
-      
-        Rsum += radius->data.F32[i];
-        Rnpt ++;
-    }
+    // find the median radius of the points in the flux range Fm - Fp:
+    { 
+        // storage vector for stats
+        psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
+        psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
 
-    // rebin with bins 1/2 the size of the mean radius
-    int Rbin = 0;
-    if (Rnpt == 0) {
-        Rbin = 1;
-    } else {
-        Rbin = MAX(1, 0.5*(Rsum / Rnpt));
+        for (int i = 0; i < flux->n; i++) {
+            if (!isfinite(flux->data.F32[i])) continue;
+            if (flux->data.F32[i] < Fm) continue;
+            if (flux->data.F32[i] > Fp) continue;
+            
+            psVectorAppend (values, radius->data.F32[i]);
+        }
+        psVectorStats (stats, values, NULL, NULL, 0);
+
+        // if we have a valid range, rebin with bin size 1/2 of median radius
+        if (isfinite(stats->sampleMedian)) {
+            Rbin = MAX(1, 0.5*stats->sampleMedian);
+        }
     }
 
@@ -94 +98 @@
                 psStatsInit(stats);
             }
+            if (!isfinite(flux->data.F32[i])) continue;
             psVectorAppend (values, flux->data.F32[i]);
         }
@@ -105 +110 @@
     bool above = true;
     for (int i = 0; i < fluxBinned->n; i++) {
+
+        if (!isfinite(fluxBinned->data.F32[i])) continue;
 
         // find the largest radius that matches the flux transition