Tuesday, February 9, 2010

Simple Stupid Convex Hull

Sometimes you just need it, you know just to make sure your 7 point polygon is actually convex. Calculating convex hull would be awesome solution, but you kinda don't have an implementation at hand and you really would not like to write one as you are solving that other problem. You fire up Google, and end up browsing stupid animating Java implementations and get frustrated.

After browsing through dozen of sources using atan and the friends and that unsuccessful adventure of trying to remove the global variables from the Computational Geometry in C example code, and few attempts to find qsort_r that is actually portable, I landed on the following code with some help from Wikipedia.

It is Gift wrap (or Jarvis March), O(nh), does not need qsort_r, has pretty compact implementation, and simple enough to debug and fix too when it fails.

[Edit] Got bitten by floating points precision again. Updated the code so that it works better when you have large numbers. Updated formatting.

 // Returns true if 'c' is left of line 'a'-'b'.  
inline bool left(const float* a, const float* b, const float* c)
{
const float u1 = b[0] - a[0];
const float v1 = b[2] - a[2];
const float u2 = c[0] - a[0];
const float v2 = c[2] - a[2];
return u1 * v2 - v1 * u2 < 0;
}
// Returns true if 'a' is more lower-left than 'b'.
inline bool cmppt(const float* a, const float* b)
{
if (a[0] < b[0]) return true;
if (a[0] > b[0]) return false;
if (a[2] < b[2]) return true;
if (a[2] > b[2]) return false;
return false;
}
// Calculates convex hull on xz-plane of points on 'pts',
// stores the indices of the resulting hull in 'out' and
// returns number of points on hull.
int convexhull(const float* pts, int npts, int* out)
{
// Find lower-leftmost point.
int hull = 0;
for (int i = 1; i < npts; ++i)
if (cmppt(&pts[i*3], &pts[hull*3]))
hull = i;
// Gift wrap hull.
int endpt = 0;
int i = 0;
do
{
out[i++] = hull;
endpt = 0;
for (int j = 1; j < npts; ++j)
if (hull == endpt || left(&pts[hull*3], &pts[endpt*3], &pts[j*3]))
endpt = j;
hull = endpt;
}
while (endpt != out[0]);
return i;
}

11 comments:

  1. Thanks for sharing the code !

    I used StanHull in the past as drop-in solution for convex hull wihtout too much hassle:
    http://www.bulletphysics.org/Bullet/phpBB3/viewtopic.php?p=756&f=&t=

    ReplyDelete
  2. Isn't StanHull in 3D? This one does 2D only.

    ReplyDelete
  3. Melkman's Convex Hull Algorithm did it for me
    http://www.ams.sunysb.edu/~jsbm/courses/345/melkman.pdf

    Here you have one implementation in Scala
    http://www.box2d.org/forum/viewtopic.php?f=3&t=3298

    There is one bug in the Scala implementation I noticed when porting it to Java.

    ReplyDelete
  4. IIRC Melkman's algo needs the input points to be in such order that they for non-intersecting polyline. Easiest way to ensure that is to sort the verts and we kinda need that qsort again :)

    I have used Andrew's algo earlier a lot. It also needs sorting and the implementation is more complex.

    ReplyDelete
  5. For the record, Realtime Collision Detection introduces Andrews's algorithm and also Quickhull, see http://www.qhull.org/

    ReplyDelete
  6. Ah I thought you had the points ordered when you said you wanted to check if your 7 point polygon was convex ;)

    ReplyDelete
  7. Guardian, that's right, but it does not contain the sources, so no cut'n'paste :) Qhull is really good, but I don't like the api and it is a bit too big to include when you could do with less. Yes, I'm picky! ;)

    Obi, my bad, missleading example :) Actually I was using this for a tool which allows you to create convex shapes just by click points on an object (see some of the previous area screenshots).

    ReplyDelete
  8. I believe there was a novel convex hull processing class in a Gems 3 article.

    ReplyDelete
  9. Graphics Gems or GPG? TO my surprise, I could find any convex hull entry in GG.

    ReplyDelete
  10. In the cmppt and left methods you reference the a[0] a[2] but never a[1]. Why is this?

    The code works great, thank you for sharing. I've been trying to get my head wrapped around the cmppt and left methods because I am trying to convert this to work with b2vec2s but it simply outputs the strangest things after I've tried to convert it. :/

    ReplyDelete
  11. There is a hint in that pts[hull*3] :)

    Yeah, it is 3d points on XZ plane. If you want to use 2D points, replace [2] with [1] and *3 with *2.

    ReplyDelete