Delaunay and Voronoi on a sphere | Hacker News Reader

Delaunay and Voronoi on a sphere

redblobgames.com

250 pointsguiambros35 comments

> 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

[1] https://cs.nyu.edu/~ajsecord/msc_thesis/ajsecord_msc_thesis....

[2] https://beta.observablehq.com/@mbostock/voronoi-stippling

[3] https://beta.observablehq.com/@jobleonard/untitled (Yes, that still needs a lot of optimization. I know)

It's similar, but the purpose is different. With Lloyd Relaxation you move the Voronoi seed points (input), and then regenerate Voronoi. What I'm doing here is moving the Voronoi vertices (output) from the circumcenters to the centroids. There's no iteration here because the input points don't move. Instead, I am constructing a non-Voronoi tesselation. It has the same structure as the Voronoi cells but the cells do not have the Voronoi property. For my map generator I don't need the Voronoi property; I needed other properties, which the centroid approach seems to give me.

> He is basically doing a form of relaxatation similar or equivalent to Lloyd's Algorithm[0], no?

Yes, this is one iteration of Lloyd’s algorithm.

> Would iterating a few times maybe improve results even more?

Yes, it would. Even one more iteration will help, the first steps tend to be the largest, and they usually get smaller and smaller until nearly converged.

One interesting property of Lloyd’s algorithm, though, is that steps don’t get monotonically smaller. If you try to wait for true convergence, sometimes you’ll see large jumps many hundreds of iterations into the simulation. It can be a bit surprising to watch a point jump after all the points have been sitting very still for a while.

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

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

[2]https://projet.liris.cnrs.fr/mi2/posts/2018/10/05/voronoi-pr...

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

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.

Thanks!

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

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

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.

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

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

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.

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.

It might be that the triangle is big enough (relative to, say, the radius of the sphere) that that approximation doesn't apply.

I've used CSSGRID for this http://www.ncarg.ucar.edu/ngmath/cssgrid/contents.html

It's solid and easy to bind from anything with an FFI.

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

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 (click </> to see the code). Probably someone did it before already :)

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

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.

See section 8.5 of 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

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.

Yeah I used regular as a shorthand but nostly I meant that it would flip edges to avoid thin triangles.

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.

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

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à!

The Delaunay triangulation of a set of points on the sphere is given by the convex hull. The Voronoi vertices on the sphere are given by the facet normals and vice versa. The convex hull is quite fast to compute.

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

[1] https://arxiv.org/abs/0901.4107

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.

If I understand correctly, the stereographic projection does handle Delaunay+Voronoi on a sphere, for all points except the south pole. See section 8.5 of http://www.cis.upenn.edu/~cis610/convex8.pdf

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/

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.

Not sure why you need the point at the origin. In fact I think it fails when you e.g. only have points on one half of the sphere.

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

Reply on news.ycombinator.com