Solved
convex hull and point in convex hull algorithms
Posted on 2004-09-07
given a vector<POINT> of points and another POINT i need to generate the convex hull of points and then determine point is inside or on the line of the convex hull. i have tried to implement grahams scan and then use the pnpoly algorithm from the comp.graphics faq but something is wrong. after finding the convex hull i use pnpoly. if pnpoly returns 1 the points is in the hull else if it is 0 it is possibly on the border so i go through every line possible line segment and check to see if the point is on the line. if it is i return true. i would prefer an algorithm that gives a definite answer to if the point is outside, inside, or on the border. for the ConvexHull algorithm i have attempted to implement Graham's Scan from CLR but I dont think I have done it correctly. I think my sorting may be off. i use sort() with
bool AngleCompare( POINT p1, POINT p2 )
{
p1.x -= p0ref.x;
p1.y -= p0ref.y;
p2.x -= p0ref.x;
p2.y -= p0ref.y;
//return CrossProd( p1, p2 ) > 0;
return atan2( (double) p1.y, p1.x ) < atan2( (double) p2.y, p2.x );
}
the algorithm says to sort the points in counterclockwise order around p0ref. i didnt know how to do this exactly so i searched online and round that it could be done based on the crossproduct but that didnt seem to work so i switched to using atan2 after translating p0ref ot the origin.
here is the rest of the code
bool Equal(double x, double y, double epsilon)
{
//PRECONDITION(epsilon >= 0.0);
// Get the exponent of the smaller of the two numbers (using the
// smaller of the two gives us a tighter epsilon value).
int exponent;
(void) frexp(fabs(x) < fabs(y) ? x : y, &exponent);
// Scale epsilon by the exponent.
double delta = ldexp(epsilon, exponent);
// If the difference between the numbers is smaller than the
// scaled epsilon we'll consider the numbers to be equal.
double difference = fabs(x - y);
bool equal = difference <= delta;
return equal;
}
int pnpoly(vector<POINT> vPoints, POINT pt)
{ int i, j, c = 0;
for (i = 0, j = vPoints.size()-1; i < vPoints.size(); j = i++) {
if ((((vPoints[i].y<=pt.y) && (pt.y<vPoints[j].y)) || ((vPoints[j].y<=pt.y)
&& (pt.y<vPoints[i].y))) &&
(pt.x < (vPoints[j].x - vPoints[i].x) * (pt.y - vPoints[i].y) /
(vPoints[j].y - vPoints[i].y) + vPoints[i].x))
c = !c; }
return c;
}
bool PointInConvexPoly( vector<POINT> vPoints, POINT ptThePoint )
{
int i;
double m, b, y;
bool online;
if( pnpoly( vPoints, ptThePoint ) )
return true;
else
{
vPoints.push_back( vPoints[0] );
online = false;
for( i = 0; i < vPoints.size() - 1; ++i )
{
m = (vPoints[i + 1].y - vPoints[i].y) / (double) (vPoints[i + 1].x - vPoints[i].x);
b = vPoints[i].y - m * vPoints[i].x;
y = m * ptThePoint.x + b;
if( abs( y - ptThePoint.y ) < 1E-6 )
online = true;
}
if( online )
return true;
else
return false;
}
}
bool ConvexHull ( vector<POINT> vPoints, POINT ptThePoint )
{
int size, i, minYIndex;
POINT ptMinY;
stack<POINT> theStack;
POINT p0, p1, p2;
size = (int) vPoints.size ();
for ( i = 0, minYIndex = 0; i < size; ++i )
{
if ( vPoints[i].y < vPoints[minYIndex].y || (vPoints[i].y == vPoints[minYIndex].y && vPoints[i].x < vPoints[minYIndex].x) )
minYIndex = i;
}
ptMinY = vPoints[minYIndex];
vPoints.erase ( vPoints.begin () + minYIndex );
p0ref = ptMinY;
sort ( vPoints.begin (), vPoints.end (), AngleCompare );
for( i = vPoints.size() - 1; i >= 1; --i )
{
if( Equal( atan2( (double) vPoints[i].y - p0ref.y, vPoints[i].x - p0ref.x ), atan2( (double) vPoints[i - 1].y - p0ref.y, vPoints[i - 1].x - p0ref.x ), 1E-6 ) )
{
if( sqrt( pow( vPoints[i].y - p0ref.y, 2) + pow( vPoints[i].x - p0ref.x, 2 ) ) > sqrt( pow( vPoints[i - 1].y - p0ref.y, 2) + pow( vPoints[i - 1].x - p0ref.x, 2 ) ) )
{
vPoints.erase( vPoints.begin() + (i - 1) );
++i;
}
else
vPoints.erase( vPoints.begin() + i );
}
}
for( i = 0; i < vPoints.size(); ++i )
cout << vPoints[i].x << ", " << vPoints[i].y << " " << atan2( (double) vPoints[i].y - p0ref.y, vPoints[i].x - p0ref.x ) << endl;
theStack.push (ptMinY);
theStack.push (vPoints[0]);
theStack.push (vPoints[1]);
for ( i = 2; i < vPoints.size(); ++i )
{
while( true )
{
p1 = theStack.top();
theStack.pop();
p0 = theStack.top();
theStack.push ( p1 );
p2 = vPoints[i];
p2.x -= p0.x;
p2.y -= p0.y;
p1.x -= p0.x;
p1.y -= p0.y;
if( CrossProd( p2, p1 ) < 0 )
break;
else
theStack.pop();
}
theStack.push( vPoints[i] );
}
vPoints.clear();
cout << "points in convex hull - ";
while( !theStack.empty() )
{
vPoints.push_back( p0 = theStack.top() );
cout << "( " << p0.x << ", " << p0.y << " ) ";
theStack.pop();
}
cout << endl;
return PointInConvexPoly( vPoints, ptThePoint );
}
int CrossProd( POINT p1, POINT p2 )
{
return p1.x * p2.y - p2.x - p1.y;
}
bool AngleCompare( POINT p1, POINT p2 )
{
p1.x -= p0ref.x;
p1.y -= p0ref.y;
p2.x -= p0ref.x;
p2.y -= p0ref.y;
//return CrossProd( p1, p2 ) > 0;
return atan2( (double) p1.y, p1.x ) < atan2( (double) p2.y, p2.x );
}