Solved

convex hull and point in convex hull algorithms

Posted on 2004-09-07
12
621 Views
Last Modified: 2008-01-09
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 );
}
0
Comment
Question by:Kitty__Kong
[X]
Welcome to Experts Exchange

Add your voice to the tech community where 5M+ people just like you are talking about what matters.

  • Help others & share knowledge
  • Earn cash & points
  • Learn & ask questions
  • 5
  • 2
12 Comments
 
LVL 1

Author Comment

by:Kitty__Kong
ID: 12002191
the other possible problem is with the part where the algorithm says while not making a left turn pop from the stack. i didnt know how to determine a left turn so i used google and found something that said a left turn occurrs if the cross product is < 0 that is the

               if( CrossProd( p2, p1 ) < 0 )
                    break;
               else
                    theStack.pop();
0
 
LVL 1

Author Comment

by:Kitty__Kong
ID: 12002200
i forgot

struct POINT
{
      int      x, y;
};

and p0ref is a global

POINT      p0ref;
0
 
LVL 2

Expert Comment

by:sneeuw_chan
ID: 12005182
First of all, IIRC, sort() requires a function that returns -1, 0 or 1, where 0 means 'equal', and + or -1 means 'larger' and 'smaller' respectively.  Your function seems to be returning a boolean.
- You've got a small error in your CrossProd there, should be p1.x * p2.y - p2.x * p1.y.
I haven't checked, but if you do it with arctangens, you may have problems around the X or Y axes.

Also, the problem with sorting is that you're sorting something that's in a circle, so you need to define an arbitrary direction where you 'open up' the circle, so to speak.  For example, your angle comparison could first check the X coordinates' signs against each other, and return +1 or -1 depending on that, and if they have the same sign (the angles are in the same hemisphere), then you check the angles against each other.  (Or maybe first check the signs of the Y coordinates also, and check the actual angles only when they're in the same quadrant.)
0
Technology Partners: We Want Your Opinion!

We value your feedback.

Take our survey and automatically be enter to win anyone of the following:
Yeti Cooler, Amazon eGift Card, and Movie eGift Card!

 
LVL 1

Author Comment

by:Kitty__Kong
ID: 12009670
according to msdn sort takes a function that returns a bool.

http://msdn.microsoft.com/library/default.asp?url=/library/en-us/vcstdlib/html/vclrfalgorithmsort.asp

how would i change AngleCompare to properly sort with CrossProd instead of using atan?
0
 
LVL 2

Accepted Solution

by:
sneeuw_chan earned 500 total points
ID: 12010657
Ah, yes.  I see.  That's different from the C sort() function, then.

Lessee, what I would do (off the top of my head, untested):

bool AngleCompare( POINT p1, POINT p2 )
{
     p1.x -= p0ref.x;
     p1.y -= p0ref.y;
     p2.x -= p0ref.x;
     p2.y -= p0ref.y;
     if ((p1.x < 0) && (p2.x > 0)) {
          return true;
     }
     if ((p1.x >= 0) && (p2.x =< 0)) {
          return false;
     }
     return CrossProd( p1, p2 ) > 0;
}

This way you make sure that the angles you're comparing are always less than 180 degrees.
If you order the points in a circle, you have to mark a point on the circle as the start and endpoint, otherwise you don't get a stable ordering.
0
 
LVL 1

Author Comment

by:Kitty__Kong
ID: 12579692
i object. i dont think what sneeuw_chan works. there are probably other problems in my algorithm but i dont know what so i need someone to point them out
0
 
LVL 1

Author Comment

by:Kitty__Kong
ID: 12600964
what does force accept mean? can i get a point refund or something?
0

Featured Post

Free Tool: SSL Checker

Scans your site and returns information about your SSL implementation and certificate. Helpful for debugging and validating your SSL configuration.

One of a set of tools we are providing to everyone as a way of saying thank you for being a part of the community.

Question has a verified solution.

If you are experiencing a similar issue, please ask a related question

Suggested Solutions

Title # Comments Views Activity
convert char array to number in c 5 93
How to convert MFC::CString to UTF8 wchar_t* 10 460
Indy 10 not Receiving UDP broadcast 3 56
print bytes of an integer 6 48
In days of old, returning something by value from a function in C++ was necessarily avoided because it would, invariably, involve one or even two copies of the object being created and potentially costly calls to a copy-constructor and destructor. A…
Written by John Humphreys C++ Threading and the POSIX Library This article will cover the basic information that you need to know in order to make use of the POSIX threading library available for C and C++ on UNIX and most Linux systems.   [s…
The viewer will learn additional member functions of the vector class. Specifically, the capacity and swap member functions will be introduced.
The viewer will learn how to clear a vector as well as how to detect empty vectors in C++.

751 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question