|
8 | 8 | *
|
9 | 9 | *
|
10 | 10 | * IDENTIFICATION
|
11 |
| - * $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.95 2007/02/27 23:48:08 tgl Exp $ |
| 11 | + * $PostgreSQL: pgsql/src/backend/utils/adt/geo_ops.c,v 1.96 2007/03/05 23:29:14 momjian Exp $ |
12 | 12 | *
|
13 | 13 | *-------------------------------------------------------------------------
|
14 | 14 | */
|
@@ -5063,128 +5063,126 @@ poly_circle(PG_FUNCTION_ARGS)
|
5063 | 5063 | ***********************************************************************/
|
5064 | 5064 |
|
5065 | 5065 | /*
|
5066 |
| - * Test to see if the point is inside the polygon. |
5067 |
| - * Code adapted from integer-based routines in WN: A Server for the HTTP |
| 5066 | + * Test to see if the point is inside the polygon, returns 1/0, or 2 if |
| 5067 | + * the point is on the polygon. |
| 5068 | + * Code adapted but not copied from integer-based routines in WN: A |
| 5069 | + * Server for the HTTP |
5068 | 5070 | * version 1.15.1, file wn/image.c
|
5069 |
| - * GPL Copyright (C) 1995 by John Franks |
5070 | 5071 | * http://hopf.math.northwestern.edu/index.html
|
5071 | 5072 | * Description of algorithm: http://www.linuxjournal.com/article/2197
|
| 5073 | + * http://www.linuxjournal.com/article/2029 |
5072 | 5074 | */
|
5073 | 5075 |
|
5074 |
| -#define HIT_IT INT_MAX |
| 5076 | +#define POINT_ON_POLYGON INT_MAX |
5075 | 5077 |
|
5076 | 5078 | static int
|
5077 | 5079 | point_inside(Point *p, int npts, Point *plist)
|
5078 | 5080 | {
|
5079 | 5081 | double x0,
|
5080 | 5082 | y0;
|
5081 |
| - double px, |
5082 |
| - py; |
5083 |
| - int i; |
| 5083 | + double prev_x, |
| 5084 | + prev_y; |
| 5085 | + int i = 0; |
5084 | 5086 | double x,
|
5085 | 5087 | y;
|
5086 |
| - int cross, |
5087 |
| - crossnum; |
5088 |
| - |
5089 |
| -/* |
5090 |
| - * We calculate crossnum, which is twice the crossing number of a |
5091 |
| - * ray from the origin parallel to the positive X axis. |
5092 |
| - * A coordinate change is made to move the test point to the origin. |
5093 |
| - * Then the function lseg_crossing() is called to calculate the crossnum of |
5094 |
| - * one segment of the translated polygon with the ray which is the |
5095 |
| - * positive X-axis. |
5096 |
| - */ |
| 5088 | + int cross, total_cross = 0; |
5097 | 5089 |
|
5098 |
| - crossnum = 0; |
5099 |
| - i = 0; |
5100 | 5090 | if (npts <= 0)
|
5101 | 5091 | return 0;
|
5102 | 5092 |
|
| 5093 | + /* compute first polygon point relative to single point */ |
5103 | 5094 | x0 = plist[0].x - p->x;
|
5104 | 5095 | y0 = plist[0].y - p->y;
|
5105 | 5096 |
|
5106 |
| - px = x0; |
5107 |
| - py = y0; |
| 5097 | + prev_x = x0; |
| 5098 | + prev_y = y0; |
| 5099 | + /* loop over polygon points and aggregate total_cross */ |
5108 | 5100 | for (i = 1; i < npts; i++)
|
5109 | 5101 | {
|
| 5102 | + /* compute next polygon point relative to single point */ |
5110 | 5103 | x = plist[i].x - p->x;
|
5111 | 5104 | y = plist[i].y - p->y;
|
5112 | 5105 |
|
5113 |
| - if ((cross = lseg_crossing(x, y, px, py)) == HIT_IT) |
| 5106 | + /* compute previous to current point crossing */ |
| 5107 | + if ((cross = lseg_crossing(x, y, prev_x, prev_y)) == POINT_ON_POLYGON) |
5114 | 5108 | return 2;
|
5115 |
| - crossnum += cross; |
5116 |
| - |
5117 |
| - px = x; |
5118 |
| - py = y; |
| 5109 | + total_cross += cross; |
| 5110 | + |
| 5111 | + prev_x = x; |
| 5112 | + prev_y = y; |
5119 | 5113 | }
|
5120 |
| - if ((cross = lseg_crossing(x0, y0, px, py)) == HIT_IT) |
| 5114 | + |
| 5115 | + /* now do the first point */ |
| 5116 | + if ((cross = lseg_crossing(x0, y0, prev_x, prev_y)) == POINT_ON_POLYGON) |
5121 | 5117 | return 2;
|
5122 |
| - crossnum += cross; |
5123 |
| - if (crossnum != 0) |
| 5118 | + total_cross += cross; |
| 5119 | + |
| 5120 | + if (total_cross != 0) |
5124 | 5121 | return 1;
|
5125 | 5122 | return 0;
|
5126 | 5123 | }
|
5127 | 5124 |
|
5128 | 5125 |
|
5129 | 5126 | /* lseg_crossing()
|
5130 |
| - * The function lseg_crossing() returns +2, or -2 if the segment from (x,y) |
5131 |
| - * to previous (x,y) crosses the positive X-axis positively or negatively. |
5132 |
| - * It returns +1 or -1 if one endpoint is on this ray, or 0 if both are. |
5133 |
| - * It returns 0 if the ray and the segment don't intersect. |
5134 |
| - * It returns HIT_IT if the segment contains (0,0) |
| 5127 | + * Returns +/-2 if line segment crosses the positive X-axis in a +/- direction. |
| 5128 | + * Returns +/-1 if one point is on the positive X-axis. |
| 5129 | + * Returns 0 if both points are on the positive X-axis, or there is no crossing. |
| 5130 | + * Returns POINT_ON_POLYGON if the segment contains (0,0). |
| 5131 | + * Wow, that is one confusing API, but it is used above, and when summed, |
| 5132 | + * can tell is if a point is in a polygon. |
5135 | 5133 | */
|
5136 | 5134 |
|
5137 | 5135 | static int
|
5138 |
| -lseg_crossing(double x, double y, double px, double py) |
| 5136 | +lseg_crossing(double x, double y, double prev_x, double prev_y) |
5139 | 5137 | {
|
5140 | 5138 | double z;
|
5141 |
| - int sgn; |
5142 |
| - |
5143 |
| - /* If (px,py) = (0,0) and not first call we have already sent HIT_IT */ |
| 5139 | + int y_sign; |
5144 | 5140 |
|
5145 | 5141 | if (FPzero(y))
|
5146 |
| - { |
5147 |
| - if (FPzero(x)) |
5148 |
| - { |
5149 |
| - return HIT_IT; |
5150 |
| - |
5151 |
| - } |
| 5142 | + { /* y == 0, on X axis */ |
| 5143 | + if (FPzero(x)) /* (x,y) is (0,0)? */ |
| 5144 | + return POINT_ON_POLYGON; |
5152 | 5145 | else if (FPgt(x, 0))
|
5153 |
| - { |
5154 |
| - if (FPzero(py)) |
5155 |
| - return FPgt(px, 0) ? 0 : HIT_IT; |
5156 |
| - return FPlt(py, 0) ? 1 : -1; |
5157 |
| - |
| 5146 | + { /* x > 0 */ |
| 5147 | + if (FPzero(prev_y)) /* y and prev_y are zero */ |
| 5148 | + /* prev_x > 0? */ |
| 5149 | + return FPgt(prev_x, 0) ? 0 : POINT_ON_POLYGON; |
| 5150 | + return FPlt(prev_y, 0) ? 1 : -1; |
5158 | 5151 | }
|
5159 | 5152 | else
|
5160 |
| - { /* x < 0 */ |
5161 |
| - if (FPzero(py)) |
5162 |
| - return FPlt(px, 0) ? 0 : HIT_IT; |
| 5153 | + { /* x < 0, x not on positive X axis */ |
| 5154 | + if (FPzero(prev_y)) |
| 5155 | + /* prev_x < 0? */ |
| 5156 | + return FPlt(prev_x, 0) ? 0 : POINT_ON_POLYGON; |
5163 | 5157 | return 0;
|
5164 | 5158 | }
|
5165 | 5159 | }
|
5166 |
| - |
5167 |
| - /* Now we know y != 0; set sgn to sign of y */ |
5168 |
| - sgn = (FPgt(y, 0) ? 1 : -1); |
5169 |
| - if (FPzero(py)) |
5170 |
| - return FPlt(px, 0) ? 0 : sgn; |
5171 |
| - |
5172 |
| - if (FPgt((sgn * py), 0)) |
5173 |
| - { /* y and py have same sign */ |
5174 |
| - return 0; |
5175 |
| - |
5176 |
| - } |
5177 | 5160 | else
|
5178 |
| - { /* y and py have opposite signs */ |
5179 |
| - if (FPge(x, 0) && FPgt(px, 0)) |
5180 |
| - return 2 * sgn; |
5181 |
| - if (FPlt(x, 0) && FPle(px, 0)) |
5182 |
| - return 0; |
5183 |
| - |
5184 |
| - z = (x - px) * y - (y - py) * x; |
5185 |
| - if (FPzero(z)) |
5186 |
| - return HIT_IT; |
5187 |
| - return FPgt((sgn * z), 0) ? 0 : 2 * sgn; |
| 5161 | + { /* y != 0 */ |
| 5162 | + /* compute y crossing direction from previous point */ |
| 5163 | + y_sign = FPgt(y, 0) ? 1 : -1; |
| 5164 | + |
| 5165 | + if (FPzero(prev_y)) |
| 5166 | + /* previous point was on X axis, so new point is either off or on */ |
| 5167 | + return FPlt(prev_x, 0) ? 0 : y_sign; |
| 5168 | + else if (FPgt(y_sign * prev_y, 0)) |
| 5169 | + /* both above or below X axis */ |
| 5170 | + return 0; /* same sign */ |
| 5171 | + else |
| 5172 | + { /* y and prev_y cross X-axis */ |
| 5173 | + if (FPge(x, 0) && FPgt(prev_x, 0)) |
| 5174 | + /* both non-negative so cross positive X-axis */ |
| 5175 | + return 2 * y_sign; |
| 5176 | + if (FPlt(x, 0) && FPle(prev_x, 0)) |
| 5177 | + /* both non-positive so do not cross positive X-axis */ |
| 5178 | + return 0; |
| 5179 | + |
| 5180 | + /* x and y cross axises, see URL above point_inside() */ |
| 5181 | + z = (x - prev_x) * y - (y - prev_y) * x; |
| 5182 | + if (FPzero(z)) |
| 5183 | + return POINT_ON_POLYGON; |
| 5184 | + return FPgt((y_sign * z), 0) ? 0 : 2 * y_sign; |
| 5185 | + } |
5188 | 5186 | }
|
5189 | 5187 | }
|
5190 | 5188 |
|
|
0 commit comments