6 ms·

Delaunay and Voronoi on a sphere

250 points by guiambros 8y ago

techeng 8y ago
Oh nice, this used to be my go-to interview question for a numerics team at <generic-tech-giant>. Solid write up!
gus_massa 8y ago
> But how do you computer a circumcenter for a triangle on a sphere? I found a solution on stackoverflow and implemented it. Why is it necessary to calculate the circumcenters in 3D? Why not calculating them in the projection to the plane and then back-projecting them to the sphere? Perhaps it's necessary to calculate the Voronoi regions in the plane too.
- gugagore 8y ago
  It might be that the triangle is big enough (relative to, say, the radius of the sphere) that that approximation doesn't apply.
- slavik81 8y ago
  It's pretty easy to do it in 3D. Doing anything on a 2D projection is hard, because the space is so distorted. When you can easily measure distances, this is a lot easier. It's a little weird that their results are inside the sphere. I don't fully understand the equation, but it looks to me like it's using angles/rotations from a to the centre, so I would expect it to remain on the surface. Maybe the problem is just numerical error? Of course, if they're projecting the solution back to the surface anyways, the simplest way would be to find the centre of the triangle in xyz space, ignoring the sphere, then projecting the result back to the surface. As long as the triangle isn't larger than half the circumference, the shortest distance through the sphere would correspond to the shortest distance on the surface. To be honest, I'm more surprised they projected to a plane in the first place. If he was doing it all from scratch, it would probably be simpler to do the triangulation in 3D.
  - frenchie14 8y ago
    I'm assuming the loss of height comes from measuring distances directly between points instead of along the surface of the sphere. Much cheaper, but by design you'll cut through the sphere to connect each point
    - slavik81 8y ago
      I think you're right. The circumcentre equation is just more complicated than the centroid calculation I was imagining. They're doing exactly what I had suggested in the third paragraph. There may also be a nice derivation of the triangle centre on the sphere coming from quaternion math that would avoid renormalization. I wonder...
  - gus_massa 8y ago
    I was thinking about having a true Voronoi diagram in the distorted version of the plane and a distorted Vornoi diagram in the sphere. I'm not sure if the result would be good enough to have all the expected properties. I'm not even sure that it's ok to do the triangulation in the plane as described in the article, perhaps mapping the triangles in the plane to distorted triangles in the sphere may case them to be too curved.
- amitp 8y ago
  In the projection, some triangles are rather spread out, and some are even going to cross from -Infinity to +Infinity. I didn't want to handle that wrap-around-infinity case and thought it'd be easier to calculate them on the sphere instead of on the projection.
tame73_binds 8y ago
I've used CSSGRID for this http://www.ncarg.ucar.edu/ngmath/cssgrid/contents.html http://www.ncarg.ucar.edu/ngmath/cssgrid/contents.html It's solid and easy to bind from anything with an FFI.
- Crespyl 8y ago
  It took me a bit to realize that your link has nothing to do with the other CSS Grid: https://developer.mozilla.org/en-US/docs/Web/CSS/CSS_Grid_Layout https://developer.mozilla.org/en-US/docs/Web/CSS/CSS_Grid_La...
icoder 8y ago
Nice! I was looking for this the other day, would have been nice to find a writeup on this site about it (since its writeups are always so clear). I ended up doing hexagons on a sphere (with 12 pentagons). There's a writeup on this site about that too. It goes to great length going around the pentagons, but for my purpose they're not an issue.
vanderZwan 8y ago
> For most of my Voronoi map generation projects, I don’t actually use Voronoi, but something similar. I move the points from the circumcenters to the centroids. Not only is this cheaper to calculate, it ensures that all the region edges have a reasonable length. With Voronoi, some of the edges have small or zero length, and that makes it harder for me to place roads and rivers using those edges. I have written more about this topic here. He is basically doing a form of relaxatation similar or equivalent to Lloyd's Algorithm[0], no? Adrian Secord used this way back in 2002 already to create stippling images[1], see also Mike Bostock's notebook[2]. Would iterating a few times maybe improve results even more? Tangent: I have been experimenting on and off with implementing more advanced variants of Secord's stippling algorithm in the past two weeks[3]. It's been a lot of fun! A bit too much at times, Observables makes the interactive feedback loop a bit too easy. [0] https://beta.observablehq.com/@mbostock/lloyds-algorithm https://beta.observablehq.com/@mbostock/lloyds-algorithm [1] https://cs.nyu.edu/~ajsecord/msc_thesis/ajsecord_msc_thesis.pdf https://cs.nyu.edu/~ajsecord/msc_thesis/ajsecord_msc_thesis.... [2] https://beta.observablehq.com/@mbostock/voronoi-stippling https://beta.observablehq.com/@mbostock/voronoi-stippling [3] https://beta.observablehq.com/@jobleonard/untitled https://beta.observablehq.com/@jobleonard/untitled (Yes, that still needs a lot of optimization. I know)
- thrower123 8y ago
  Yep, he is using a form of Lloyd relaxation, as detailed in the original polygonal map generator article http://www-cs-students.stanford.edu/~amitp/game-programming/polygon-map-generation/ http://www-cs-students.stanford.edu/~amitp/game-programming/...
- mxfh 8y ago
  Philippe Rivière did quite a number of iterations[1] on d3/Bostock's geo-voronoi, culminating in a proposed Voronoi Projection method[2]. [1]https://beta.observablehq.com/collection/@fil/voronoi https://beta.observablehq.com/collection/@fil/voronoi [2]https://projet.liris.cnrs.fr/mi2/posts/2018/10/05/voronoi-projection-en.html https://projet.liris.cnrs.fr/mi2/posts/2018/10/05/voronoi-pr...
  - 8y ago
blt 8y ago
Red blob is one of the best programming sites! I love the interactive demos and the clear writing style. They always strike a good balance between rigor and intuition.
- amitp 8y ago
  Thanks!
  - ChrisClark 8y ago
    Your articles are so good I've felt the need to actually build a game around one of them. Because unfortunately I don't need any of the articles for my current game. :(
vbarrielle 8y ago
While the stereographic prohection trick sounds neat, it looks like it could produuce a bad triangulation if the input points are distributed non-uniformly. Indeed the 2d Delaunay triangulation algorithm tries to build the most uniform triangulation, but the projection step creates distortion which means a uniform projected triangulation might map to a non-uniform 3d triangulation. I'm not sure it is an issue because some other factor may compensate this issue, but there should be a rigorous analysis of the stability of this approach to various inputs.
- amelius 8y ago
  I'm not sure if delaunay by itself would produce "regular" triangulations. To get some kind of regularity, you'd generally need a refinement step where more points are inserted. Of course, this could be difficult with the given projection.
  - vbarrielle 8y ago
    Yeah I used regular as a shorthand but nostly I meant that it would flip edges to avoid thin triangles.
- amitp 8y ago
  See section 8.5 of http://www.cis.upenn.edu/~cis610/convex8.pdf http://www.cis.upenn.edu/~cis610/convex8.pdf for the analysis showing that the stereographic projection is the projection that preserves the properties needed for Delaunay triangulation. I learned this technique from Philippe Riviere, who uses it in the d3-geo-voronoi library - https://beta.observablehq.com/collection/@fil/voronoi https://beta.observablehq.com/collection/@fil/voronoi
hamoid 8y ago
By accident I figured out how to draw a simple Voronoi pattern: draw lots of planes at the same distance from a center point, each plane perpendicular to that center point. https://www.openprocessing.org/sketch/611667 https://www.openprocessing.org/sketch/611667 (click </> to see the code). Probably someone did it before already :)
- kaffeemitsahne 8y ago
  Sounds similar to drawing intersecting cones at each point, which is a pretty common method I think. For example: http://alexbeutel.com/webgl/voronoi.html http://alexbeutel.com/webgl/voronoi.html A nice thing about that is if you draw cones with just 4 sides (so, pyramids..) you get a manhattan distance voronoi pattern (or something that looks very similar at least).
PezzaDev 8y ago
Tessellating a sphere is one of those never ending areas of study because no solution feels 100% satisfying
rosshemsley 8y ago
Perhaps an easier (and more robust) way to do this is as follows: 1) Take your points in 3D, and add the point at the origin 2) Compute Delaunay in 3D 3) The 2D faces not containing the origin is your triangulation. To get the Voronoi tessellation: 1) Circulate the faces of each vertex in order (use orientation test if they are not ordered yet) 2) Connect the circumcenters of each 2D face you visit. Voilà!
- vanderZwan 8y ago
  Well, conceptually maybe, but I have some doubts regarding: > 2) Compute Delaunay in 3D How much slower is 3D triangulation compared to 2D? How does it scale with more points? And just as importantly: how complex does the algorithm get when we want to improve performance? I'm not saying you're wrong - if anything I hope I'm wrong, but I suspect the answers to my questions are "significant", "badly", and "very", which would mean this is not easier/more robust in practice (or it comes at the price of being unacceptably slow). The reason I think so is Julia's Voronoi/Delaunay package[0]. It implements an algorithm from an extremely dense (to non-astronomers) astronomy paper[1]. Even the name is intimidating: "E pur si muove: Galiliean-invariant cosmological hydrodynamical simulations on a moving mesh" Page five and six of the paper lists a short history of triangulation algorithms, making it clear that Delaunay triangulation so important for scientific research that many people have given it considerable thought. Sadly, I cannot really follow the algorithm explained in that paper. It seems to fall for the inevitable "algorithm by PhDs with complex enough math that you need a comparable PhD to being to grok it"-curse (and/or I am just too dumb). Of course, if someone already did the hard work of implementing really fast and robust 3D Delaunay libraries to steal from... ;) (the Julia package only implements the 2D algorithm. Sadly) [0] https://github.com/JuliaGeometry/VoronoiDelaunay.jl https://github.com/JuliaGeometry/VoronoiDelaunay.jl [1] https://arxiv.org/abs/0901.4107 https://arxiv.org/abs/0901.4107
  - contravariant 8y ago
    Computing Delaunay on a sphere isn't exactly the easiest either, distances get rather tricky and there's no simple projection that preserves all of them.

250 points by guiambros 8y ago

Oh nice, this used to be my go-to interview question for a numerics team at <generic-tech-giant>. Solid write up!

> But how do you computer a circumcenter for a triangle on a sphere? I found a solution on stackoverflow and implemented it. Why is it necessary to calculate the circumcenters in 3D? Why not calculating them in the projection to the plane and then back-projecting them to the sphere? Perhaps it's necessary to calculate the Voronoi regions in the plane too.

gugagore 8y ago
It might be that the triangle is big enough (relative to, say, the radius of the sphere) that that approximation doesn't apply.
slavik81 8y ago
It's pretty easy to do it in 3D. Doing anything on a 2D projection is hard, because the space is so distorted. When you can easily measure distances, this is a lot easier. It's a little weird that their results are inside the sphere. I don't fully understand the equation, but it looks to me like it's using angles/rotations from a to the centre, so I would expect it to remain on the surface. Maybe the problem is just numerical error? Of course, if they're projecting the solution back to the surface anyways, the simplest way would be to find the centre of the triangle in xyz space, ignoring the sphere, then projecting the result back to the surface. As long as the triangle isn't larger than half the circumference, the shortest distance through the sphere would correspond to the shortest distance on the surface. To be honest, I'm more surprised they projected to a plane in the first place. If he was doing it all from scratch, it would probably be simpler to do the triangulation in 3D.
- frenchie14 8y ago
  I'm assuming the loss of height comes from measuring distances directly between points instead of along the surface of the sphere. Much cheaper, but by design you'll cut through the sphere to connect each point
  - slavik81 8y ago
    I think you're right. The circumcentre equation is just more complicated than the centroid calculation I was imagining. They're doing exactly what I had suggested in the third paragraph. There may also be a nice derivation of the triangle centre on the sphere coming from quaternion math that would avoid renormalization. I wonder...
- gus_massa 8y ago
  I was thinking about having a true Voronoi diagram in the distorted version of the plane and a distorted Vornoi diagram in the sphere. I'm not sure if the result would be good enough to have all the expected properties. I'm not even sure that it's ok to do the triangulation in the plane as described in the article, perhaps mapping the triangles in the plane to distorted triangles in the sphere may case them to be too curved.
amitp 8y ago
In the projection, some triangles are rather spread out, and some are even going to cross from -Infinity to +Infinity. I didn't want to handle that wrap-around-infinity case and thought it'd be easier to calculate them on the sphere instead of on the projection.

tame73_binds 8y ago

I've used CSSGRID for this http://www.ncarg.ucar.edu/ngmath/cssgrid/contents.html http://www.ncarg.ucar.edu/ngmath/cssgrid/contents.html It's solid and easy to bind from anything with an FFI.

Crespyl 8y ago
It took me a bit to realize that your link has nothing to do with the other CSS Grid: https://developer.mozilla.org/en-US/docs/Web/CSS/CSS_Grid_Layout https://developer.mozilla.org/en-US/docs/Web/CSS/CSS_Grid_La...

icoder 8y ago

Nice! I was looking for this the other day, would have been nice to find a writeup on this site about it (since its writeups are always so clear). I ended up doing hexagons on a sphere (with 12 pentagons). There's a writeup on this site about that too. It goes to great length going around the pentagons, but for my purpose they're not an issue.

vanderZwan 8y ago

> For most of my Voronoi map generation projects, I don’t actually use Voronoi, but something similar. I move the points from the circumcenters to the centroids. Not only is this cheaper to calculate, it ensures that all the region edges have a reasonable length. With Voronoi, some of the edges have small or zero length, and that makes it harder for me to place roads and rivers using those edges. I have written more about this topic here. He is basically doing a form of relaxatation similar or equivalent to Lloyd's Algorithm[0], no? Adrian Secord used this way back in 2002 already to create stippling images[1], see also Mike Bostock's notebook[2]. Would iterating a few times maybe improve results even more? Tangent: I have been experimenting on and off with implementing more advanced variants of Secord's stippling algorithm in the past two weeks[3]. It's been a lot of fun! A bit too much at times, Observables makes the interactive feedback loop a bit too easy. [0] https://beta.observablehq.com/@mbostock/lloyds-algorithm https://beta.observablehq.com/@mbostock/lloyds-algorithm [1] https://cs.nyu.edu/~ajsecord/msc_thesis/ajsecord_msc_thesis.pdf https://cs.nyu.edu/~ajsecord/msc_thesis/ajsecord_msc_thesis.... [2] https://beta.observablehq.com/@mbostock/voronoi-stippling https://beta.observablehq.com/@mbostock/voronoi-stippling [3] https://beta.observablehq.com/@jobleonard/untitled https://beta.observablehq.com/@jobleonard/untitled (Yes, that still needs a lot of optimization. I know)

thrower123 8y ago
Yep, he is using a form of Lloyd relaxation, as detailed in the original polygonal map generator article http://www-cs-students.stanford.edu/~amitp/game-programming/polygon-map-generation/ http://www-cs-students.stanford.edu/~amitp/game-programming/...
mxfh 8y ago
Philippe Rivière did quite a number of iterations[1] on d3/Bostock's geo-voronoi, culminating in a proposed Voronoi Projection method[2]. [1]https://beta.observablehq.com/collection/@fil/voronoi https://beta.observablehq.com/collection/@fil/voronoi [2]https://projet.liris.cnrs.fr/mi2/posts/2018/10/05/voronoi-projection-en.html https://projet.liris.cnrs.fr/mi2/posts/2018/10/05/voronoi-pr...
- 8y ago

blt 8y ago

Red blob is one of the best programming sites! I love the interactive demos and the clear writing style. They always strike a good balance between rigor and intuition.

amitp 8y ago
Thanks!
- ChrisClark 8y ago
  Your articles are so good I've felt the need to actually build a game around one of them. Because unfortunately I don't need any of the articles for my current game. :(

vbarrielle 8y ago

While the stereographic prohection trick sounds neat, it looks like it could produuce a bad triangulation if the input points are distributed non-uniformly. Indeed the 2d Delaunay triangulation algorithm tries to build the most uniform triangulation, but the projection step creates distortion which means a uniform projected triangulation might map to a non-uniform 3d triangulation. I'm not sure it is an issue because some other factor may compensate this issue, but there should be a rigorous analysis of the stability of this approach to various inputs.

amelius 8y ago
I'm not sure if delaunay by itself would produce "regular" triangulations. To get some kind of regularity, you'd generally need a refinement step where more points are inserted. Of course, this could be difficult with the given projection.
- vbarrielle 8y ago
  Yeah I used regular as a shorthand but nostly I meant that it would flip edges to avoid thin triangles.
amitp 8y ago
See section 8.5 of http://www.cis.upenn.edu/~cis610/convex8.pdf http://www.cis.upenn.edu/~cis610/convex8.pdf for the analysis showing that the stereographic projection is the projection that preserves the properties needed for Delaunay triangulation. I learned this technique from Philippe Riviere, who uses it in the d3-geo-voronoi library - https://beta.observablehq.com/collection/@fil/voronoi https://beta.observablehq.com/collection/@fil/voronoi

hamoid 8y ago

By accident I figured out how to draw a simple Voronoi pattern: draw lots of planes at the same distance from a center point, each plane perpendicular to that center point. https://www.openprocessing.org/sketch/611667 https://www.openprocessing.org/sketch/611667 (click </> to see the code). Probably someone did it before already :)

kaffeemitsahne 8y ago
Sounds similar to drawing intersecting cones at each point, which is a pretty common method I think. For example: http://alexbeutel.com/webgl/voronoi.html http://alexbeutel.com/webgl/voronoi.html A nice thing about that is if you draw cones with just 4 sides (so, pyramids..) you get a manhattan distance voronoi pattern (or something that looks very similar at least).

PezzaDev 8y ago

Tessellating a sphere is one of those never ending areas of study because no solution feels 100% satisfying

rosshemsley 8y ago

Perhaps an easier (and more robust) way to do this is as follows: 1) Take your points in 3D, and add the point at the origin 2) Compute Delaunay in 3D 3) The 2D faces not containing the origin is your triangulation. To get the Voronoi tessellation: 1) Circulate the faces of each vertex in order (use orientation test if they are not ordered yet) 2) Connect the circumcenters of each 2D face you visit. Voilà!

vanderZwan 8y ago
Well, conceptually maybe, but I have some doubts regarding: > 2) Compute Delaunay in 3D How much slower is 3D triangulation compared to 2D? How does it scale with more points? And just as importantly: how complex does the algorithm get when we want to improve performance? I'm not saying you're wrong - if anything I hope I'm wrong, but I suspect the answers to my questions are "significant", "badly", and "very", which would mean this is not easier/more robust in practice (or it comes at the price of being unacceptably slow). The reason I think so is Julia's Voronoi/Delaunay package[0]. It implements an algorithm from an extremely dense (to non-astronomers) astronomy paper[1]. Even the name is intimidating: "E pur si muove: Galiliean-invariant cosmological hydrodynamical simulations on a moving mesh" Page five and six of the paper lists a short history of triangulation algorithms, making it clear that Delaunay triangulation so important for scientific research that many people have given it considerable thought. Sadly, I cannot really follow the algorithm explained in that paper. It seems to fall for the inevitable "algorithm by PhDs with complex enough math that you need a comparable PhD to being to grok it"-curse (and/or I am just too dumb). Of course, if someone already did the hard work of implementing really fast and robust 3D Delaunay libraries to steal from... ;) (the Julia package only implements the 2D algorithm. Sadly) [0] https://github.com/JuliaGeometry/VoronoiDelaunay.jl https://github.com/JuliaGeometry/VoronoiDelaunay.jl [1] https://arxiv.org/abs/0901.4107 https://arxiv.org/abs/0901.4107
- contravariant 8y ago
  Computing Delaunay on a sphere isn't exactly the easiest either, distances get rather tricky and there's no simple projection that preserves all of them.