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

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

@@ -4 +4 @@
 // separations
 
-// XXX choices for rMax and Nsec?
+float psphotMeanSectorValue (psImage *image, int x, int y, float dL, float dW, float theta);
+psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta);
+psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW);
+psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords);
 
 bool psphotRadialProfilesByAngles (pmPetrosian *petrosian, pmSource *source, int Nsec, float Rmax) {
@@ -28 +31 @@
         psVector *flux   = psVectorAllocEmpty(Rmax, PS_TYPE_F32);
 
-        // start at Xo,Yo and find the x,y locations for r_i, theta where r_i increments by 1 pixel
-        // XXX at large radii ( > 10-15) use stats in a patch rather than sub-pixel interpolation
-        for (float r = 0; r < Rmax; r += 1.0) {
+        // Start at Xo,Yo and find the x,y locations for r_i, theta where r_i initially
+        // increments by 1 pixel.  At large radii (r*dtheta > 2) use stats in a box rather than
+        // sub-pixel interpolation
+
+        int dR = 1.0;
+        for (float r = 0; r < Rmax; r += dR) {
 
             float Xo = source->peak->xf + 0.5;
@@ -38 +44 @@
             float x = r * cos (theta) + Xo;
             float y = r * sin (theta) + Yo;
-
-            if (x < 0) continue;
-            if (y < 0) continue;
-            if (x >= source->pixels->parent->numCols) continue;
-            if (y >= source->pixels->parent->numRows) continue;
-
-            // value is NAN if we run off the image
-            float value = psImageInterpolatePixelBilinear(x, y, source->pixels);
-
-            // keep the nan values so all vectors are matched
-            // if (isnan(value)) continue;
-
+            dR = 2*(int)(0.5*r*sin(dtheta)) + 1;
+
+            if (x < 0) goto badvalue;
+            if (y < 0) goto badvalue;
+            if (x >= source->pixels->parent->numCols) goto badvalue;
+            if (y >= source->pixels->parent->numRows) goto badvalue;
+
+            float value = NAN;
+            if (dR < 2) {
+                // value is NAN if we run off the image
+                value = psImageInterpolatePixelBilinear(x, y, source->pixels);
+            } else {
+                value = psphotMeanSectorValue(source->pixels, x, y, dR, dR, theta);
+            }
+
+            // keep the all values (even NAN) so all vectors are matched in length
             psVectorAppend (radius, r);
             psVectorAppend (flux, value);
+            continue;
+            
+        badvalue:
+            psVectorAppend (radius, r);
+            psVectorAppend (flux, NAN);
         }
 
@@ -103 +118 @@
     return true;
 }
+
+float psphotMeanSectorValue (psImage *image, int x, int y, float dL, float dW, float theta) {
+
+    psVector *values = psphotBoxValues (image, x, y, dL, dW, theta);
+    
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+    psVectorStats (stats, values, NULL, NULL, 0);
+
+    float value = stats->sampleMedian;
+
+    psFree (stats);
+    psFree (values);
+    
+    return value;
+}
+
+psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta) {
+
+    // extract pixels from a series of lines (from -0.5*dW to +0.5*dW) of length dL, 
+    // centered on x0, y0 in parent pixel coordinates (not pixel indicies)
+
+    float xs = x0 - image->col0 - 0.5*dL*cos(theta);
+    float ys = y0 - image->row0 - 0.5*dL*sin(theta);
+
+    float xe = xs + 0.5*dL*cos(theta);
+    float ye = ys + 0.5*dL*sin(theta);
+
+    psVector *values = psphotLineValues (image, xs, ys, xe, ye, (int) dW);
+    return values;
+}
+
+/**
+ * identify the quadrant and draw the correct line
+ */
+psVector *psphotLineValues (psImage *image, double x1, double y1, double x2, double y2, int dW) {
+
+  int FlipDirect, FlipCoords;
+  int X1, Y1, X2, Y2, dX, dY;
+
+  /* rather than draw the line from float positions, we find the closest
+     integer end-points and draw the line between those pixels */
+
+  X1 = ROUND(x1);
+  Y1 = ROUND(y1);
+  X2 = ROUND(x2);
+  Y2 = ROUND(y2);
+
+  dX = X2 - X1;
+  dY = Y2 - Y1;
+
+  FlipCoords = (abs(dX) < abs(dY));
+  FlipDirect = FlipCoords ? (y1 > y2) : (x1 > x2);
+
+  psVector *values = NULL;
+  if (!FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X1, Y1, X2, Y2, dW, FALSE);
+  if ( FlipDirect && !FlipCoords) values = psphotLineValuesBresen (image, X2, Y2, X1, Y1, dW, FALSE);
+  if (!FlipDirect &&  FlipCoords) values = psphotLineValuesBresen (image, Y1, X1, Y2, X2, dW, TRUE);
+  if ( FlipDirect &&  FlipCoords) values = psphotLineValuesBresen (image, Y2, X2, Y1, X1, dW, TRUE);
+
+  return values;
+}
+
+/**
+ * use the Bresenham line drawing technique
+ * integer-only Bresenham line-draw version which is fast
+ */
+psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords) {
+
+    int X, Y, dX, dY;
+    int e, e2;
+
+    psVector *values = psVectorAllocEmpty(100, PS_TYPE_F32);
+
+    dX = X2 - X1;
+    dY = Y2 - Y1;
+
+    Y = Y1;
+    e = 0;
+    for (X = X1; X <= X2; X++) {
+        if (X > 0) {
+            if (swapcoords) {
+                if (X >= image->numRows) continue;
+                for (int y = Y - dW; y <= Y + dW; y++) {
+                    if (y < 0) continue;
+                    if (y >= image->numCols) continue;
+                    psVectorAppend(values, image->data.F32[X][y]);
+                }
+            } else {
+                if (X >= image->numCols) continue;
+                for (int y = Y - dW; y <= Y + dW; y++) {
+                    if (y < 0) continue;
+                    if (y >= image->numRows) continue;
+                    psVectorAppend(values, image->data.F32[y][X]);
+                }
+            }
+        }
+        e += dY;
+        e2 = 2 * e;
+        if (e2 > dX) {
+            Y++;
+            e -= dX;
+        }
+        if (e2 < -dX) {
+            Y--;
+            e += dX;
+        }
+    }
+    return values;
+}
+