Fortune's Algorithm is pretty simple, but it's tricky to get right, particularly if you're also maintaining the information in a form that lets you do something useful with it besides spitting out the points of the Voronoi centers.
I agree, Voronoi diagrams are so much fun! For that map project, at the last minute I switched to something slightly different. The page still says "Voronoi" but the actual demo doesn't use Voronoi. The Voronoi diagram uses the circumcenters of the triangles produced from Delaunay triangulation. I found that for my needs, the centroids worked better ("barycentric dual mesh"). Some comparison pictures here: https://www.redblobgames.com/x/1721-voronoi-alternative/ The graph connectivity is the same but the shapes are rounder.
I also love using Voronoi diagrams as well as Delaunay diagrams (via Mathematica).
My most recent exploration with them was to understand and visualize the different structures and characteristics of various low discrepancy quasirandom point distributions in two dimensions.
I also consider Voronoi diagrams fun, and in my life I've implemented three different algorithms for generating 2D Voronoi diagrams (or Delaunay triangulations, which are the duals of Voronoi diagrams):
1. Fortune's algorithm [0]
2. The Bowyer-Watson algorithm for generating the Delaunay triangulation incrementally [1]
3. Quickhull (which generates the 3D convex hull of points, which is the 2D delaunay triangulation of points). [2]
Of these, Quickhull is the simplest, and Bowyer-Watson is by far the hardest. Bowyer-Watson seems very simple, but there are demons hiding in that algorithm, and almost every implementation you can find online is incorrect (for instance, this one [3]). You can tell that these implementations of Bowyer-Watson are incorrect, because every Delaunay triangulation contains the convex hull of the points, and it's easy to tell when they don't. In fact, even the description in Computational Geometry: Algorithms and Applications (a standard textbook) is subtly wrong.
The reason Bowyer-Watson is hard is this: it's an incremental algorithm where you add the points one-by-one to an already constructed Delaunay triangulation, you make sure all edges are flipped right, and at the end you have the full triangulation. However: you have to have a triangle to start with. In order for the algorithm to work, this initial "super-triangle", has to contain all the points in the set, and also be made up of points that are not contained in any circumcircle of any combination of three points in the set. However: if three points are colinear (or very close, which it is almost guaranteed some points are going to be), the circumcircle of those three points is MASSIVE (essentially "infinite", covering the entire half-plane). But the super-triangle points still have to be outside of it. This means that the super-triangle points have to be essentially "symbolic" points (not normal points with coordinates, but special magic points), and you have to hard-code special rules for them. It is extremely difficult to implement these rules correctly.
Fortune's algorithm is on the surface more complex, but there's less goblins hiding in that algorithm. There's some tricky data-structure stuff, but it's not too bad. Quickhull is fairly straight-forward, and is the algorithm I would recommend if you want to give this whole Voronoi adventure a go (also has the benefit of giving you a 3D convex hull algorithm, which you can have all sorts of fun with!).
I spent a couple weekends trying to implement Voronoi via Fortune's algorithm in PHP, since there is no existing library for generating them.
I couldn't get past the point of creating a seemingly random mass of intersecting lines. My math knowledge is lacking (I'm not good at groking geometry for some reason), so I wasn't sure what was wrong or how to approach fixing it.
Maybe I should give Quickhull a go. I think I understand what it's doing.
To be clear: quickhull is not "easy" compared to other things you might program, it's an involved algorithm. But honestly, all Voronoi/Delaunay algorithms are difficult to implement in performant fashion, there's no free lunch.
The hardest part of quickhull in 3D is the "horizon" calculation: basically, every time you add a point to the hull, you have to find all the faces that the point can "see", and the horizon of all those faces, and replace them with new faces. This page describes the horizon part well: http://algolist.manual.ru/maths/geom/convhull/qhull3d.php
When I did it, I basically had to break down to the point of making a WinForms app that I could step through the algorithm iteration by iteration, drawing out the beach-line and and all the circles and edges as it progressed; I had to have that visualization to see that it was correct.
I've not implemented Fortune's myself. I use the Delaunator library, which maintains all the additional information I wanted in the form of a half-edge triangle mesh: https://github.com/mapbox/delaunator ; Mike Bostock wrote a Voronoi wrapper on top it: https://beta.observablehq.com/@mbostock/the-delaunays-dual