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.
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.
Before explaining the implementation, it should explain how are the data structures.
vector3fis a 3
floatvector with some geometrical functions (dotProduct, crossProduct, etc…)
facetis a 3
intvector 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
- 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.
- A points collection is a
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.
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.
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