He is basically doing a form of relaxatation similar or equivalent to Lloyd's Algorithm, no? Adrian Secord used this way back in 2002 already to create stippling images, see also Mike Bostock's notebook. 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. It's been a lot of fun! A bit too much at times, Observables makes the interactive feedback loop a bit too easy.
 https://beta.observablehq.com/@jobleonard/untitled (Yes, that still needs a lot of optimization. I know)
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.
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, how appropriate! :D
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 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.
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'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.
It's solid and easy to bind from anything with an FFI.
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.
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.
> 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. It implements an algorithm from an extremely dense (to non-astronomers) astronomy paper. 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)
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. But that is a personal thing - it may also just be sign I'm bad at maths in 3D space.
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.