Monogon's blog

By Monogon, history, 4 years ago, In English

Warning: The following contains graphic depictions of geometry, precision errors, and degenerate cases. Viewer discretion is advised.

Prerequisites

I assume the reader is familiar with:

  • 2D Convex Hulls
  • 3D Vector Operations (dot and cross products)

Introduction

Recall that in the 2D convex hull problem, you are given a set of 2D points, and you must compute the smallest convex polygon containing all the given points. By convex, we mean that for any two points $$$A$$$ and $$$B$$$ inside the polygon, the entire line segment $$$AB$$$ is also inside the polygon.

The problem in 3D is completely analogous. You are given $$$n$$$ 3D points, and you must compute the smallest convex polyhedron containing all the given points. Similarly, a polyhedron is convex if for any two points $$$A$$$ and $$$B$$$ inside the polyhedron, the line segment $$$AB$$$ is also inside the polyhedron.

In the 2D case, it is more obvious what the output is. We can simply output a circular list of vertices on the boundary of the polygon. But how do we represent a polyhedron? Well, we can represent it as a list of triangular faces. A non-triangular face can always be divided into triangles, so this requirement isn't limiting.

In the 2D case, there are some degenerate cases that are easy to miss if you're not careful. For example, if there are 3 collinear points, the polygon might have two adjacent sides of the same slope, or you might combine them into one segment. Or if two points have the same $$$x$$$ coordinate, you need to handle this carefully when sorting the points by $$$x$$$ coordinate. Or if all the points are on the same line, then the convex hull isn't even a non-degenerate polygon!

Now, degenerate cases are bad enough with two dimensions. When you enter the third dimension, it only gets worse. What if four points lie on the same plane? Since we are requiring triangular faces, we could triangulate this in multiple ways. Or maybe we could choose to have one triangle of three points and the forth point lies inside this triangle, so we ignore it. What if all the points are on the same line? Or on the same plane even? Then the convex hull isn't a non-degenerate polyhedron. I will ignore such degenerate cases for now, and revisit them when applicable.

Brute Force Algorithm in $$$O(n^4)$$$

Suppose that $$$n$$$ is very small, and we are guaranteed that no four points are coplanar. Then how can we compute the 3D convex hull in the dumbest way possible?

We could simply consider every triplet of three points $$$(\vec{a}, \vec{b}, \vec{c})$$$, and check if they create a triangular face on the convex hull. In order for it to be a face, the remaining points must "see" the same side of the triangle. In other words, if we consider the plane containing this triangle, the remaining points should lie on the same side of the plane. To compute which side of a plane a point $$$\vec{p}$$$ is on, we can first take a vector $$$\vec{q}=(\vec{b}-\vec{a})\times (\vec{c}-\vec{a})$$$ orthogonal to the plane.

Then the sign of $$$(\vec{p}-\vec{a})\cdot \vec{q}$$$ tells us the side of the plane. In particular, a result of $$$0$$$ tells us that $$$\vec{p}$$$ lies on the plane.

For each triplet, we perform this check with all points, so the total time complexity is $$$O(n^4)$$$. Because of its simplicity, this should be the approach when the constraints allow it. And by examining the brute force case, we learned how to perform the most fundamental operation in any 3D convex hull algorithm: checking which side of a plane a point is on.

Incremental Algorithm in $$$O(n^2)$$$

Now we want a more practical algorithm. The strategy of the incremental algorithm is to build the convex hull for the first $$$i$$$ points, as $$$i$$$ goes from $$$1$$$ to $$$n$$$. All we have to do is figure out how to update the convex hull in order to accommodate one new point.

Let's start by making an analogy with the 2D convex hull. Suppose we currently have the convex hull of the first $$$i-1$$$ points, and we wish to add the $$$i$$$-th point. Well, if the point is already inside the hull, we should do nothing. Otherwise, let's pretend we are looking at the polygon from the perspective of the $$$i$$$-th point. We should delete all line segments it can see, and add two new line segments: one from the new point to each of the extreme points.

A similar thing happens in the 3D case. If a point is already inside the polyhedron, we simply do nothing. Otherwise, we delete all faces the new point can see. But what's less obvious is what we should add. Well, we've left the polyhedron with a connected set of faces removed. This exposes a cycle of edges. For each of these exposed edges, we should create a new face with the new point and that edge. Effectively, we create a cone of faces to repair the opening.

Now, let's analyze the time complexity. For each iteration, we process all faces of the hull in that iteration. And it can be proven that the number of faces is at most $$$O(n)$$$. For the proof, we use Euler's Formula. It states that for any polyhedron with $$$V$$$ vertices, $$$E$$$ edges, and $$$F$$$ faces, $$$V-E+F=2$$$. All faces are triangles, so we can substitute $$$E=3F/2$$$ since each face has $$$3$$$ edges, and we count each edge twice for the $$$2$$$ faces it touches. Then we have $$$V-F/2=2$$$ and $$$F=2V-4=O(n)$$$. Since we process $$$O(n)$$$ faces in $$$O(n)$$$ iterations, the overall time complexity is $$$O(n^2)$$$.

Implementation

Clarkson-Shor Algorithm in $$$O(n\log n)$$$

There exist deterministic algorithms to compute the 3D convex hull in $$$O(n\log n)$$$ time. Most famously, there is a divide-and-conquer algorithm that solves it, but it is notoriously difficult to implement correctly. But don't lose hope! There is also a randomized algorithm based on the same incremental strategy, which takes $$$O(n\log n)$$$ time in expectation if the points are added in random order.

Unfortunately, we can't just shuffle the points, run the usual incremental algorithm, and call it a day. In order to achieve $$$O(n\log n)$$$, we should eliminate extra work of processing faces invisible to the new point. How can we do this? Well, before we checked for each point, which faces it can see. Instead, we will maintain for each point a list of faces it can see. Then when a face is created, we add it to the lists of all points that will be able to see it.

"Wait a minute," I hear you ask, "isn't this still $$$O(n^2)$$$ because you still process the same thing for each face-point pair, just in a different order?" The answer is yes, we're still in $$$O(n^2)$$$ territory. But, using some cleverness we can avoid checking with every future point every time. The key idea is that a new face $$$F$$$ is attached to a preexisting edge $$$e$$$, and a point $$$p$$$ can see $$$F$$$ only if it was able to see one of the two faces of $$$e$$$. Instead of checking with all future points, we will only check the points in the lists associated with these two faces. Unfortunately, this requires more work to implement. This time, we need to start caring about how the faces touch each other in order to find these two lists.

Precision Errors and Degenerate Cases

The only potential source for precision error is our check of which side of a plane a point is on. If all coordinates are integers, then our computations will stay as integers. If it overflows on 64-bit integer types, we can convert to floating points for computations, and set $$$\epsilon=0.5$$$. Note that if the coordinates are at most $$$X,Y,Z$$$ in absolute value (in the $$$x$$$, $$$y$$$, and $$$z$$$ directions, respectively), our computations are on the order of $$$XYZ$$$. Use this to determine whether it will overflow on integers or find what you should set $$$\epsilon$$$ to. Or just try every possible $$$\epsilon$$$ until it works.

To handle degenerate cases, we should make sure that the first $$$4$$$ points are not coplanar. If all the $$$n$$$ points do not lie in the same plane, we can find $$$4$$$ such points and move them to the beginning. Once this is taken care of, the only real degenerate case remaining is when we need to decide if a face is visible to a point in the same plane. I made the decision to consider the face invisible in such a case. This ensures that the algorithm won't unnecessarily subdivide faces. However, it is still possible that the algorithm can have non-unique output depending on the order the points are processed. This is inevitable since there is not a unique way to triangulate faces in degenerate cases. If we want to, though, we can combine adjacent coplanar faces afterward.

To my knowledge, the Clarkson-Shor algorithm doesn't make any guarantees about $$$O(n\log n)$$$ performance when there are degenerate cases. I think by considering coplanar faces invisible, we are safe. But if anyone has a nasty test case against it, I'm interested to hear it.

Here is my implementation of the $$$O(n\log n)$$$ algorithm with full handling of degenerate cases (I hope). If anyone finds bugs or ways to improve the efficiency and/or readability of this implementation, please let me know!

Implementation

Application: Delaunay Triangulation and Voronoi Diagrams

Suppose we have $$$n$$$ 2D points $$$p_1,\ldots,p_n$$$. We may want to answer queries of the form "Which point is closest to $$$q_i$$$?" for $$$m$$$ points $$$q_1,\ldots,q_m$$$.

The Voronoi diagram divides the plane into $$$n$$$ regions in a way that can answer these kinds of queries. It is defined so that the region corresponding to the point $$$i$$$ consists of all points closer to $$$p_i$$$ than any other of the $$$n$$$ points. That is, it's the region of points whose nearest neighbor is $$$p_i$$$. The Voronoi diagram can be represented by $$$O(n)$$$ line segments that border the regions (some line segments will extend to a point at infinity). We can answer our queries by sorting all line segments and queries, and performing a sweep line.

The dual of the Voronoi diagram is called the Delaunay Triangulation. It is the graph that connects two points if their Voronoi cells share a boundary edge. Conversely, the points of a Voronoi diagram are the circumcenters of Delaunay triangles, and an edge in the Voronoi diagram connects the circumcircles of two adjacent Delaunay triangles.

The Delaunay triangulation has several nice properties that come in handy. Most notably, the Euclidean minimum spanning tree consists of a subset of Delaunay edges. A nice consequence of implementing 3D convex hull is that we get Delaunay triangulation for free. We can simply map each point $$$(x,y)$$$ into a 3D point $$$(x,y,x^2+y^2)$$$. Then the downward-facing triangles of the 3D convex hull are precisely the Delaunay triangles. The proof is left as an exercise to the reader.

Resources

Practice Problems

  • Vote: I like it
  • +393
  • Vote: I do not like it

| Write comment?
»
4 years ago, # |
  Vote: I like it +20 Vote: I do not like it

dang it, i didn't see the warning, now i will have nightmare about 3d convex hulls D:

»
4 years ago, # |
  Vote: I like it +96 Vote: I do not like it

4D convexhull tutorial when

»
4 years ago, # |
  Vote: I like it +35 Vote: I do not like it

Whoa, this is cool. What about speed? Does it work quickly for around 500,000 points?

Btw. it's worth noting that there exists an easy $$$O(n^3)$$$ brute force. Instead of testing triples of vertices as sides, test pairs of vertices as possible edges of the convex hull.

If all coordinates are integers, then our computations will stay as integers. If it overflows on 64-bit integer types, we can convert to floating points

I wouldn't expect coordinates bigger than $$$10^9$$$ in the input so __int128 should be enough to handle the cross product. If this type isn't available, we have a choice between floating values and implementing an easy version of big integers — represent a number as $$$x \cdot 10^9 + y$$$ where $$$y < 10^9$$$ and $$$x$$$ fits in long long. This is enough for numbers up to around $$$10^{28}$$$ and addition/subtraction/multiplication are easy.

  • »
    »
    4 years ago, # ^ |
    Rev. 2   Vote: I like it +10 Vote: I do not like it

    What about speed? Does it work quickly for around 500,000 points?

    I tested on 500,000 random points, and it seems to take between 5 and 8 seconds (on my own computer) (varying by choice of integers or floating points). As for the speed for non-random points, I have no idea.

    About integer overflow, specifically, when computing Delaunay triangulation, the z-coordinate can be much larger than the input coordinates. If you want to compute Delaunay triangulation with coordinates larger than $$$10^5$$$, you'll probably need floating points or bigger integers.

  • »
    »
    4 years ago, # ^ |
    Rev. 2   Vote: I like it +8 Vote: I do not like it

    Can you elaborate a bit more about easy $$$\mathcal{O}(n^3)$$$ bruteforce?

    How do you actually test a pair of points $$$(a, b)$$$ to be an edge of convex hull? For it to be true there must exist a plane, that contains both $$$a$$$ and $$$b$$$ and splits the space into halfspaces in such way that all other points belong to one of them. But with $$$3$$$ points this plane is well defined, with $$$2$$$ points it is not.

    I can only think of this algorithm:

    For each other point $$$p$$$ find point $$$p'$$$ that belongs to line $$$ab$$$ and is closest to $$$p$$$. Line $$$pp'$$$ is perpendicular to $$$ab$$$. Consider all $$$pp'$$$ as free vectors, they all belong to one plane. Now we have to check that all of these vectors can be covered with a single halfplane, or speaking in terms of polar angles all these vectors belong to sector of angle less than $$$\pi$$$. To achieve that we can add vectors one by one keeping two of them (lets say $$$i$$$, $$$j$$$) with greatest angle between. When we add next vector $$$k$$$ we check that $$$angle(i,j) + angle(j,k) + angle(k,i) < 2\pi$$$ and then keep as new $$$i, j$$$ pair with $$$max(angle(i,j), angle(j,k), angle(k,i))$$$.

    This check is $$$\mathcal{O}(n)$$$ but it still requires some mess with floating point angles and doesn't seem to be reliable and simple. What am I missing here?

    Also, is seems for me that $$$3D$$$ convex hull is not really useful in form of set of all edges. Form of set of all faces allows checking weather point lies inside convex hull, decomposing hull into tetrahedrons to compute volume or perform other manipulations. How do you use hull in form of edges? Or do you always convert it into set of faces form?

    To convert set of edges to set of faces we do something like: iterate through all points as $$$p$$$ that belong to convex hull, find all points that share an edge with $$$p$$$, these points must form a cycle: every point must share an edge with exactly $$$2$$$ other points from this set of neighbours. Recover this cycle and add all faces $$$(p, u, v)$$$ where $$$u$$$ and $$$v$$$ are points from the cycle that share an edge. This seems to simple enough and able to recover all faces that touch $$$p$$$ in $$$\mathcal{O}(n)$$$.

    • »
      »
      »
      4 years ago, # ^ |
        Vote: I like it +8 Vote: I do not like it

      I agree that finding the set of edges doesn't seem to be very useful on its own. And the $$$O(n^2)$$$ incremental algorithm is simple enough to implement anyway.

    • »
      »
      »
      4 years ago, # ^ |
        Vote: I like it +8 Vote: I do not like it

      It's easier than that and doesn't require floating values. It's just like testing whether a point $$$A$$$ belongs to a convex hull in 2D. We only use sign of cross product in both versions (2D and 3D).

      Slow version: sort the other $$$n-1$$$ points by angle with respect to $$$A$$$. Check if the first point in this order is indeed on the left from each of remaining $$$n-2$$$ points.

      Fast version: instead of sorting, find the leftmost-by-angle point in $$$O(n)$$$. Then again check if it's on the left from remaining $$$n-2$$$ points.

      Similarly, let's check in $$$O(n)$$$ whether an edge $$$AB$$$ belongs to CH of given $$$n$$$ points in 3D. Sort the other $$$n-2$$$ points by angle with respect to $$$AB$$$ or just linearly find the leftmost point, say $$$C$$$. Check if indeed $$$C$$$ is on the left from everything. In other words, check if $$$ABC$$$ is a side of CH. So by testing edges we got sides anyway.

      • »
        »
        »
        »
        4 years ago, # ^ |
          Vote: I like it 0 Vote: I do not like it

        I think I figured out all the moments and I hope I am understanding things correctly.

        The part that confused me most and pushed to think about things like atan2 was this problem: how do we find this leftmost-by-angle vector relying only on sign of cross product in case when all parts of circle are populated. How do we even define leftmost-by-angle vector when for each vector there we can find another vector that is located to the left from it, due to their circular nature. And correct answer is that in this case we can leave leftmost vector undefined and immediately say that current $$$2D$$$ point or $$$3D$$$ edge does not belong to convex hull.

        In case when there exists leftmost-by-angle vector, such vector that all other vectors are located to right from this, we can for sure can find it. Slow but straightforward $$$\mathcal{O}(n^2)$$$ algorithm will be to calculate located-left relation between all pairs. That will clearly find leftmost-by-angle vector or say that it does not exist. Boosting this to $$$\mathcal{O}(n)$$$ or $$$\mathcal{O}(n \cdot \log n)$$$, we can find "left maximum" or sort vectors with "$$$a$$$ is located left to $$$b$$$" comparator. In case when leftmost-by-angle vector does not exist both of these approaches will leave us with completely undefined behavior, to fix this at the end we can perform complete linear check wheather this vector is really leftmost-by-angle.

        It is also it is clear that leftmost-by-angle vector will give us $$$2D$$$ edge / $$$3D$$$ face of convex hull due to the fact that borders of some shadow can only be projected from borders of original solid body. But, out of curiousity I will also ask, are there some tricks or theory that will guarantee us visiting all faces of convex hull just by picking leftmost-by-angle point $$$C$$$ for each edge $$$AB$$$? Or is it still better to do fair recovery from set of all edges?

        Thank you for clarification.

        • »
          »
          »
          »
          »
          4 years ago, # ^ |
            Vote: I like it +8 Vote: I do not like it

          Yup, but keep in mind that sort() assumes that comparator satisfies some properties (like $$$a < b$$$ and $$$b < c$$$ implies $$$a < c$$$) and you might get RTE if you use sort() with respect to some point inside, I think. For both time and correctness, better use the version of linearly finding the leftmost-by-angle and checking it.

          are there some tricks or theory that will guarantee us visiting all faces of convex hull just by picking leftmost-by-angle point C for each edge AB

          Every triangular face will actually be found three times. Just make sure that you iterate over all $$$n \cdot (n-1)$$$ ordered edges.

»
4 years ago, # |
  Vote: I like it +127 Vote: I do not like it

Live a healthy life , stay away from Geometry

»
4 years ago, # |
  Vote: I like it +11 Vote: I do not like it

An absolutely orgasmic tutorial! Thanks!

»
4 years ago, # |
  Vote: I like it 0 Vote: I do not like it

Can anybody share me the correct code of the problem Stars in a Can?I got a wrong answer.I have spent many time in it but I still don't know where my code is wrong.And I can't find other people's correct code in the Internet.

  • »
    »
    4 years ago, # ^ |
      Vote: I like it 0 Vote: I do not like it

    Oh,thank you! I get a correct code and find my error.