|
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 |
| -#defineHIT_IT INT_MAX |
| 5076 | +#definePOINT_ON_POLYGON INT_MAX |
5075 | 5077 |
|
5076 | 5078 | staticint
|
5077 | 5079 | point_inside(Point*p,intnpts,Point*plist)
|
5078 | 5080 | {
|
5079 | 5081 | doublex0,
|
5080 | 5082 | y0;
|
5081 |
| -doublepx, |
5082 |
| -py; |
5083 |
| -inti; |
| 5083 | +doubleprev_x, |
| 5084 | +prev_y; |
| 5085 | +inti=0; |
5084 | 5086 | doublex,
|
5085 | 5087 | y;
|
5086 |
| -intcross, |
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 | +intcross,total_cross=0; |
5097 | 5089 |
|
5098 |
| -crossnum=0; |
5099 |
| -i=0; |
5100 | 5090 | if (npts <=0)
|
5101 | 5091 | return0;
|
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 | return2;
|
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 | return2;
|
5122 |
| -crossnum+=cross; |
5123 |
| -if (crossnum!=0) |
| 5118 | +total_cross+=cross; |
| 5119 | + |
| 5120 | +if (total_cross!=0) |
5124 | 5121 | return1;
|
5125 | 5122 | return0;
|
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 | staticint
|
5138 |
| -lseg_crossing(doublex,doubley,doublepx,doublepy) |
| 5136 | +lseg_crossing(doublex,doubley,doubleprev_x,doubleprev_y) |
5139 | 5137 | {
|
5140 | 5138 | doublez;
|
5141 |
| -intsgn; |
5142 |
| - |
5143 |
| -/* If (px,py) = (0,0) and not first call we have already sent HIT_IT */ |
| 5139 | +inty_sign; |
5144 | 5140 |
|
5145 | 5141 | if (FPzero(y))
|
5146 |
| -{ |
5147 |
| -if (FPzero(x)) |
5148 |
| -{ |
5149 |
| -returnHIT_IT; |
5150 |
| - |
5151 |
| -} |
| 5142 | +{/* y == 0, on X axis */ |
| 5143 | +if (FPzero(x))/* (x,y) is (0,0)? */ |
| 5144 | +returnPOINT_ON_POLYGON; |
5152 | 5145 | elseif (FPgt(x,0))
|
5153 |
| -{ |
5154 |
| -if (FPzero(py)) |
5155 |
| -returnFPgt(px,0) ?0 :HIT_IT; |
5156 |
| -returnFPlt(py,0) ?1 :-1; |
5157 |
| - |
| 5146 | +{/* x > 0 */ |
| 5147 | +if (FPzero(prev_y))/* y and prev_y are zero */ |
| 5148 | +/* prev_x > 0? */ |
| 5149 | +returnFPgt(prev_x,0) ?0 :POINT_ON_POLYGON; |
| 5150 | +returnFPlt(prev_y,0) ?1 :-1; |
5158 | 5151 | }
|
5159 | 5152 | else
|
5160 |
| -{/* x < 0 */ |
5161 |
| -if (FPzero(py)) |
5162 |
| -returnFPlt(px,0) ?0 :HIT_IT; |
| 5153 | +{/* x < 0, x not on positive X axis */ |
| 5154 | +if (FPzero(prev_y)) |
| 5155 | +/* prev_x < 0? */ |
| 5156 | +returnFPlt(prev_x,0) ?0 :POINT_ON_POLYGON; |
5163 | 5157 | return0;
|
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 |
| -returnFPlt(px,0) ?0 :sgn; |
5171 |
| - |
5172 |
| -if (FPgt((sgn*py),0)) |
5173 |
| -{/* y and py have same sign */ |
5174 |
| -return0; |
5175 |
| - |
5176 |
| -} |
5177 | 5160 | else
|
5178 |
| -{/* y and py have opposite signs */ |
5179 |
| -if (FPge(x,0)&&FPgt(px,0)) |
5180 |
| -return2*sgn; |
5181 |
| -if (FPlt(x,0)&&FPle(px,0)) |
5182 |
| -return0; |
5183 |
| - |
5184 |
| -z= (x-px)*y- (y-py)*x; |
5185 |
| -if (FPzero(z)) |
5186 |
| -returnHIT_IT; |
5187 |
| -returnFPgt((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 | +returnFPlt(prev_x,0) ?0 :y_sign; |
| 5168 | +elseif (FPgt(y_sign*prev_y,0)) |
| 5169 | +/* both above or below X axis */ |
| 5170 | +return0;/* 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 | +return2*y_sign; |
| 5176 | +if (FPlt(x,0)&&FPle(prev_x,0)) |
| 5177 | +/* both non-positive so do not cross positive X-axis */ |
| 5178 | +return0; |
| 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 | +returnPOINT_ON_POLYGON; |
| 5184 | +returnFPgt((y_sign*z),0) ?0 :2*y_sign; |
| 5185 | +} |
5188 | 5186 | }
|
5189 | 5187 | }
|
5190 | 5188 |
|
|