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 );

}