Convex Hull, one algorithm implementation

I talked about convex hulls some time ago in an alpha shape post.

The convex hull is probably one of the most basic computational geometry algorithms, and because of that it is present in almost, if not all, geometry/cad/gis libraries and software packages. In this post you will find an explanation of one of the existing algorithms to compute it, an implementation with C++, plus a set of scripts to generate various point clouds and the corresponding hulls.

Objective

• Describe one of the possible convex hull algorithms
• Implement the algorithm in C++
• Provide various scripts to generate random point clouds and compute its convex hulls

I won’t lie to you, this post might be boooring ;-).

Convex hull

For me, the best description is: Imagine a board with nails and an elastic string. Place the elastic string covering all the nails and you have a convex hull. But if you want a more intelligent definition, check the wolfram alpha page for convex hull.

The algorithm

I am aiming to a usable algorithm. At least with a nice cost. I read in wikipedia that there are a lot of algorithms, and that the bottom complexity is of the order O(nlogn). With this in mind I will use an incremental algorithm to build the upper convex hull and the lower convex hull.

The trick here is: when walking the boundary of a polygon on a clockwise direction, on each vertex there is a turn left, or right. On the convex hull polygon, this turn will always be a right turn.

The presented algorithm is an incremental algorithm that will contain the upper hull for all the points treated so far.

n = number of points.
CHULLU = list of ordered points forming the upper hull.
CHULLL = list of ordered points forming the lower hull.
CHULL = list of points forming the convex hull.

• order the points by x coordinate. (if equal x, order by y)
• CHULLU = {} (a)
• for i=3 to n
• append to CHULLU
• while |CHULLU| > 2 and last 3 points of CHULLU make a left turn (b)
• Delete the penultimate point from CHULLU

The invariant in this case is that CHULLU contains the upper hull of the points processed so far. At the beginning there are only 3 points so the upper hull is the line joining these 3 points. We can add because it is the rightmost point seen so far, and this point has to be on the hull. After adding it may happen that CHULLU contains a left turn, in this case, keep deleting points from the CHULL until one of the two options is met: either the CHULLU does not contain a left turn, or there are only 2 points left.

For the lower hull, the process is the same, but traversing the point list from n-3 to 0.

Finally, CHULLU and CHULLL are joined into a single list CHULL.

Complexity

For the sake of doing it, here is the “proof” of the algorithm complexity.

First step, sorting the points can be achieved in O(nlogn) with any of the famously known algorithms: Quicksort or Mergesort for example.

The external for loop (a) is executed a specific number of times n. The difference will be in how many times the internal loop (b) is executed.
For each execution of loop (a), loop (b) is executed at least once, and each time we execute it a point is deleted.
Each point can only be deleted once, so this computation is also bounded by n.
Note that this argument also holds for the lower hull.

The hull construction loops have a complexity of O(n). Because of the sorting step, the algorithm has a time complexity of O(nlogn).

Right Turn

We need to check if 3 points make a right turn. Left turns and collinear points are not accepted.

You may want to review your vector skills, WikiHow has a nice refresh on finding angles between vectors with some nice images 🙂 even if you don’t need it, I think they did a cute job in that page.

The problem with the classic approach is that the result domain of is in not giving knowledge of the direction. And on the other hand, for this case I don’t care for the angle at all, just for the turn: left or right (counter clockwise or clockwise).

Having the two vectors with components , I can just calculate the determinant formed by those two vectors:

The sign will indicate the orientation. Negative indicates a clockwise direction, positive a counter clockwise direction. Check the Wikipedia determinants article.

Source

With the basic idea of the convex hull algorithm and the knowledge on how to check if 3 points (2 vectors) perform a right or left turn, the implementation should not be difficult.

I will not include all the sources in this post, you can visit my github for it. I will just explain the various source files and show the convex hull algorithm.

My solution has various files:

• Point.cpp/h
  Class defining point. Just x, y doubles , << operator and read from cin.
• Segment.cpp/h
  Class defining a segment from 2 points.
• geom.cpp/h
  Source file with various geometry functions needed to compute the convex hull: chull, iscw (is clockwise). It also contains other non used functions that I added while playing.
• main_chullgnuplot.cpp
  Read points from stdin (in gnuplot style) and prints the chull points.

The C++ function for computing the convex hull:

vector<Point> chull(const vector<Point>& pst){
    vector<Point> points(pst);
    vector<Point> lupper,llower,edges;
    typedef vector<Point>::size_type vcs;

    sort(points.begin(),points.end(),point_order_x);
    vector<Point>::const_iterator it;

    //Constructing the upper hull. 
    //The point with lowest x will be part of the hull
    lupper.push_back(points[0]);
    lupper.push_back(points[1]);

    //Loop the rest of the points to obtain the upper hull
    for(it=points.begin()+2;it!=points.end();++it){
        lupper.push_back(*it);

        //while size>2 and not a right turn
        while ((lupper.size() > 2) && !iscw(lupper))
            lupper.erase(lupper.end() - 2);
    }

    //Constructing the lower hull.
    it = points.end()-1;
    llower.push_back(*it);
    llower.push_back(*(it-1));

    for(it=points.end()-3;it>=points.begin();--it){
        llower.push_back(*it);

        //while size>2 and not a right turn
        while ((llower.size() > 2) && !iscw(llower))
            llower.erase(llower.end() - 2);

    }

    //First llower is already in lupper
    copy(lupper.begin(),lupper.end(),back_inserter(edges));
    copy(llower.begin()+1,llower.end(),back_inserter(edges)); 
    return edges;
}

I happily use C++ STD facilities like: sort, copy and back_inserter.

For the implementation I used vectors, but note how for the lupper and llower vectors I remove items from the middle of the vector. That does not play very nice with the STD vector implementation because every time I erase something it has to relocate all the elements. So if I want it to scale to work more efficiently I should reimplement it with another container, for example, lists.

Example execution

Once compiled the main file expects points from stdin and prints the chull to stdout.

10 10
20 20
30 30
5 6 

main_chullgnuplot.out < input.dat > input.dat.chull

5 6
30 30
10 10
5 6

Extra scripts

To test the algorithm I wrote some extra scripts:

• pointgen.py
  Generate a set of random points. (int and float versions)
• plot.p
  gnuplot script to plot the points and the convex hull. It expects pcloud, chull
  and outfile variables to be set.
• gnuplot.sh
  For each of the data files under the data directory, runs gnuplot with the plot.p script, and joins all the results in a single ps file.
• Makefile:
  make plots will drive all this.

Results

Finally. We discussed an algorithm, presented a C++ implementation for it, tested for a simple case and prepared a set of scripts to run a lot of times with different cases. Probably not degenerate cases, but I am happy with that right now :-). Next there is the result postscript after running make plots

Conclusions and Remarks

Not a very fancy conclusion here. I’ve been reading a little about a convex hull algorithm and applied it :-).

References

Latex in wordpress!⇒GO
Computational Geometry ⇒GO
Computational Geometry in C⇒GO
Java applet with convex hull algorithms example ⇒GO
Source code ⇒GO