QuickHull 3D

At the Computational Geometry course I have implemented the Quickhull algorithm for its application in 3D, following the paper “Barber, C. B., Dobkin, D. P., & Huhdanpaa, H. T. (1996). The Quickhull algorithm for convex hulls. ACM Trans. on Mathematical Software, 22(4), 469—483”. And following the example of this java applet developed by UNSW School of Computer Science and Engineering.

Algorithm

The algorithm proposed by the paper says:

create a simplex of d+1 points

for each facet F

for each unassigned point p

if p is above F

assign p to F’s outside set

for each facet F with a non-empty outside set

select the furthest point p of F’s outside set

initialize the visible set V to F

for all unvisited neighbors N of facets in V

if p is above N

add N to V

the boundary of V is the set of horizon ridges H

for each ridge R in H

create a new facet from R and p

link the new facet to its neighbors

for each new facet F’

for each unassigned point q in an outside set of a facet in V

if q is above F’

assign q to F”s outside set

delete the facets in V

I made some changes on my implementation but the essence is the same.

quickhull-initial.jpg

Data Structures

Before explaining the implementation, it should explain how are the data structures.

  • vector3f is a 3 float vector with some geometrical functions (dotProduct, crossProduct, etc…)
  • facet is a 3 int vector and a std::Vector<int>, to store the indices of the 3 points of the triangle and the indices of the outside points.
  • All the collection used are class std::Vector
  • The QuickHull class has:
    • A points collection is a std::Vector <vector3f>, where there are stored all the points. This should not be changed any time.
    • A facets collection is a std::Vector <facet>, where there are stored all the facets.

Implementation

The first thing to do is calculate the simplex. The simplex is a triangle that separates points in two areas. This will take 3 points, two with maximum X and Y coordinate, and another with the minimum Z coordinate of the whole set of points.

quickhull-simplex.jpg

After the simplex is created, split the remaining points that are visible for the two facets from the simplex (a triangle in 3D has 2 faces). To know if a point is visible for a facet, we have to calculate if this point is above the facet. And to know if a point is above a plane, we evaluate the point on the plane normal.

Now the main loop begins. In each iteration and for each facet, the facet’s furthest point is calculated. This point p is checked whether is visible from the facets neighbourhood, the neighbourhood facets are the ones that have at least one vertex in common, then all ridges and outside points from the neighbourhood facets are saved in different collections.

The ridge collection determine the boundary from the hole that is created when the neighbourhood facets are deleted. Using this boundary and the furthest point p will create new triangles using the Delaunay algorithm.

Finally once created the new triangles will split the points of collection points on the new triangles created. This loop is executed while there are facets with outside points.

quickhull-final.jpg

Solved Facts

Triangle orientation

One of the problems that can happen when you are adding new triangles on a closed polygon, is that according to the order of the vertices of the triangles, the normal of the triangle is pointing inside the polygon. This is a problem because to perform calculations to find the location of the points on the triangle requires that the normal is facing outside the polygon.

The solution implemented, was using a reference point. This point must always remain within the polygon, even in the first iteration. Therefore, is calculated the baricenter of the initial simplex.

Here is the source code: SRC