I wrote a moderately popular Delaunay/Voronoi library for Unity a few years back [1] with a neat little demo video [2]. I implemented the incremental Bowyer-Watson algorithm for generating the triangulations, and then took the dual to generate the Voronoi tesselation (I also added a "clipper" that clips the voronoi diagram to a convex outline, which was fun, I haven't seen that anywhere else before and had to figure out how to do it myself).
Bowyer-Watson (which is described in this article) seems very simple to implement when you start: just start with a "big triangle" and then add points iteratively to it, and perform the flips you need to do. To do it performant, you have to build up a tree structure as you go, but it's not very tricky.
However: almost every description (including on Wikipedia) and implementation of Bowyer-Watson I could find was wrong. There's an incredibly subtle and hard to deal with issue with the algorithm, which is the initial "big triangle". Most people who implement it (and indeed I did the same initally) just make the triangle big enough to contain all the points, but that's not enough: it needs to be big enough to contain all the points in all the circumcircles of the triangles. These circumcircles can get ENORMOUS: in the limit of three points on a line, it's an infinite half-plane. Even if they aren't literally collinear, just close, the triangle becomes way to huge to deal with.
The answer is that you have to put the points "at infinity", which is a very weird concept. Basically, you have to have special rules for these points when doing comparisons, it's really tricky and very hard to get right.
If I were doing this again, I wouldn't use Bowyer-Watson, this subtlety is too tricky and hard to get right. Fortune's sweep-line is more complex on the surface, but that's the one I would go with. Or the 3D convex hull technique (which I also wrote a library for, by the way [3]).
I tried doing this by just assigning the points to [-inf, -inf, +inf, +inf, +inf, -inf], and I quickly ran into many, many NaNs, and was quite sad that it didn't work. It would be cool if floating point infinite numbers actually worked like you expect infinity to - although, perhaps that's just my shoddy understanding of infinity talking.
Bowyer-Watson (which is described in this article) seems very simple to implement when you start: just start with a "big triangle" and then add points iteratively to it, and perform the flips you need to do. To do it performant, you have to build up a tree structure as you go, but it's not very tricky.
However: almost every description (including on Wikipedia) and implementation of Bowyer-Watson I could find was wrong. There's an incredibly subtle and hard to deal with issue with the algorithm, which is the initial "big triangle". Most people who implement it (and indeed I did the same initally) just make the triangle big enough to contain all the points, but that's not enough: it needs to be big enough to contain all the points in all the circumcircles of the triangles. These circumcircles can get ENORMOUS: in the limit of three points on a line, it's an infinite half-plane. Even if they aren't literally collinear, just close, the triangle becomes way to huge to deal with.
The answer is that you have to put the points "at infinity", which is a very weird concept. Basically, you have to have special rules for these points when doing comparisons, it's really tricky and very hard to get right.
If I were doing this again, I wouldn't use Bowyer-Watson, this subtlety is too tricky and hard to get right. Fortune's sweep-line is more complex on the surface, but that's the one I would go with. Or the 3D convex hull technique (which I also wrote a library for, by the way [3]).
[1]: https://github.com/OskarSigvardsson/unity-delaunay
[2]: https://www.youtube.com/watch?v=f3T5jtsokz8
[3]: https://github.com/OskarSigvardsson/unity-quickhull