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