Computational Geometry primer
In this series of tutorials we will try to cover some basic mathematic algorithms needed in computer graphics and games. No need to be afraid if math is not your strong point. We will not delve with high level calculus and any scary science. We will instead use some high school math and geometry in order to explain the algorithms and then the code we will develop will be ready to use as is. Just copy and paste into your projects and you are done. Here we go then.
Go to part 1, Line segment intersection
2. Point in polygon
This article addresses the point in polygon problem. This problem arises many times in computer programs that have to deal with geometric problems. The situation arises quite often in games. Suppose we have two people chasing each other as in the diagram. At some time the blue player decides to shoot at the red, it's up to the game engine to tell if the shot has any chance of succeeding. Or even in a puzzle game, you are asked to determine if the red dot can move out of the enclosed area. Knowing that the area is closed (easy to tell: the last and the first points are the same), if the point in INSIDE the polygon there is no way out. I bet you have many examples of your own when a point in polygon algorithm is needed.
The algorithm we are going to develop here relies on the line segment intersection developed in the previous tutorial.
Determining weather the red point lies inside the polygon is not so hard. Just imagine a horizontal line from that point to the right long enough to go past the right edge of a bounding rectangle around the polygon. Notice the number of intersections that this segment has with the polygon. If the point is inside the number of intersections is odd, otherwise it is even. There is one situation we have to take into account, which might lead us to error results. That is what do we do if the intersection point is one of the points determining the polygon.
The biggest problem though is that we might count the intersection twice, once for each edge. So we must ignore either intersections at the starting point of an edge or at the end point. When the intersection point is a point of the polygon there are two cases. First the two edges cross the horizontal line, or they are on the same side. In the first case we have an intersection, but in the later we don't. The solution is simple. Assume that the intersection point is on segment (p1,p2) and on point p2. The next segment is (p2,p3). In order to check if p1 and p3 lie on either side of the horizontal line at p2 we check
(p2.y-p1.y) * (p3.y-p2.y)
If the result is greater than 0 the boundary crosses the line and the intersection is valid. Otherwise the three points are on the same side and the intersection is invalid.
For the need of this tutorial I created a simple polygon class that calculates its boundaries as points are added. The code uses some simple trick I learned from the older gurus to speed up the calculations by avoiding time consuming calculations when possible.
class clf_point2D{
public:
float x;
float y;
clf_point2D():x(0.f), y(0.f){
}
clf_point2D(float _x, float _y):x(_x), y(_y){
}
clf_point2D(int _x, int _y):x((float)_x), y((float)_y){
}
clf_point2D(clf_point2D& p) : x(p.x), y(p.y){
}
~clf_point2D(){
}
};
class clf_rect2D{
public:
float left;
float top;
float right;
float bottom;
clf_rect2D() : left(0.f),top(0.f),right(0.f),bottom(0.f) {
}
clf_rect2D(float l, float t, float r, float b) : left(l),top(t),right(r),bottom(b) {
}
clf_rect2D(clf_rect2D& p) : left(p.left),top(p.top),right(p.right),bottom(p.bottom) {
}
~clf_rect2D(){
}
};
class polygon : public std::vector< clf_point2D>{
public:
clf_rect2D bounds;
polygon();
~polygon();
void add_point(const CPoint& pt);
void draw(CDC* pDC);
};
polygon::polygon()
{
}
polygon::~polygon()
{
}
void polygon::add_point(const CPoint& pt)
{
if (size() == 0)
{
bounds = clf_rect2D((float)pt.x, (float)pt.y, (float)pt.x, (float)pt.y);
}
else
{
if (bounds.left > (float)pt.x) bounds.left = (float)pt.x;
if (bounds.right < (float)pt.x) bounds.right = (float)pt.x;
if (bounds.top < (float)pt.y) bounds.top = (float)pt.y;
if (bounds.bottom > (float)pt.y) bounds.bottom = (float)pt.y;
}
push_back(clf_point2D(pt.x, pt.y));
}
void polygon::draw(CDC* pDC)
{
if (!size()) return;
std::vector<clf_point2D>::iterator it;
it=begin();
pDC->MoveTo((int)it->x, (int)it->y);
for (; it!=end(); ++it)
{
pDC->LineTo((int)it->x, (int)it->y);
}
it=begin();
pDC->LineTo((int)it->x, (int)it->y);
}
bool intersect(const clf_point2D& p1, const clf_point2D& p2, const clf_point2D& p3, const clf_point2D& p4, clf_point2D* p0)
{
// first check horizontal and vertical gap
if (min(p1.x, p2.x) > max(p3.x,p4.x))
return false;
if (max(p1.x, p2.x) < min(p3.x,p4.x))
return false;
if (min(p1.y, p2.y) > max(p3.y,p4.y))
return false;
if (max(p1.y, p2.y) < min(p3.y,p4.y))
return false;
// firsty calculate the common denominator
float d = (p4.y-p3.y)*(p2.x-p1.x)-(p4.x-p3.x)*(p2.y-p1.y);
if (fabs(d) < 0.0000001) // almost 0 (compensate for computer accuracy)
return false;
// calculate fractional part for first segment
float na = ((p4.x-p3.x)*(p1.y-p3.y)-(p4.y-p3.y)*(p1.x-p3.x))/d;
if (na < 0.f || na > 1.f) // beyond limits of first segment, search no more
return false;
// calculate fractional part for second segment
float nb = ((p2.x-p1.x)*(p1.y-p3.y)-(p2.y-p1.y)*(p1.x-p3.x))/d;
if (nb < 0.f || nb > 1.f) // beyond limits of second segment, search no more
return false;
// we have intersection, calculate actual coordinates (if required)
if (p0)
{
p0->x = p1.x+na*(p2.x-p1.x);
p0->y = p1.y+na*(p2.y-p1.y);
// or using the second segment endpoints
// p0->x = p3.x+nb*(p4.x-p3.x);
// p0->y = p3.y+nb*(p4.y-p3.y);
}
// return true
return true;
}
bool point_in_polygon(polygon& poly, clf_point2D pt)
{
// first check if the point is out of the bounds of the polygon
if (pt.x < poly.bounds.left || pt.x > poly.bounds.right || pt.y < poly.bounds.bottom || pt.y > poly.bounds.top)
return false;
// generate a point out of the polygon with the same y coordinate as the check point
clf_point2D right(poly.bounds.right+1.f, pt.y);
clf_point2D isect;
// now iterate through the points of the polygon
polygon::iterator next = poly.begin();
polygon::iterator it;
int count = 0;
for (it=poly.begin() ; it!=poly.end(); ++it)
{
// -next- points to the 'next' point while -it- is the current
++next;
if (next == poly.end())
next = poly.begin();
if (intersect(*it, *next, pt, right, &isect))
{
if (fabs(it->y-isect.y) < 0.001) // the benefit of horizontal lines is that we only need to check the y coordinate
{
// ignore it, we will catch it as ending point later
}
// if the ending point of the segment is the intersection
else if (fabs(next->y-isect.y) < 0.001) // the benefit of horizontal lines is that we only need to check the y coordinate
{
polygon::iterator next_next = next;
++next_next;
if (next_next == poly.end())
next_next = poly.begin();
// check if the two segments are on the same side of the horizontal line
if ((next_next->y-next->y)*(next->y-it->y) > 0) // dy1*dy2>0 ==> the polygon intersects
{
count++; // boundary crossing, count the intersection
}
else
{
// not an actual intersection, do not count it
}
}
else
{ // intersection point, count it
count ++;
}
}
}
if (count%2) // odd number of intersections
return true;
return false;
}