Voronoi-Delaunay mesh

 from Red Blob Games
DRAFT
15 Jun 2018
  1. Delaunator manual https://mapbox.github.io/delaunator/[1]
  2. Map generator that uses the ideas on this page https://twitter.com/redblobgames/status/1024002191978160130[2]

–

People know how to work with square grids; I have a page about hexagonal grids. Right now I have a map generation page with voronoi/delaunay. But it's voronoi + how to use voronoi for maps + specific algorithms for rotmg maps. Could the non-rotmg topics be a topic as detailed as the hexagon page?

–

{ use two column layout? }

Square grids and hex grids (“tesselations[4]”) are common in games. We can generalize from these to Voronoi tesselations[5]. When the diagram's red points are equally spaced, we get a grid. When we use an irregular spacing, we can get a cool organic look. Voronoi can produce both. On this page I'll describe how to work with Voronoi and related structures.

{ canvas might be better than svg for this page; need to make an inventory of diagrams – do I need hover etc. more often, or do I want larger draw counts? it would also be nice to produce images for search engines }

Given these red points, we can construct a Voronoi diagram. {khan academy animation here?}

{ should be possible to smoothly vary from voronoi to centroid in the same diagram, by using x to decide on the mix }

we lose

However, we retain neighbors and neighbor-related algorithms like pathfinding, and we gain a cool organic look. We can vary the density and shapes of the cells.

I'll start with Voronoi, but will also describe a related system that has even nicer properties for game maps.

{ might be nice for interactive diagram allow moving points around }

{ show diagrams for square, hex, triangle, hybrid, curved, blue noise, and random }

Put a bunch of examples here

Example
Variable density

{ this might look cool, especially monokot_root.jpg , but it's not the primary focus }

http://byronknoll.blogspot.com/2013/02/smooth-voronoi-diagrams.html[6]

https://twitter.com/HumanEarthPig/status/985993311835090944[7]

http://mkmra2.blogspot.com/2016/10/[8] "smooth voronoi"

 2  Dual mesh#

There are two meshes here.

  1. Voronoi: there are the red polygons with blue corners, and white sides
  2. Delaunay: there are the blue triangles with red corners, and black sides

This is a dual mesh, so the central points in one mesh also act as the vertices for the other mesh. The red points are labels for the regions, and vertices for the triangles. The blue points are labels for the triangles and vertices for the regions. This is weird and it will take some time to get used to.

Centroid (Barycentric) ← → Circumcenter (Voronoi)
dual mesh

{ this diagram should not have the centroid/circumcenter slider, but instead a voronoi/delaunay slider }

What should I call the red and blue points? The names “corners” and “vertices” are ambiguous. The word “polygons” also includes triangles. The word “tiles” suggests polygons with the same shape. I've decided to call them regions and triangles.

red points
These are the inputs to Delaunay triangulation, and will get integers 0 ≤ r < numRegions.
delaunay triangles
These are the outputs of Delaunay triangulation, and will get integers 0 ≤ t < numTriangles.
black lines
These are also the outputs of the Delaunay triangulation, and get integers 0 ≤ s < numSides. There are three sides per triangle, so numSides = 3 * numTriangles. They're represented as half-edges, so adjacent triangles have separate sides instead of sharing a side ID.
blue points
There's one corner of a Voronoi polygons for each of the triangles, so they'll get the same IDs.
voronoi regions
There's one region for each of the red points, so they'll get the same IDs.
white lines
There's one side of a Voronoi regions for each side of the triangles (it may not be obvious), so they'll get the same IDs. These too are half-edges, so adjacent regions will have separate sides instead of sharing a side ID.

The Voronoi regions do not require any additional IDs beyond what the Delaunay triangulation uses.

To summarize, there are three types of ids: regions, triangles, and sides. Region IDs and triangle IDs are unrelated.

NameVoronoiDelaunayId
faces/cells/tiles/polygons0 ≤ r < numRegions
vertices/triangles0 ≤ t < numTriangles
edges/borders/sides0 ≤ s < numSides

{ this section could be condensed - seems a little redundant - but it's a key concept so repetition may be useful }

 2.1. Parts

There's no coordinate system for an irregular grid, so instead we can either assign IDs to each part or make objects for each part. I prefer to assign numeric IDs. This lets me store additional data in arrays. For example if I want to store the elevation with each region, I can have a separate array elevation[r] indexed by the region number.

{ are these diagrams still useful? could they be combined into one that shows the labels on mouseover? or can I use the left/right trick from above to show all three in the same diagram? }

2.1.1. Regions

"tiles", "regions", "polygons", "cells", "faces", ("vertices" in dual)

Regions

2.1.2. Corners

"vertices", "corners", ("triangles" in dual)

Corners

2.1.3. Sides

"edges", "sides", "borders"

Sides split into two half-edges

 2.2. Delaunator library

I'll describe the Delaunator library's data structures here. Other libraries will have different representations. We need ways to get from triangles to sides, sides to triangles, triangles to regions, regions to triangles, regions to sides, and sides to regions. Delaunator provides the building blocks, and then I'll build the rest in the next section.

{{ this section is obsolete now that delaunator has a guide, right? }}

2.2.1. Triangle and Side IDs

Delaunator uses side (edge) IDs, and puts them in groups of three for each triangle. This means triangleId = floor(sideId / 3), and sideID = 3*triangleId + {0,1,2}. It also means numSides = 3 * numTriangles.

sideId = triangleId + {0, 1, 2}

This lets us iterate the sides of a triangle. Given a side ID, we can find the triangle it's part of. Given a triangle ID, we can find the three sides of it. Given a side, we can find the next and previous sides in the triangle.

function s_prev_s(s) { return (s % 3 == 0) ? s+2 : s-1; }
function s_next_s(s) { return (s % 3 == 2) ? s-2 : s+1; }

2.2.2.triangles array

Delaunator gives us a triangles[sideId] array, returning the region (vertex) id for the beginning of that side (half edge).

triangles: side id → region id

For example triangles[s1] is r1. (The numbers being the same is a coincidence.) This array lets us get from sides to regions. Since we can already get from triangles to sides, we also can use this to get from triangles to regions.

2.2.3.starts array

For some algorithms I need the inverse of the triangles array. Delaunator doesn't provide this but it provides everything I need to build it, so I build starts[regionId] to tell me any side (half edge) that starts at a given region (vertex).

starts: region id → start id
let starts = new Int32Array(numRegions);
for (let s = 0; s < triangles.length; s++) {
    starts[triangles[s]] = s;
}

For example starts[r2]] can be s14. It could've been s0 or s5 or s9 or s13. It doesn't matter which (unless you don't use ghosts, which I'll explain later). This gets us from regions to sides and triangles.

  1. TODO non-ghost version

    There's a problem on the outer hull of the triangulation. We're going to use the starts array to iterate through the halfedges, but some halfedges may be -1. The solution is rewind each element by iterating backwards until we hit a -1.

    { need to write code for this and test it thoroughly }

    function s_rewind_s(start_s) {
  // TODO: TEST THIS!! I don't think it's right
  let s = null,
      next_s = start_s;
  do {
      s = next_s;
      let opposite_s = halfedges[s];
      next_s = s_prev_s(opposite_s);
  } while (s != start_s && next_s != -1);
  return s;
}

    or maybe

    function s_rewind_s(start_s) {
  let s = start_s;
  while (true) {
    let back_s = halfedges[s_prev_s(s)];
    if (back_s === -1 || back_s === start_s) {
      break;
    }
    s = back_s;
  }
  return s;
}

    and then use

    let starts = new Int32Array(numRegions);
for (let s = 0; s < triangles.length; s++) {
    starts[triangles[s]] = s_rewind_s(s);
}

2.2.4.halfedges array

Delaunator provides a halfedges[sideId] array, returning another sideId for the half edge's “opposite pair”.

halfedges data structure goes from side id to opposite side id

For example halfedges[s2] is s5, and halfedges[s5] is s2. This is how we'll get from one triangle to adjacent triangles.

 3  Traversal algorithms#

Delaunator provides the building blocks. Here I'll describe how I use its data structures for traversing the dual mesh. 9 of these, and my names for them

 3.1. TODO Halfedge properties

Halfedge id to start, end r

Halfedge id to inner, outer t

 3.2. Delaunay triangle circulation

There are two main patterns to look at. Inside a triangle, we can circulate around the sides (halfedges):

circulate inside a triangle

This function gives us the three sides of a triangle, which have consecutive IDs starting at 3 times the triangle ID:

function t_circulate_s(t) {
  let out_s = [];
  for (let i = 0; i < 3; i++) {
    let s = 3*t + i;
    out_s.push(s);
  }
  return out_s;
}

Given a side we can find the region point by looking in the triangles array. This lets us get the three region points that make up a triangle:

function t_circulate_r(t) {
  return t_circulate_s(t).map(s => triangles[s]);
}

To go between triangles, we can jump between a side (halfedge) and its opposite:

halfedges[s2] == s5; halfedges[s5] == s2

This lets us find the adjacent triangles around a given triangle:

find adjacent triangles
function t_circulate_t(t) {
  return t_circulate_s(t).map(s => floor(halfedges[s]/3));
}

{ note that this doesn't work without ghosts }

If we start with a triangle ID we can use this loop to find sides, triangles, or regions.

 3.3. Voronoi region circulation

Another useful way to combine these two operations is to circulate one step, jump to the opposite, circulate one step, jump to the opposite, repeat:

find adjacent polygons

This is tricky, so let's list the steps to get the sides leading out of the center region point:

  1. Find any halfedge leading out from the starting point. One such halfedge is s14. (use the starts array for this)
  2. Find s14's opposite, which is s2. (use the halfedges array)
  3. Circulate around s2's triangle, finding s0. (use the s_next_s helper function)
  4. Find s0's opposite, which is s5.
  5. Circulate around s5's triangle, finding s3.
  6. Find s3's opposite, which is s8.
  7. Circulate around s8's triangle, finding s6.
  8. Find s6's opposite, which is s10.
  9. Circulate around s10's triangle, finding s11.
  10. Find s11's opposite, which is s13.
  11. Circulate around s13's triangle, finding s14.
  12. We've come back to s14, so we're done. The outgoing halfedges were [s14, s0, s3, s6, s11]. The incoming halfedges were [s2, s5, s8, s10, s13].

Here's what the code looks like for finding outgoing sides from a region point:

function r_circulate_s(r) {
  let out_s = [];
  let start_s = starts[r];
  let s = start_s;
  do {
    out_s.push(s);
    let opposite_s = halfedges[s];
    s = s_next_s(opposite_s);
  } while (s != start_s);
  return out_s;
}

The same loop works for finding the triangles touching a region point, or finding the vertices to draw the voronoi polygon:

function r_circulate_t(r) {
  return r_circulate_s(r).map(s => floor(s/3));
}

or finding the neighboring regions adjacent to a given region:

function r_circulate_r(r) {
  return r_circulate_s(r).map(s => triangles[s]);
}

That's it! Those are the two basic loops we need to find our way around the dual mesh. In both cases we focus on the halfedges and then construct the region or triangle data from them.

{ animation might be useful here }

{ to test: if we don't have ghosts, then the loop might work nicely if we stop at halfedges[start_s] >= 0? s_prev_s(start_s) : -1 }

 3.4. TODO Drawing

Show sample code for drawing the delaunay triangles

Show sample code for drawing the voronoi regions

{ diagrams corresponding to sample code }

 3.5. TODO Distances

Have to explore the graph with breadth first search; no coordinates means no fast formula.

{ diagram where you mouseover any region and it shows distances to all other regions }

 4  TODO Generation#

 4.1. Three ways to generate Voronoi

  1. Starting from points, find areas that are closer to one point than any other point.
  2. Starting from Delaunay edges, use perpendicular bisectors to construct new edges.
  3. Starting from triangle areas, find center points using circumcenter.

 4.2. TODO point selection

{ this seems like an independent problem }

4.2.1. grid

4.2.2. jittered grid

4.2.3. random

4.2.4. poisson disc

4.2.5. lloyd relaxation

 4.3. TODO delaunay

I often want to work with the centroid version, not the circumcenters. It's simpler for me to start with Delaunay and then calculate the centroids than for me to use a Voronoi library.

Circumcenters and orthocenters are not guaranteed to be inside the triangle. Centroids and incenters are.

{ should this explain what to look for in a delaunay lib? need the half edge and triangle data }

look for a library that gives you the IDs instead of x,y!

does it matter if all the triangles are oriented the same way? this might be a cute animation, showing not only the circulation order but have that sprite spin when you mouse over a triangle

Some libraries use iterators or arrays so that you can get an integer id.

Some libraries don't give you all the information you need, but you might be able to reconstruct the dual graph.

There are other formats too like http://www.gradientspace.com/tutorials/dmesh3[18]

 4.4. TODO triangle centers

4.4.1. circumcenter (voronoi)

{formula for calculating; link to wikipedia}

4.4.2. centroid (barycentric mesh)

{formula for calculating; link to wikipedia}

{ also link to all the other candidates for triangle centers }

4.4.3. incenter

 4.5. TODO Ghost and solid

Problem: halfedge could be -1. This complicates some of the traversals. {but maybe this is after traversal}

Solution: add "ghost" region on the outside of the mesh. Then add sides and triangles between the ghost region and the solid regions.

Solution: add checks to all the code.

 5  TODO Appendix: Relationship between parts#

squares:

hexagons:

triangles:

voronoi:

So let's look at our diagram again. I said that square grid corners always have 4 neighbors. But I also said voronoi corners always have 3 neighbors. And I also said square grids have the voronoi property.

All three of these can't be true.

Overlap, maybe centroid slider to show this.

 6  TODO Appendix: Centroids and Incenters#

For my map projects I want to put things at the blue points, like roads or castles. When two blue points have the same coordinate, or even are close, things can go wrong.

Let's separate the blue points. We lose the voronoi property and can no longer generate squares. It keeps the same mesh structure and it simplifies other things. More on this later.

I should also try incenter because it seems to have some good properties too (e.g. it's the angle bisector for the wedge of a polygon in a triangle's direction, and the biggest circle inside the triangle is drawn there) and I don't actually know which will be best

Centroid is useful in simulations - if you have f(r) for all the vertices r of a triangle, then the best estimate of f(centroid) is the weighted average of the three f(r) values. And if we build half wedges for polygons then we are bisecting the white edge

Not sure if this will go here, or if it's up near the top.

The square grid is a clear example of the problem. Square grid voronoi produces are square grid. But the triangles tell you connectivity between adjacent polygons. And the triangles tell us that each square has six neighbors, not four!! If we look closely we'll see there are circumcenters that overlap. Two triangles have the same circumcenter. This is messy.

{ diagram with a tiny square grid showing the centroid points joining together when they become a circumcenter }

 7  TODO Appendix: Drawing with bezier curves#

Related: https://sighack.com/post/chaikin-curves[19]

 8  TODO Appendix: Weighted regions#

Instead of calculating the centroid, calculated the weighted center based on weights assigned to the region centers

 9  TODO Appendix: Pixel to polygon#

spatial hash, point-in-poly test

{ diagrams showing what the spatial hash matches, and how we need the point-in-poly to be sure }

 10  TODO Appendix: Decomposing into triangles#

{ diagram }

 11  TODO Appendix: Local support#

In bezier curves we say "local support" means only nearby control points matter. Something similar happens here. Only nearby points affect the triangulation nearby.

{ is this true? need diagram }

 12  TODO Appendix: Subdivision#

Subdivide the triangle, not the regions. Then construct new regions. Not sure if I should have this section.

 13  TODO Appendix: Cellular automata#

https://www.jasss.org/4/4/6.html[20]

https://link.springer.com/article/10.1007/s00158-016-1614-z[21]

https://www.researchgate.net/publication/224644505_Using_Rule-based_Fuzzy_Cognitive_Maps_to_Model_Dynamic_Cell_Behavior_in_Voronoi_Based_Cellular_Automata?_sg=R4LC83_UdXvKMkST3L0JWFP3rgOA_8o-q29eEjo6_8PGshHyoIIEXzE58w32fvd70Ydd1R2NZg[22]

 14  TODO Appendix: Voronoi property#

Voronoi finds all regions closer to a given point than to any other point.

Radius:
how voronoi works

Animation idea from Khan Academy[23] / Peter Collingridge.

 15  TODO Appendix: How to tell if it's Voronoi#

{ diagram }

  1. Take two adjacent (red) points
  2. Draw the black edge between them
  3. The white edge should be along the perpendicular bisector

{ should this be a quiz? }

 16  TODO NOTES#

What to call this?

Twitter suggestions: chickenwire, cellular grid, baked mudflat, a pretty butterfly, foam, bubbles, cellgrid, polygrid, polymesh, convex irregular polygon tiling, tiling, net, voronoi grid, cell structure ngram viewer suggest unstructured grids beats irregular grids, but polygon mesh is also popular.

https://books.google.com/ngrams/graph?content=irregular+grid%2Cunstructured+grid%2Cpolygon+mesh%2Cirregular+grids%2Cunstructured+grids[27]

I decided on "voronoi" because it's what gamedevs will use (especially for procedural maps, which is my main application). https://en.wikipedia.org/wiki/Tessellation#Voronoi_tilings[28]

"Dirichlet tessellations"

"Thiessen polygons"