Changeset 25128 for branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.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/psphotRadialProfileByAngle.c (modified) (3 diffs)

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

@@ -5 +5 @@
 
 // XXX choices for rMax and Nsec?
-// XXX where does the output go?
-// XXX how does this fit in context?
 
 bool psphotRadialProfilesByAngles (pmPetrosian *petrosian, pmSource *source, int Nsec, float Rmax) {
 
+    // we want to have an even number of sectors so we can do 180 deg symmetrizing
+    Nsec = (Nsec % 2) ? Nsec + 1 : Nsec;
     float dtheta = 2.0*M_PI / Nsec;
 
@@ -46 +46 @@
             // value is NAN if we run off the image
             float value = psImageInterpolatePixelBilinear(x, y, source->pixels);
-            if (isnan(value)) continue;
+
+            // keep the nan values so all vectors are matched
+            // if (isnan(value)) continue;
 
             psVectorAppend (radius, r);
@@ -62 +64 @@
     }
 
+    for (int i = 0; i < Nsec / 2; i++) {
+
+        psVector *r1 = petrosian->radii->data[i];
+        psVector *r2 = petrosian->radii->data[i+Nsec/2];
+
+        psVector *f1 = petrosian->fluxes->data[i];
+        psVector *f2 = petrosian->fluxes->data[i+Nsec/2];
+
+        psAssert (r1->n == r2->n, "mis-matched vectors");
+        psAssert (f1->n == f2->n, "mis-matched vectors");
+
+        // we have a pair of vectors i, i+Nsec/2; replace them with the finite minimum of the pair
+        for (int j = 0; j < r1->n; j++) {
+            
+            float flux;
+
+            if (!isfinite(f1->data.F32[j]) && !isfinite(f2->data.F32[j])) {
+                flux = NAN;
+                goto setflux;
+            }
+
+            if (!isfinite(f1->data.F32[j])) {
+                flux = f2->data.F32[j];
+                goto setflux;
+            }
+            if (!isfinite(f2->data.F32[j])) {
+                flux = f1->data.F32[j];
+                goto setflux;
+            }
+
+            flux = PS_MIN(f1->data.F32[j], f2->data.F32[j]);
+
+        setflux:
+            f1->data.F32[j] = flux;
+            f2->data.F32[j] = flux;
+        }
+    }    
     return true;
 }