点在多边形内的判断

geos-3.0.3中的代码，扫描线法:

 

/*public static*/boolCGAlgorithms::isPointInRing(const Coordinate& p, const Coordinate::ConstVect& ring){ double xInt; // x intersection of segment with ray int crossings = 0; // number of segment/ray crossings double x1; // translated coordinates double y1; double x2; double y2; /* * For each segment l = (i-1, i), see if it crosses ray from * test point in positive x direction. */ for(size_t i=1, nPts=ring.size(); i<nPts; ++i) { const Coordinate *p1=ring[i]; const Coordinate *p2=ring[i-1]; x1 = p1->x - p.x; y1 = p1->y - p.y; x2 = p2->x - p.x; y2 = p2->y - p.y; if (((y1 > 0) && (y2 <= 0)) || ((y2 > 0) && (y1 <= 0))) { /* * segment straddles x axis, so compute intersection. */ xInt = RobustDeterminant::signOfDet2x2(x1, y1, x2, y2) / (y2 - y1); /* * crosses ray if strictly positive intersection. */ if (0.0 < xInt) crossings++; } } /* * p is inside if number of crossings is odd. */ if ((crossings % 2) == 1) return true; return false;}

 

附:

int RobustDeterminant::signOfDet2x2(double x1,double y1,double x2,double y2) {// returns -1 if the determinant is negative,// returns 1 if the determinant is positive,// retunrs 0 if the determinant is null.int sign=1;double swap;double k;long count=0;/** testing null entries*/if ((x1==0.0) || (y2==0.0)) {if ((y1==0.0) || (x2==0.0)) {return 0;} else if (y1>0) {if (x2>0) {return -sign;} else {return sign;}} else {if (x2>0) {return sign;} else {return -sign;}}}if ((y1==0.0) || (x2==0.0)) {if (y2>0) {if (x1>0) {return sign;} else {return -sign;}} else {if (x1>0) {return -sign;} else {return sign;}}}/** making y coordinates positive and permuting the entries* so that y2 is the biggest one*/if (0.0<y1) {if (0.0<y2) {if (y1<=y2) {;} else {sign=-sign;swap=x1;x1=x2;x2=swap;swap=y1;y1=y2;y2=swap;}} else {if (y1<=-y2) {sign=-sign;x2=-x2;y2=-y2;} else {swap=x1;x1=-x2;x2=swap;swap=y1;y1=-y2;y2=swap;}}} else {if (0.0<y2) {if (-y1<=y2) {sign=-sign;x1=-x1;y1=-y1;} else {swap=-x1;x1=x2;x2=swap;swap=-y1;y1=y2;y2=swap;}} else {if (y1>=y2) {x1=-x1;y1=-y1;x2=-x2;y2=-y2;} else {sign=-sign;swap=-x1;x1=-x2;x2=swap;swap=-y1;y1=-y2;y2=swap;}}}/** making x coordinates positive*//** if |x2|<|x1| one can conclude*/if (0.0<x1) {if (0.0<x2) {if (x1 <= x2) {;} else {return sign;}} else {return sign;}} else {if (0.0<x2) {return -sign;} else {if (x1 >= x2) {sign=-sign;x1=-x1;x2=-x2;} else {return -sign;}}}/** all entries strictly positive x1 <= x2 and y1 <= y2*/while (true) {count=count+1;k=floor(x2/x1);x2=x2-k*x1;y2=y2-k*y1;/** testing if R (new U2) is in U1 rectangle*/if (y2<0.0) {return -sign;}if (y2>y1) {return sign;}/** finding R'*/if (x1>x2+x2) {if (y1<y2+y2) {return sign;}} else {if (y1>y2+y2) {return -sign;} else {x2=x1-x2;y2=y1-y2;sign=-sign;}}if (y2==0.0) {if (x2==0.0) {return 0;} else {return -sign;}}if (x2==0.0) {return sign;}/** exchange 1 and 2 role.*/k=floor(x1/x2);x1=x1-k*x2;y1=y1-k*y2;/** testing if R (new U1) is in U2 rectangle*/if (y1<0.0) {return sign;}if (y1>y2) {return -sign;}/** finding R'*/if (x2>x1+x1) {if (y2<y1+y1) {return -sign;}} else {if (y2>y1+y1) {return sign;} else {x1=x2-x1;y1=y2-y1;sign=-sign;}}if (y1==0.0) {if (x1==0.0) {return 0;} else {return sign;}}if (x1==0.0) {return -sign;}}}

 

 

ogre在1.6中亦有实现，采用方法是计算该点与多边形边的角之和是否等于360度。

//-----------------------------------------------------------------------bool Polygon::isPointInside(const Vector3& point) const{// sum the angles Real anglesum = 0;size_t n = getVertexCount();for (size_t i = 0; i < n; i++) {const Vector3& p1 = getVertex(i);const Vector3& p2 = getVertex((i + 1) % n);Vector3 v1 = p1 - point;Vector3 v2 = p2 - point;Real len1 = v1.length();Real len2 = v2.length();if (Math::RealEqual(len1 * len2, 0.0f, 1e-4f)){// We are on a vertex so consider this insidereturn true; }else{Real costheta = v1.dotProduct(v2) / (len1 * len2);anglesum += acos(costheta);}}// result should be 2*PI if point is inside polyreturn Math::RealEqual(anglesum, Math::TWO_PI, 1e-4f);}