Friday, February 26, 2010

Determining if a point lies on the interior of a polygon


 

Solution 1 (2D)

Determining whether or not a point (x,y) lies inside or outside a 2D polygonally bounded plane. This is necessary for example in applications such as polygon filling on raster devices.

Consider a polygon made up of N vertices (xi,yi) where i ranges from 0 to N-1. The last vertex (xN,yN) is assumed to be the same as the first vertex (x0,y0), that is, the polygon is closed. To determine the status of a point (xp,yp) consider a horizontal ray emanating from (xp,yp) and to the right.

If the number of times this ray intersects the line segments making up the polygon is even then the point is outside the polygon.

Whereas if the number of intersections is odd then the point (xp,yp) lies inside the polygon.

The following shows the ray for some sample points and should make the technique clear.


0 will be considered even, the test for even or odd will be based on modulus 2, that is, if the number of intersections modulus 2 is 0 then the number is even, if it is 1 then it is odd.

The only trick is what happens in the special cases when an edge or vertex of the polygon lies on the ray from (xp,yp). The possible situations are illustrated below.


The thick lines above are not considered as valid intersections, the thin lines do count as intersections. Ignoring the case of an edge lying along the ray or an edge ending on the ray ensures that the endpoints are only counted once.

Note that this algorithm also works for polygons with holes as illustrated below


The following C function (Taken from website)

returns INSIDE or OUTSIDE indicating the status of a point P with respect to a polygon with N points.

#define MIN(x,y) (x < y ? x : y)

#define MAX(x,y) (x > y ? x : y)

#define INSIDE 0

#define OUTSIDE 1


 

typedef struct {

double x,y;

} Point;


 

int InsidePolygon(Point *polygon,int N,Point p)

{

int counter = 0;

int i;

double xinters;

Point p1,p2;


 

p1 = polygon[0];

for (i=1;i<=N;i++) {

p2 = polygon[i % N];

if (p.y > MIN(p1.y,p2.y)) {

if (p.y <= MAX(p1.y,p2.y)) {

if (p.x <= MAX(p1.x,p2.x)) {

if (p1.y != p2.y) {

xinters = (p.y-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x;

if (p1.x == p2.x || p.x <= xinters)

counter++;

}

}

}

}

p1 = p2;

}


 

if (counter % 2 == 0)

return(OUTSIDE);

else

return(INSIDE);

}


 

Another : It returns 1 for interior points and 0 for exterior points.

int pnpoly(int npol, float *xp, float *yp, float x, float y)

{

int i, j, c = 0;

for (i = 0, j = npol-1; i < npol; j = i++) {

if ((((yp[i] <= y) && (y < yp[j])) ||

((yp[j] <= y) && (y < yp[i]))) &&

(x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))

c = !c;

}

return c;

}

No comments:

Post a Comment