Classification Of Point In Polygon

ptInPolygon - Point classification

This is an algorithm to test if a point is inside, outside or on (according to an input tolerance) a closed polygon definite with an ordinate set of points.
The classical algorithm is "Ray casting algorithm" (http://en.wikipedia.org/wiki/Point_in_polygon) and counts the number of intersections between the polygon and a segment that connect the point to a point that is out for sure. It has a lot of drawbacks and special cases.
The algorithm I propose counts how many times the polygon cross from one quadrant to another. It can be considered a variant of the "Winding number algorithm", described in the same article on wikipedia, but does not uses trigonometric functions.
It is more stable than the standard one and very fast, at least as fast as the Ray casting algorithm.
All the code is in a template class:

   template < class Pt, class PtIterator = Pt* > class ptInPolygon;

that must be used in as follow:

#include "ptInPolygon.h"
....
Pt2d ptest;
double tol = 0.1;
ptInPolygon < Pt2d, ptsarr::iterator > xx(tol);
GeoStuff::Locate r = xx(ptest, pts.begin(), pts.end()); // operator()
if(r == GeoStuff::In) {..}
else if(r == GeoStuff::Out) {..}
else if(r == GeoStuff::On) {..}
else {/*impossible*/}


The attached example uses this algorithm to classify a grid of points with a polygon. In the image you can see red out points, green inside points and blue points found on polygon according to the input tolerance.

All the code is in ptInPolygon.h, that I copy here, if you do not want to download the example.

// ==== Start =====

// ptInPolygon.h
#pragma once

#include "math.h"

// GeoMetric Stuff (Vect, tol, and other definitions).

class GeoStuff {
public:
GeoStuff(double tol) : m_tol2(sq(tol)){}
enum Locate{In=0, On, Out};
struct Vect {
Vect(double x=0, double y=0) { m_v[0] = x; m_v[1] = y;}
double m_v[2];
double & operator[](int i) {return m_v[i];}
double operator[](int i) const {return m_v[i];}
};
static inline double dot(Vect& a, Vect & b) {
return a[0]*b[0]+a[1]*b[1];
}
static inline double cross(Vect& a, Vect & b) {
return a[0]*b[1]-a[1]*b[0];
}
static inline double sq(double a) {return a*a;}

template < class Pt >
static Vect diff(const Pt& a, const Pt & b) {
return Vect(a[0]-b[0],a[1]-b[1]);
}
template < class Pt >
bool equiv(const Pt & a, const Pt & b){
return sq(a[0]-b[0])+sq(a[1]-b[1]) <= m_tol2;
}
protected:
double m_tol2;
static char quadrante(Vect & v) {
static char qq[2][2]={{0,1},{3,2}};
return qq[v[1] < 0][v[0] < 0];
}
};

// The class ptInPolygon is able to calssify points as inside, outside or on a polygon
// defined by a set of successve points.
// The point class must only have operator []
// to access x (coordinate 0) and y (coordinate 1).

template < class Pt, class PtIterator = Pt* >
class ptInPolygon : public GeoStuff {
private:
public:
ptInPolygon(double tol): GeoStuff(tol){}
Locate operator()(const Pt & pt, const PtIterator & begin, const PtIterator & end)
{
int npts = 0;
for(PtIterator it = begin; it != end; ++it)
{
if(equiv((*it), pt))
return On;
++npts;
}
if(npts < 2)
return Out;
Vect vp=diff(*begin,pt);
Vect va;
int qp = quadrante(vp), qa;
int quad_crosses = 0;
for(PtIterator it = end-1; ;--it) {
va = diff((*it), pt);
qa = quadrante(va);
int diff = qa - qp;
if(diff != 0) {
double vv = dot(vp,vp), ww = dot(va,va), vw = dot(va,vp);
double den = vv+ww-vw*2;
if((vv*ww-sq(vw)) <= m_tol2 * den) {// mindist < tol
if(fabs(ww-vv) < den) // projection inside segment
return On;
}
switch(diff) {
case -3:
case 1:
++quad_crosses;
break;
case -1:
case 3:
--quad_crosses;
break;
case -2:
case 2:{
if(cross(vp,va) < 0)
quad_crosses -= 2;
else
quad_crosses += 2;
break;
}
}
}
if(it == begin)
break;
vp = va;
qp = qa;
}
return abs(quad_crosses) > 2 ? In : Out;
}
};

// ==== End =====



About the Author

Marco Amagliani

Geometric kernel developer for a CAD industry

Downloads

Comments

  • Polygon triangulation into an uniform rectangular grid

    Posted by rraallvv on 06/05/2012 12:24pm

    I've searched the entire web for something like this, fortunately you got the code, thanks a lot. I also wanted to use it to triangulate a polygon into a uniform rectangular grid, any suggestion on this one. Regards.

    Reply
Leave a Comment
  • Your email address will not be published. All fields are required.

Top White Papers and Webcasts

  • Live Event Date: August 20, 2014 @ 1:00 p.m. ET / 10:00 a.m. PT When you look at natural user interfaces as a developer, it isn't just fun and games. There are some very serious, real-world usage models of how things can help make the world a better place – things like Intel® RealSense™ technology. Check out this upcoming eSeminar and join the panel of experts, both from inside and outside of Intel, as they discuss how natural user interfaces will likely be getting adopted in a wide variety …

  • Managing your company's financials is the backbone of your business and is vital to the long-term health and viability of your company. To continue applying the necessary financial rigor to support rapid growth, the accounting department needs the right tools to most efficiently do their job. Read this white paper to understand the 10 essentials of a complete financial management system and how the right solution can help you keep up with the rapidly changing business world.

Most Popular Programming Stories

More for Developers

Latest Developer Headlines

RSS Feeds