Changeset 25178 for branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c

Timestamp:

Aug 24, 2009, 8:40:34 AM (17 years ago)

Author:

eugene

Message:

various petrosian / radial profile analysis improvements

File:

• branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c (modified) (7 diffs)

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

@@ -12 +12 @@
 // finely spaced than r_{i+1} = r_i * \beta / \alpha.  for the integration, we need to
 // track the non-overlapping radius values.
-
-# define PETROSIAN_ALPHA 0.8
-# define PETROSIAN_BETA 1.25
 
 bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax) {
@@ -33 +30 @@
     psVector *radBet  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
 
-    psVector *fluxBin = psVectorAllocEmpty(nMax, PS_TYPE_F32);
-    psVector *radBin  = psVectorAllocEmpty(nMax, PS_TYPE_F32);
-    psVector *area    = psVectorAllocEmpty(nMax, PS_TYPE_F32);
+    psVector *binSB      = psVectorAllocEmpty(nMax, PS_TYPE_F32); // surface brightness of radial bin
+    psVector *binSBstdev = psVectorAllocEmpty(nMax, PS_TYPE_F32); // surface brightness of radial bin
+    psVector *binRad     = psVectorAllocEmpty(nMax, PS_TYPE_F32); // mean radius of radial bin
+    psVector *binArea    = psVectorAllocEmpty(nMax, PS_TYPE_F32); // area of radial bin (contiguous, non-overlapping)
 
-    psVectorInit (fluxBin, 0.0);
-    psVectorInit (radBin, 0.0);
+    psVectorInit (binSB, 0.0);
+    psVectorInit (binSBstdev, 0.0);
+    psVectorInit (binRad, 0.0);
 
     // generate radial bin bounds
@@ -56 +55 @@
     radBet->data.F32[2] = 2.0;
     
-    int nPts = 3;
-    for (int i = 3; i < radiusMax; i++) {
-        radMin->data.F32[nPts] = (i - 1);
-        radMax->data.F32[nPts] = i;
-        nPts++;
+# define PETROSIAN_ALPHA 0.8
+# define PETROSIAN_BETA 1.25
+# define POWER_LAW_SPACING true
+    
+    // power-law spacing with overlapping boundaries at the geometric mid-points
+    float rBeta = sqrt(PETROSIAN_BETA);
+    for (int i = 3; radBet->data.F32[i-1] < radiusMax; i++) {
+        if (POWER_LAW_SPACING) {
+            radMin->data.F32[i] = radMax->data.F32[i-1];
+            radMax->data.F32[i] = radMin->data.F32[i] * PETROSIAN_BETA;
+            radAlp->data.F32[i] = radMin->data.F32[i] / rBeta;
+            radBet->data.F32[i] = radMax->data.F32[i] * rBeta;
+        } else {
+            radMin->data.F32[i] = radMax->data.F32[i-1];
+            radMax->data.F32[i] = radMin->data.F32[i] + 1;
+            float rMid = 0.5*(radMin->data.F32[i] + radMax->data.F32[i]);
+            radAlp->data.F32[i] = rMid * PETROSIAN_ALPHA;
+            radBet->data.F32[i] = rMid * PETROSIAN_BETA;
+        }
+        radMin->n = radMax->n = radAlp->n = radBet->n = i + 1;
     }
-    radMin->n = radMax->n = radAlp->n = radBet->n = nPts;
 
     // generate radial area-weighted mean radius & non-overlapping areas
@@ -78 +91 @@
         
         // XXX calculate area-weighted radius rather than asserting?
-        radBin->data.F32[i] = rBin;
-        area->data.F32[i] = M_PI * (rMax2 - rMin2);
+        binRad->data.F32[i] = rBin;
+        binArea->data.F32[i] = M_PI * (rMax2 - rMin2);
 
-        if (i > 2) {
-            radAlp->data.F32[i] = rBin*PETROSIAN_ALPHA;
-            radBet->data.F32[i] = rBin*PETROSIAN_BETA;
+        if (0) {
+            fprintf (stderr, "%3d  %5.1f %5.1f : %5.1f : %5.1f %5.1f\n", 
+                 i, radAlp->data.F32[i], radMin->data.F32[i], binRad->data.F32[i],
+                 radMax->data.F32[i], radBet->data.F32[i]);
         }
     }
@@ -89 +103 @@
     // storage vector for stats
     psVector *values = psVectorAllocEmpty (flux->n, PS_TYPE_F32);
-    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
 
     // integrate flux, radius for each of these bins.  since flux is sorted by radius, 
@@ -107 +121 @@
             // calculate the value for the nOut bin
             psVectorStats (stats, values, NULL, NULL, 0);
-            fluxBin->data.F32[nOut] = stats->sampleMedian;
+            // binSB->data.F32[nOut] = stats->sampleMedian;
+            binSB->data.F32[nOut] = stats->robustMedian;
+            binSBstdev->data.F32[nOut] = stats->robustStdev / sqrt(values->n);
+
+            if (1) {
+                fprintf (stderr, "%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 ++;
-            if (nOut >= nMax) break;
+            if (nOut >= radAlp->n) break;
             Rmin = radAlp->data.F32[nOut];
             Rmax = radBet->data.F32[nOut];
@@ -122 +144 @@
         psVectorAppend (values, flux->data.F32[i]);
     }
-    fluxBin->n = radBin->n = area->n = nOut;
+    binSB->n = binSBstdev->n = binRad->n = binArea->n = nOut;
     // XXX I think this misses the last radial bin -- do we care?
 
     // save the vectors
-    petrosian->radialBins = radBin;
-    petrosian->area = area;
-    petrosian->binnedFlux = fluxBin;
+    petrosian->radialBins = binRad;
+    petrosian->area       = binArea;
+    petrosian->binSB      = binSB;
+    petrosian->binSBstdev = binSBstdev;
 
-    psphotPetrosianVisualProfileRadii (radius, flux, radBin, fluxBin, 0.0);
+    psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, 0.0);
 
     psFree(radMin);