
Delaunay and Voronoi on a sphere - guiambros
https://www.redblobgames.com/x/1842-delaunay-voronoi-sphere/
======
vanderZwan
> _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....](https://cs.nyu.edu/~ajsecord/msc_thesis/ajsecord_msc_thesis.pdf)

[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)

~~~
mxfh
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-
pr...](https://projet.liris.cnrs.fr/mi2/posts/2018/10/05/voronoi-projection-
en.html)

~~~
vanderZwan
Oh wow, thanks! I missed those. Will look into it - maybe he has optimizations
I missed and vice versa.

If only Observables had something like pull requests and a way to have meta-
discussions on notebooks... but it's addictive enough an environment as is ;)

EDIT: He even has a Spherical Lloyds Relaxation sketch[0], how appropriate! :D

[0] [https://beta.observablehq.com/@fil/spherical-lloyds-
relaxati...](https://beta.observablehq.com/@fil/spherical-lloyds-relaxation)

------
blt
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
Thanks!

~~~
ChrisClark
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. :(

------
gus_massa
> _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.

~~~
slavik81
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
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
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...

------
tame73_binds
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
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_La...](https://developer.mozilla.org/en-
US/docs/Web/CSS/CSS_Grid_Layout)

------
hamoid
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
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).

------
vbarrielle
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
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
Yeah I used regular as a shorthand but nostly I meant that it would flip edges
to avoid thin triangles.

------
icoder
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.

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

------
rosshemsley
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
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
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.

~~~
vanderZwan
True, but it can be decomposed into different steps at least (trying to find
less-wrong projections, effective 2D triangulation).

The complexity of working with 3D spaces typically feels worse and harder to
untangle to me.

It's like.. I dunno, having two jugglers juggling ten balls each and passing
them between each other, or one being asked to juggle twenty[0]. But that is a
personal thing - it may also just be sign I'm bad at maths in 3D space.

[0]
[http://www.juggling.org/papers/limits/](http://www.juggling.org/papers/limits/)

~~~
contravariant
In the end it's a trade-off between having an extra dimension to deal with, or
having to deal with a non-euclidean space. Between the two it's nearly always
easier to go for the higher dimensional euclidean one, but exceptions might
exist.

Of course if you want an approximate solution then projecting the sphere to 2D
might be an option, but glueing together the different projections required
could be a bit of a hassle.

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

