[Programming tools]: Constrained Delaunay Triangulation

Hi everybody!

I'm creating a 2D NavMesh from scratch for my game, with the intention of automatically generate the navigable areas of the level tilemaps. The first step has been the implementation of the triangulation algorithm whose result will be used as the navigation mesh. I want to share with you my implementation for Unity and a tutorial that explains both the constrained Delaunay triangulation algorithm and the implementation. For me, it was a bit hard to code it by just reading the papers of Sloan, so I thought it could be useful for other people to count on an alternative and more detailed explanation. I spent a lot of time writing and drawing this week to make it as easy to understand as possible, please don't hesitate to give me your opinion or ask any question. Hope it helps!

I re-posted this thread in my website:
https://twitter.com/JailbrokenGame/status/1641540424195280897
https://www.jailbrokengame.com/2023/03/30/constrained-delaunay-triangulation/

You can find other of my shared implementations here:

• LineRenderer2D: GPU pixel-perfect 2D line renderer for Unity URP (2D Renderer)
• Script for generating ShadowCaster2Ds for Tilemaps
• Target sorting layers as assets

You can follow me (and my game) on Twitter: @JailbrokenGame

Constrained Delaunay Triangulation

Implementation (Unity / C#): https://github.com/QThund/ConstrainedDelaunayTriangulation
PDF version (better formatted, more images): Download

Documentation references
These are the 2 only papers I have used as a base for my implementation:

• A fast algorithm for constructing Delaunay triangulations in the plane (Advanced Engineering Software, Vol. 9 No. 1, 1987), by S. W. Sloan.
  
  
• A fast algorithm for generating constrained Delaunay triangulations (Computers and Structures Vol. 47 No. 3, 1993), by S. W. Sloan.

Inputs and outputs
The algorithm expects 2 inputs, the main point cloud to triangulate and, optionally, the shape of the holes, non-connected areas of the triangulated surface that must not contain triangles. All the points must belong to the XY plane. The holes are determined by sorted point sets that form closed polygon outlines (constrained edges), with a minimum of 3 edges. It is not necessary that polygons are convex.

The result of the algorithm is a list of triangles that belong to the XY plane and fulfill the Delaunay rule. All the points provided in the main point cloud will be present in the resultant triangles unless a hole prevents them from being there. Points provided as the outline of a hole may be included in the result.

Spatial partitioning optimization
Sloan proposes in his papers the construction of a rectangular grid that divides the space occupied by the point cloud in cells in an early stage of the process, with the aim of optimizing a later step, concretely the search-triangle-by-adjacency step. Every cell contains a numerated “bin” in which we deposit points. What this process achieves is to sort all the points by proximity, all the points that belong to the same region of the space are stored together, and the next cell contains a set of points that are as far as the width of a cell.

This is how it works:

1. Calculate the boundaries of the point cloud.
2. Determine the amount of columns and rows of the grid. Sloans proposes a MxM grid where M is the fourth root of the amount of points in the main point cloud.
3. Reserve a memory buffer to store the M^2 cells.
4. Every time a point P is added to the grid, calculate the row and column that corresponds to its position relative to the bottom-left corner (0, 0).
   
   Row = 0.99 x M x Py / GridHeight
   Column = 0.99 x M x Px / GridWidth
   
   
5. According to the row and column of the cell that contains the point, calculate the bin number B.
   
   B = Row x M + Column + 1, when Row is even
   B = (Row + 1) x M - Column, when Row is odd
   
   
6. Add the point to the bin (a point list).

As you will read later, when searching for the triangle that contains a point, we use the last created triangle as the first triangle to check; in case that triangle does not contain the new point, we jump to its adjacent; since we will be adding the points to the triangulation in the same order they are stored in the grid, we know that the next point will be near to the last created triangle.

The class PointBinGrid implements this mechanism. It stores an array of point lists (bins).

The triangle set
The algorithm requires to keep track of the vertices in the triangulation, the triangles formed by those vertices and the spatial relation between the triangles, i. e. what triangles are adjacent. All that data is managed by the same entity, the DelaunayTriangleSet class. It also provides some functions to obtain useful processed data, like searching for the triangles that share an edge or calculating which triangles intersect a line. Despite its name, this class knows nothing about the Delaunay condition or triangulation algorithms.

There are 3 arrays in this entity: adjacent triangles (indices), triangle vertices (indices) and points (2D vectors). There is no “array of triangles”, triangles are defined by groups of 3 consecutive indices; for example, the first triangle is defined by the vertices in positions 0, 1 and 2 of the array, whereas the second triangle is defined by the vertices in the positions 3, 4 and 5. So, when we want to access the data of the triangle with number T (starting at zero), we read the array elements T x 3 + 0, T x 3 + 1 and T x 3 + 2.

• Adjacent triangles: It stores the 3 adjacent triangles of a triangle, by their triangle index, in counter-clockwise (CCW) order. The adjacent triangle at the position 0 shares the edge which starts at the vertex 0 and ends at the vertex 1.
• Triangle vertices: It stores the indices of the points in the Points array that form every triangle, in CCW order.
• Points: The actual position of the points in the space, in the order they were added. A point may belong to many triangles. Points are never removed from the triangulation.

A very common operation in this implementation is wrapping around when iterating through the 3 vertices of a triangle. Imagine we are processing the second vertex V1 (1) of a triangle (0, 1, 2), and you want to modify the other two vertices of the triangle; in order to make the code as generic as possible and avoid IF blocks, you can just add 1 and 2 to the known index, and apply the modulo 3 operation afterwards; in this case, V2 = (V1 + 1) % 3 = 2, V0 = (V1 + 2) % 3 = 0.

Data structures
In order to make the code more readable and easy to work with, 2 data structures have been defined: DelaunayTriangle and DelaunayTriangleEdge. Both are aimed to be temporary (not persistent) data vehicles when moving packed triangle data among functions. The first contains all the data related to a triangle, which are its vertices (indices) and which are its adjacent triangles; the second represents a single edge of a triangle.



Algorithm

1. Normalization
This step is not required but, as Sloan points out, it reduces the lack of precision when rounding off and allows to skip some checks, as we can assume all points will be always between the coordinates [0, 0] and [1, 1].

First, we calculate the bounds of the main point cloud by iterating through all of them and getting the minimum and maximum values for X and Y. Then pick the maximum distance between the height and the width of the bounding box Dmax.

Height = Ymax - Ymin
Width = Xmax - Xmin
Dmax = Max(Height, Width)​

Finally, for each input point, move it in such a way that it keeps its relative distance to (Xmin , Ymin) but as if such bottom-left corner was displaced to [0, 0]; then scale down the point by the maximum size of the bounding box Dmax.

Xn = (X - Xmin) / Dmax
Yn = (Y - Ymin) / Dmax​

2. Addition of points to the space partitioning grid
The normalized points are sorted by relative distance using the space partitioning grid. We can add them in any order. Once the points of the grid are added to the Triangle set (in the following steps) this grid will become useless.

3. Supertriangle initialization
The first triangle in the Triangle set will be a “supertriangle”. We can choose any 3 points as long as the triangle contains all the points. It is preferable to use a very big triangle, even if a smaller one could fit, since that will assure that constrained edges added in the final steps will lay inside of the triangle. Using a supertriangle is convenient as an initialization step for the triangulation process, it avoids the need of performing some initial checks and makes sure the first point we add belongs to a known triangle. I recommend sticking to the triangle proposed in the papers ([-100, -100], [100, -100], [0, 100]). Remember that its vertices must be sorted CCW. This will create 3 points, 3 triangle vertices (0, 1 and 2) and 3 adjacent triangles (#, # and #) in the Triangle set. I will use the symbol ‘#’ to represent “no adjacent triangle” or “undefined” in this document, although in code I use constants with the value -1.

4. Adding points to the Triangle set and Triangulation
The triangulation algorithm operates with one input point at a time and modifies part of the already stored triangle data. It returns the index of the point in the Points array.

4.1. Check point existence
When attempting to add a new point, we first check whether that point already exists; if it does, we just return the index of the position where the point is.

4.2. Search containing triangle
If it is really new, we search for a triangle in the Triangle set that contains it. The first triangle we will check is the last created triangle. For the first point ever added to the set, the first triangle to check will be, obviously, the supertriangle.

For each edge in the triangle, we calculate on which side of the line the point is (in the XY plane). Being the edge formed by vertices A and B, you can use the cross-product between AB and AP. Since vertices are sorted CCW, if the point is on the right side of any of the edges then we are sure that it is not contained in the triangle.

The next triangle to check is the one that shares the same edge we just checked. Once we find a triangle in which the point is on the left side (edge included) of the 3 edges, we stop searching.

4.3. Store the point
Just add it to the Points array of the Triangle set and get its index.

4.4. Create 2 triangles
We join the 3 vertices of the containing triangle to the new point, forming 3 new sub-triangles. We have information about where all the vertices are and which triangles are adjacent to which others. We can do some assumptions too, for example, if we use the new point as the first vertex of every triangle then we can be sure about which is the opposite edge (index 1) and hence which is its adjacent triangle outside of the containing triangle. For the new triangles T, and T2, contained in T0:

T1.vertices = (P, T0v[0], T0v[1])
T1.adjacents = (T2 , T0a[0], T0 )
T2.vertices = (P, T0v[2], T0v[0])
T2.adjacents = (T0 , T0a[2], T1 )
4.5. Transform containing triangle into the third
The containing triangle T0 is not removed but transformed into the third sub-triangle. Only 3 indices require modification:

T0.vertices = (P, T0v[1], T0v[2])
T0.adjacents = (T1 , T0a[1], T2 )

Note: I’m using the triangle letters, like T, alternatively as triangle indices, triangle data and groups of 3 points, I hope it does not lead to confusion.

4.6. Add new triangles to a stack
In the next step we are going to use either 2 parallel stacks or 1 stack that contains 2 indices per item: 1) the index of the adjacent triangles that are pending to be checked and 2) the local index (0 to 2) of the edge of those triangles that was shared with a previously processed triangle at the moment they were added to the stack.
The 3 new triangles are added to the stack for a later process step, if they have an adjacent triangle that is opposite to the point P (i. e. if Tna[1] != #). It is assumed that the index of the edge shared with the adjacent triangle, in the three cases, equals 1.

Note: There is a discrepancy here among the 2 documents written by Sloan. One says the 3 triangles must be added whereas the other says only their adjacent triangles (which oppose to P) must.

4.7. Check Delaunay constraint
From Sloan’s papers: “The Delaunay constraint stipulates that the diagonal V1-V2 is replaced by the diagonal P-V3 if P lies inside the circumcircle for the triangle V1-V2-V3”. This is something we have to check every time a new triangle is created.

In the previous step we filled the stack with the new triangles. Triangles are processed one at a time. Let’s call the current triangle Ti . We know that the new point is at Tiv[0] because we forced it to be. That is going to be our point P.

We obtain the adjacent triangle TA that opposes P by reading the value of Tia[X], where X equals the index of the edge that is stored along with the index of the triangle in the stack.

Then V1, V2 and V3 in the previous wording are TAv[0], TAv[1] and TAv[2], respectively. If P is in the circumcircle of TA , the diagonal of the quadrilateral formed by both triangles, Ti and TA, has to be swapped; otherwise, we go back to the beginning of this step and take the next triangle from the stack. Please note that, even after the swap occurs, and no matters how far the swap propagates, the point P will always be Ti [0].

Before swapping the edge, we need to know which of the edges of the adjacent triangle TA coincides with it. We know the index of the current triangle so it is as simple as iterating through the adjacent triangle indices of TA until they match. That gives us the edge number EA (whose value may go from 0 to 2). We use the edge number to know the position of the other 2 adjacent triangles of TA and we add their indices to the triangle stack, propagating the Delaunay constraint check to the next adjacent triangles until no edge needs to be swapped.

Adjacent triangle 1 = TAa[(EA + 1) % 3]
Adjacent triangle 2 = TAa[(EA + 2) % 3]​

4.8. Swap edges
Basically, if we have a convex quadrilateral formed by 2 triangles T0(A, B , C) and T1(C, B, D), this process will generate 2 alternative triangles T0(A, D, C) and T1(A, B, D). We just have to perform certain changes in the vertex and adjacency data of both triangles and their adjacents, without reserving nor releasing memory. Technically we are only moving one vertex in each triangle. Before showing the operations, let’s define the local index (0 to 2) of the vertex of TA that does not belong to the shared edge EA:

VA = TA[(EA + 2) % 3]​

And the same for the index of the opposite vertex in Ti , which we will call Vi . These are the changes to make in all the implied triangles, in order:

Tiv[(Vi + 1) % 3)] = TAv[VA ]
TAv[EA ] = Tiv[Vi ]
TAa[EA ] = Tia[Vi]
T,a[Vi ] = TA
Tia[(Vi + 1) % 3] = TAa[VA ]
TAa[VA ] = Ti​

For the 2 adjacent triangles whose neighbor has changed, we need to know the index of the adjacent triangle that matched their original neighbors, either Ti or TA, and replace it with the other. The index of the triangles to update are:

TB = Tia[(Vi + 1) % 3] → Find TBa[n] that matches TA and replace with Ti
TC = TAa[EA] → Find TCa[m] that matches Ti and replace with TA​

5. Holes creation
A hole is a closed polygon (it is not required that it is convex) described by an array of points H that represent the edges of the outline, sorted CCW. So H[0] → H[1] defines an edge, H[1] → H[2] defines the next edge, and so on. Obviously the last edge is defined by H[Hcount - 1] → H[0]. These edges, also known as constrained edges, will be added to the existing triangulation creating new triangles that must satisfy the Delaunay constraint but avoiding to swap the edges that belong to the polygon outline. If there are 2 consecutive coincident points (zero-length edge) in the outline, one of them must be discarded. For each hole, the points of its outline are normalized and added to the triangulation.

5.1. Normalize
The points of the polygon have to be normalized in the same way the main point cloud was, using the same bounds we calculated for those points, so the polygon is moved and scaled accordingly.

5.2. Add the points to the Triangle set
All the points of the outline are added to the Triangle set, one by one, in exactly the same way we did with the main point cloud. The index of every added point has to be stored in an array so we know which point is connected to which other. There is one array per hole. As we add the points, new triangles are generated which may or may not have the edges of the polygon.

5.3. Create the constrained edges
For each hole, we iterate through its vertex indices, defining the constrained edges with pairs of vertices. First, check if the edge already exists, i. e. whether there is a triangle that has the current edge. If so, ignore the edge; otherwise we can proceed to create it.

5.3.1. Search for the triangle that contains the beginning of the new edge
In a later step we will calculate which triangle edges intersect the constrained edge but, in order to do that in an optimum way, we need to know which triangle to check the first and follow the trajectory of the line; otherwise we would have to calculate the intersection for every edge in the Triangle set.

The edge starts at an existing vertex ESTART so we search for the index of the vertex in the triangle vertices array of the Triangle set and store the indices of the triangles that have that vertex. Then, for each triangle (V0, V1, V2), we identify at which position the vertex Vi is (0, 1 or 2) and, using that value, we obtain the contiguous edges:

Vi1 = (Vi + 1) % 3
Vi2 = (Vi + 2) % 3
E1 = Vi → Vi1
E2 = Vi2 → Vi​

Finally we check if the other endpoint of the edge, EEND , is on the left side of both edges E1 and E2, and we stop searching if that’s the case. We have found which of all the triangles that share the same vertex contains the first endpoint, ESTART , of the future constrained edge.

5.3.2. Get all the triangle edges intersected by the constrained edge
For each triangle TN in the Triangle set, and starting with the triangle that contains ESTART , we first check that TN is not the triangle that contains the other endpoint EEND. Take into account that, normally, the line will intersect 2 edges of every triangle and that, for every triangle, only the “exit hole” is calculated; the “entry hole”, or first intersection, is calculated by the previous checked adjacent triangle, and so on. So the last intersected triangle is a special case, we do not need to check for intersections, it is obvious that the “entry hole” is somewhere and the “exit hole” is its vertex. As you may foresee, as soon as this special case occurs, the process to find intersecting edges finishes and the output list is returned.

For each of the 3 edges in TN we check if EEND is on the right side, which means the line may intersect that edge. If we find that case, we take note of the index of the current edge Ei (0, 1 or 2) and store the value in another place, let’s call it ET; this is the tentative edge of the next adjacent triangle. Now we calculate the intersection among the triangle edge and the constrained edge.

If they intersect:

• Store the index of the vertices that form the edge (VA and VB) in the output list.
• Use a flag to know when any of the edges of TN has been crossed. Set it to True.
• Set the next triangle to check
  
  TN = TN a[Ei ]

If they do not intersect, continue iterating through the edges of TN.

If none of the 3 edges of TN intersect with the constrained edge (check the flag), then

TN = TN a[ET ]​

This process is similar to the search for a triangle that contains a point, we jump from one triangle to its adjacent in the direction of EEND, instead of testing all triangles by brute force.

5.3.3. Form quadrilaterals and swap intersected edges
Once we know all the existing edges that obstruct the constrained edge, we iterate through that list with the purpose of removing them until the constrained edge can join ESTART and EEND, keeping everything triangulated.

We copy the last element (the 2 indices of the vertices that form the intersected edge, VA and VB) and remove it from the list. We use them to search by brute force for a triangle TN in the Triangle set that has that edge, in that order, and also what is the local index Ei (0, 1 or 2) of the edge in such triangle. We cannot store the triangle index of previous steps to save the cost of this search because triangles may change as we swap edges; but vertices stay the same.

We get the opposite adjacent triangle TA that shares the same intersected edge, VA → VB .

TA = TN a[Ei ]​

Then we calculate the local index (0, 1 or 2) of the same edge in TA, by iterating linearly and comparing with VB (remember that the same edge is flipped in the adjacent triangle), and we also get the vertex VP of TA that opposes that edge. If the mentioned comparison is true for the index J, then we can calculate:

EA = J
VP = TAv[(J + 2) % 3]​

We now have almost all the information we need for swapping the edge, in case it is necessary. The condition for that is that the quadrilateral formed by the vertices TNv[0], TNv[1], TNv[2] and VP is convex. If so, we only need one more data, the vertex VN that opposes the edge in TN:

VN = TN [(Ei + 2) % 3]​

We swap the edge VA → VB in the same way we did in previous steps. In case the quadrilateral was not convex, the pair VA-VB is added back to the first position of the list and the current iteration ends.

The last part of this step is to check whether the swapped edge still intersects the constrained edge. We can obtain the endpoints of the edge, VC→ VD easily using a trick. A known consequence of the swap operation is that the local position of the shared edge of the first triangle, in this case TN, is moved 2 positions forward:

Ei’ = (Ei + 2) % 3​

Therefore, the new endpoints are:

VC = TNv[Ei’ ]
VD = TNv[(Ei’ + 1) % 3]​

Be careful with not using a cached version of the information related to TN, it has changed after the swap operation. Now we calculate the intersection between ESTART → EEND and VC → VD. In order to avoid false positives, if the edges share any vertex then we consider them not intersecting. In case they do intersect, the pair of indices of the vertices VC and VD is inserted at the first position of the list of intersected edges; otherwise, we add that information to a list of new edges. The current iteration ends.

5.3.4. Check Delaunay constraint and swap edges
When the algorithm reaches this step, the constrained edge has been already added to the triangulation and there is no other edge obstructing it. For all the new triangles we created, it is necessary to check whether they still fulfill the Delaunay constraint. We iterate through every new edge of the list.

We do not want the constrained edge to be touched anymore so the first thing to do is to discard that the current edge V1 → V2 coincides with it.

As we did in a previous step, we search for a triangle TC that uses V1 → V2 and deduce the data of the adjacent triangle TA. We check if the point TCv[(V1 + 2) % 3] (the one that is not in the edge) is in the circumcircle formed by TAv[0], TAv[1] and TAv[2]. If so, we calculate the local index (0, 1 or 2) of the shared edge in TA and swap that edge; otherwise, just ignore it. The iteration finishes and we evaluate the next edge in the list.

5.4. Identify all the triangles in the polygon
This step can be split in 2 stages: identifying the triangles of the polygon outline and identifying all the triangles enclosed by them. The input L for this step is the list of vertex indices of the outline we generated at the beginning of this section. Since we order the triangle vertices CCW, and we added the vertices of the new edges in the same order of the list, we can be sure that, for each edge EP, the triangle that is inside the polygon is the one that has the 2 vertices in the order they are stored in the input list. This process will generate an output list of triangle indices.

For each triangle TN found in the mentioned way, we first check that it does not form a corner (2 consecutive edges in the outline) by comparing it to the last triangle added to the output list and also to the first one, since the polygon is closed (first edge and last edge are connected).

Then we get the 2 adjacent triangles TA of TN that do not share the current edge EP and check whether such triangles have any of the edges in L that are contiguous to EP. If K is the index of the first vertex in EP, the vertices of such 2 contiguous edges are:

PrevA = L[(K + LCOUNT - 1) % LCOUNT ]
PrevB = L[K]
NextA = L[(K + 1) % LCOUNT ]
NextB = L[(K + 2) % LCOUNT ]​

Remember that the polygon is closed so we need to wrap around L using the modulo operation. If J is the iterator (from 0 to 2) of the vertices of each triangle, then we check the following comparisons:

TAv[J] == PrevA AND TAv[(J + 1) % 3] == PrevB
OR
TAv[J] == PrevB AND TAv[(J + 1) % 3] == PrevA
OR
TAv[J] == NextA AND TAv[(J + 1) % 3] == NextB
OR
TAv[J] == NextB AND TAv[(J + 1) % 3] == NextA​

Take into account that we are comparing the vertices of the edges of TA in the opposite direction too, regarding the outline vertices. When 2 triangles are adjacent, the shared edge V1 → V2 in one triangle is V2 → V1 in the other. This occurs when the triangle is a corner and one of its adjacent triangles is outside of the polygon.

If TA does not have an edge in the outline and it has not been added to the output list yet, then we add it to a stack of adjacent triangles that we can call S, to be processed later. Consider the previous comparisons as an optimization that avoids iterating through the output list always; instead, we only do it if TA is not in the outline.

For the second stage, we propagate the quality of “being inside of the polygon” by adjacency until all the triangles that are not at the outline have been processed. While there are elements in S, we take a triangle out and check if it is already added to the output list. If that is the case, we ignore it and take the next one; otherwise, for each of its adjacent triangles we check whether it exists and whether it has not been added to the output list yet. If those conditions are met, we add it to S. Once the 3 adjacent triangles have been checked, we add the current triangle to the output list.

When this step finishes, we will have a list of indices of all the triangles that lay inside the polygon. We can use this list to remove or skip such triangles in a later step, leaving holes in our triangulation.

6. Supertriangle removal
The process ends when we remove the excess triangles that surround the main point cloud. Those triangles exist because we used a supertriangle as the first triangle of the set. It is easy to identify them, they all share at least one of the vertices of the supertriangle which, as you may deduce, are located at positions 0, 1 and 2 of the Points array. We just collect them by iterating through the triangle vertex array and store them in the list of triangles to remove.

7. Denormalization
If you normalized the input points, it is time to denormalize them. Use exactly the same bounding box calculated in the beginning for the main point cloud. For each point [X, Y] we get [XD, YD]:

XD = X * Dmax + Xmin
YD = Y * Dmax + Ymin​

8. Output filtering
We have the full triangulation, L, on one hand and the list of triangles to be removed, R, on the other. For each triangle in L, we check if it exists in R and, if it does, we remove it from R and check the next triangle; otherwise we add the triangle to the output list.

Known limitations
Although this algorithm works for the majority of the use cases, it has some limitations that are worth to mention:

• The polygons that represent the holes cannot overlap each other. If you need some figures to intersect in order to form a figure, you should merge both polygons into one and use that “superpolygon” as a hole.

• If you use holes that are partially or totally out of the main point cloud, the part of the hole that lays outside may generate undesired new triangles that do not share any of the main points.
  
  
• If your points or holes are placed in such a way that very small triangles are formed, that may produce precision loss, which may lead to exceptions or freezing. You can mitigate this by not scaling the points by DMAX during normalization and denormalization steps.
  
  
• If points are added right on an existing edge, or if a constrained edge crosses exactly though an existing vertex, the process may fail. In order to avoid this, it is recommended to apply a tiny random offset to every point added to the triangulation, specially if constrained edges form figures with many orthogonal corners.

 

Huh, you did some nice extensive description and explanations here. Well done

 

Thank you!

 

Great work.

 

Very cool, and very nice and useful of you to write up this little tutorial and share the code. Such work is a great boon to the community!

 

Thank you, you are so kind.

 

UPDATE: Added the following paragraph to the Known limitations:

If your points or holes are placed in such a way that very small triangles are formed, that may produce precision loss, which may lead to exceptions or freezing. You can mitigate this by not scaling the points by DMAX during normalization and denormalization steps.

 

New updates in the code repository:
-Added optional tesselation (CPU).
-Fixed an important bug related to the calculation of the right side of a vector that was causing a lack of precision and leading to infinite loops in some specific cases.

 

Detected a bug in the Mr Sloan's algorithm, as far as I can see. It does not take into account 2 corner cases:
1-Inserting a point in an existing edge.
2-Adding a constrained edge that intersects with a vertex, instead of an edge.

Will propose a solution soon.

 

After several days spent on diagnosing the problem and applying different patches which improved the result but were not enough, I've come to an unorthodox solution that works like a charm: add random offsets to the constrained edge vertices before adding them, as small as possible (not noticeable for a human eye) but big enough to make a difference, avoiding the 2 corner cases exposed above. It's not a perfect solution since it's not deterministic, but it reduces ENOURMOUSLY the chance of a process failure, which was already a tiny possibility, and I have no more time to dedicate to this.

 

Main post updated.

 

In case you want to see the algorithm in action, here you have some examples I created for my NavMesh2D:

https://twitter.com/SiliconHeartDev/status/1370883076772003840

 

Do you plan to make it multi threaded so it creates meshes off the main thread? That would be awesome.

Been reading the documentation confused why the grid math requires you multiply by 0.99? It didn't explain that number, my guess is to prevent it being exactly on the value of 1 in the grid bounds?

 

Hi, I'm not planning to keep improving this project, I'm using it as it is and it's enough for me.
The 0.99 is explained in Sloan's papers:
"The factor of 0.99 is required to ensure that points with the maximum coordinates do not fall outside the grid"

 

This is pure gold, man! Thanks for sharing this! It's exactly the optimization for our game I was looking for!

 

Thank you! I have to warn you though, one person supposedly found a bug in my implementation so the generated triangles do not meet 100% the rules of the original algorithm, although most of the time they will.
I've seen that Unity devs are going to add an optimization for the collision meshes that uses Delaunay triangulation too, maybe you can use it instead (I guess they tested it more than me), although I don't think they allow holes in the mesh, don't know.

Note: I'm still using my implementation in my game and haven't had any problem so far.

 

That's great, thank you for clarifying it

 

Code (CSharp):

1. Vector2[] vertices2D = new Vector2[] {
2.             new Vector2(0,0),
3.             new Vector2(0,1),
4.             new Vector2(1,1),
5.             new Vector2(1,2),
6.             new Vector2(0,2),
7.             new Vector2(0,3),
8.  
9.             new Vector2(3,2),
10.             new Vector2(2,2),
11.             new Vector2(2,1),
12.             new Vector2(3,1),
13.             new Vector2(3,0)
14.     };

Hi, how can I triangulate concave vertices saved in a vector2 field into a 2d mesh?
Or is it possible to draw a polygon simply by an editor handle?

 

Any chance you can extract the triangulation algorithm from colliders specifically and attach it also to the mesh api so we can also generate actual planes of mesh data without the need for 2D colliders stuck to the x:y plane?

On the mesh api you can have all the same methods but also accept an input for a normal, so for example if the user inputs Vector.Up normal then the mesh will triangulate on the X:Z plane etc.

Also is the algorithm run on a thread?

 

Thank you so much for all of this useful information. I have used it to make my own navmesh.

 

I don't follow how Libtess2 is strictly limited to a 2DCollider though ? But any chance you can send a memo to the teams working on mesh API tools to add a deluanay/triangulation functionality as well? Hopefully they could apply jobs/burst and all that goodness to it as well - then its not strictly tied to 2D and physics related stuff.

It does make more sense for it to belong to the mesh system than a 2D collider. I only found it by chance otherwise I was about to attempt to roll my own which was daunting to say the least.

 

Oh theres no communication between teams thats a shame I assumed devs at Unity can pass messages along, since this thread is about mesh triangulation ultimately. More so than physics. But Unity really only offers triangulation currently via the 2D physics api so there is untapped potential by extracting it from physics and having it as a generalised api - that of course the physics api can still link to anyway.

 

Haha no problem! It's been a very interesting discussion, you provided quite useful information.

 

Could you make a tutorial on how to use this in an existing project

 

Unfortunately I don't time for that now, but anyway I would suggest to investigate the way Unity generates meshes in latest versions, as stated in previous comments. See: https://docs.unity3d.com/2022.1/Documentation/ScriptReference/30_search.html?q=delaunay

 

A contributor called picolino0 has detected and fixed an important bug in the triangulation algorithm that was causing wrong triangles in some edge cases and Unity freezing in others.
https://github.com/QThund/ConstrainedDelaunayTriangulation/pull/3

 

Added a small performance improvement. Available in the code repository.

 

I re-posted this thread in my website:
https://twitter.com/JailbrokenGame/status/1641540424195280897
https://www.jailbrokengame.com/2023/03/30/constrained-delaunay-triangulation/

 

2D Delaunay Triangulation in single C# script (Bowyer–Watson algorithm):

Code (CSharp):

1. // https://en.wikipedia.org/wiki/Bowyer%E2%80%93Watson_algorithm
2. using System.Collections;
3. using System.Collections.Generic;
4. using UnityEngine;
5. public class DelaunayTriangulation : MonoBehaviour
6. {
7.     public int VertexCount = 128;
8.     Material _Material;
9.     List<Triangle> _Triangles = new List<Triangle>();
10.  
11.     struct Edge
12.     {
13.         public Vector2 a;
14.         public Vector2 b;
15.    
16.         public Edge(Vector2 a, Vector2 b)
17.         {
18.             this.a = a;
19.             this.b = b;
20.         }
21.     }
22.  
23.     struct Triangle
24.     {
25.         public Vector2 a;
26.         public Vector2 b;
27.         public Vector2 c;
28.    
29.         public Triangle(Vector2 a, Vector2 b, Vector2 c)
30.         {
31.             this.a = a;
32.             this.b = b;
33.             this.c = c;
34.         }
35.    
36.         public Vector2 this[int i]
37.         {
38.             get
39.             {
40.                 switch (i)
41.                 {
42.                     case 0:
43.                         return a;
44.                     case 1:
45.                         return b;
46.                     case 2:
47.                         return c;
48.                     default:
49.                         throw new System.ArgumentOutOfRangeException();
50.                 }
51.             }
52.         }
53.     }
54.  
55.     bool IsEdgesEqual(Edge p, Edge q)
56.     {
57.         bool x = Mathf.Abs(p.a.x * p.b.x - q.a.x * q.b.x) < 0.00001f;
58.         bool y = Mathf.Abs(p.a.y * p.b.y - q.a.y * q.b.y) < 0.00001f;
59.         return x && y;
60.     }
61.  
62.     bool IsTrianglesEqual(Triangle p, Triangle q)
63.     {
64.         bool x = Mathf.Abs(p.a.x * p.b.x * p.c.x - q.a.x * q.b.x * q.c.x) < 0.00001f;
65.         bool y = Mathf.Abs(p.a.y * p.b.y * p.c.y - q.a.y * q.b.y * q.c.y) < 0.00001f;
66.         return x && y;
67.     }
68.  
69.     bool IsPointInCircumcircleOfTriangle (Vector2 p, Triangle t)
70.     {
71.         float ax = t.b.x - t.a.x;
72.         float ay = t.b.y - t.a.y;
73.         float bx = t.c.x - t.a.x;
74.         float by = t.c.y - t.a.y;
75.         float m = t.b.x * t.b.x - t.a.x * t.a.x + t.b.y * t.b.y - t.a.y * t.a.y;
76.         float u = t.c.x * t.c.x - t.a.x * t.a.x + t.c.y * t.c.y - t.a.y * t.a.y;
77.         float s = 1.0f / (2.0f * (ax * by - ay * bx));
78.         float cx = ((t.c.y - t.a.y) * m + (t.a.y - t.b.y) * u) * s;
79.         float cy = ((t.a.x - t.c.x) * m + (t.b.x - t.a.x) * u) * s;
80.         float dx = t.a.x - cx;
81.         float dy = t.a.y - cy;
82.         float radius = dx * dx + dy * dy;
83.         float distance = (cx - p.x) * (cx - p.x) + (cy - p.y) * (cy - p.y);
84.         return (distance - radius) <= 0.00001f;  // removed square root
85.     }
86.  
87.     Triangle CreateSuperTriangle(List<Vector2> points, out Vector4 bounds)
88.     {
89.         float minX = 1e9f;
90.         float minY = 1e9f;
91.         float maxX = 1e-9f;
92.         float maxY = 1e-9f;
93.         for (int i = 0; i < points.Count; i++)
94.         {
95.             Vector2 p = points[i];
96.             minX = Mathf.Min(minX, p.x);
97.             maxX = Mathf.Max(maxX, p.x);
98.             minY = Mathf.Min(minY, p.y);
99.             maxY = Mathf.Max(maxY, p.y);
100.         }
101.         bounds = new Vector4(minX, minY, maxX, maxY);
102.         float dmax = Mathf.Max(maxX - minX, maxY - minY);
103.         float xmid = (minX + maxX) * 0.5f;
104.         float ymid = (minY + maxY) * 0.5f;
105.         Vector2 p1 = new Vector2(xmid - 20f * dmax, ymid - dmax);
106.         Vector2 p2 = new Vector2(xmid, ymid + 20f * dmax);
107.         Vector2 p3 = new Vector2(xmid + 20f * dmax, ymid - dmax);
108.         return new Triangle(p1, p2, p3);
109.     }
110.  
111.     List<Triangle> Triangulation(List<Vector2> pointList)
112.     {   // pointList is a set of coordinates defining the points to be triangulated
113.         List<Triangle> triangulation = new List<Triangle>();
114.         Triangle superTriangle = CreateSuperTriangle(pointList, out Vector4 bounds);
115.         triangulation.Add(superTriangle);
116.         for (int i = 0; i < pointList.Count; i++)
117.         {   // add all the points one at a time to the triangulation
118.             Vector2 p = pointList[i];
119.             List<Triangle> badTriangles = new List<Triangle>();
120.             for (int j = 0; j < triangulation.Count; j++)
121.             {   // first find all the triangles that are no longer valid due to the insertion
122.                 Triangle t = triangulation[j];
123.                 if (IsPointInCircumcircleOfTriangle(p, t))
124.                 {
125.                     badTriangles.Add(t);
126.                 }
127.             }
128.             List<Edge> polygon = new List<Edge>();
129.             for (int j = 0; j < badTriangles.Count; j++)
130.             {   // find the boundary of the polygonal hole
131.                 Triangle t = badTriangles[j];
132.                 for (int k = 0; k < 3; k++)
133.                 {
134.                     Edge e = new Edge(t[k], t[(k + 1) % 3]);
135.                     List<Edge> current = new List<Edge>();
136.                     for (int m = 0; m < badTriangles.Count; m++)
137.                     {
138.                         if (m == j) continue;
139.                         Triangle tr = badTriangles[m];
140.                         for (int q = 0; q < 3; q++)
141.                         current.Add(new Edge(tr[q], tr[(q + 1) % 3]));
142.                     }
143.                     bool contains = false;
144.                     for (int n = 0; n < current.Count; n++)
145.                     {
146.                         if (IsEdgesEqual(current[n], e)) contains = true;
147.                     }
148.                     if (contains == false)
149.                     {
150.                         polygon.Add(e);
151.                     }
152.                 }
153.             }
154.             for (int j = 0; j < badTriangles.Count; j++)
155.             {   // remove them from the data structure
156.                 for(int m = triangulation.Count - 1; m > -1; m--)
157.                 {
158.                     if(IsTrianglesEqual(triangulation[m], badTriangles[j]))
159.                     {
160.                         triangulation.RemoveAt(m);
161.                     }
162.                 }
163.             }
164.             for (int j = 0; j < polygon.Count; j++)
165.             {   // re-triangulate the polygonal hole
166.                 Edge e = polygon[j];
167.                 Triangle t = new Triangle(p, e.a, e.b);
168.                 triangulation.Add(t);
169.             }
170.         }
171.         for (int i = triangulation.Count - 1; i > -1; i--)
172.         {   // done inserting points, now clean up
173.             bool ax = triangulation[i].a.x < bounds.x || triangulation[i].a.x > bounds.z;
174.             bool ay = triangulation[i].a.y < bounds.y || triangulation[i].a.y > bounds.w;
175.             bool bx = triangulation[i].b.x < bounds.x || triangulation[i].b.x > bounds.z;
176.             bool by = triangulation[i].b.y < bounds.y || triangulation[i].b.y > bounds.w;
177.             bool cx = triangulation[i].c.x < bounds.x || triangulation[i].c.x > bounds.z;
178.             bool cy = triangulation[i].c.y < bounds.y || triangulation[i].c.y > bounds.w;
179.             if (ax || ay || bx || by || cx || cy) triangulation.RemoveAt(i);
180.         }
181.         return triangulation;
182.     }
183.  
184.     void Generate()
185.     {
186.         List<Vector2> points = new List<Vector2>();
187.         for (int i = 0; i < VertexCount; i++)
188.         {
189.             points.Add(Random.insideUnitCircle * new Vector2(-5f, 5f));
190.         }
191.         _Triangles = Triangulation(points);
192.     }
193.  
194.     void Start()
195.     {
196.         _Material = new Material(Shader.Find("Hidden/Internal-Colored"));
197.         Generate();
198.     }
199.  
200.     void Update()
201.     {
202.         if (Input.GetKeyDown(KeyCode.Space))
203.         {
204.             Generate();
205.         }
206.     }
207.  
208.     void OnRenderObject()
209.     {
210.         GL.PushMatrix();
211.         _Material.SetPass(0);
212.         GL.Begin(GL.LINES);
213.         for (int i = 0; i < _Triangles.Count; i++)
214.         {
215.             Triangle triangle = _Triangles[i];
216.             GL.Vertex3(triangle.a.x, triangle.a.y, 0);
217.             GL.Vertex3(triangle.b.x, triangle.b.y, 0);
218.             GL.Vertex3(triangle.b.x, triangle.b.y, 0);
219.             GL.Vertex3(triangle.c.x, triangle.c.y, 0);
220.             GL.Vertex3(triangle.c.x, triangle.c.y, 0);
221.             GL.Vertex3(triangle.a.x, triangle.a.y, 0);
222.         }
223.         GL.End();
224.         GL.PopMatrix();
225.     }
226.  
227.     void OnDestroy()
228.     {
229.         Destroy(_Material);
230.     }
231. }

 

Did you write it?

 

Yes. I translated pseudocode from Wikipedia into Unity C# and I added rest of code from myself. Algorithm seems to work properly, but opinions are welcomed.

A few months ago I experimented with techniques like Voronoi diagrams and Delaunay Triangulation and these are examples of the results.

GPU-based 2D Voronoi diagram generation, technique Cone Projection, scene with 50 000 seeds:

Code (CSharp):

1. using UnityEngine;
2. using System.Collections.Generic;
3. using System.Runtime.InteropServices;
4. #if UNITY_EDITOR
5. using UnityEditor;
6. #endif
7.  
8. public class VoronoiCones : MonoBehaviour
9. {
10.     [Header("Cone Settings")]
11.     [SerializeField] Shader _Shader;
12.     [SerializeField] int _Resolution = 4096;
13.     [SerializeField] bool _Animation = true;
14.     [SerializeField] [Range(8,  128)] int _ConeSides = 64;
15.     [SerializeField] [Range(0f, 10f)] float _ConeRadius = 1f;
16.     [SerializeField] [Range(0f, 12f)] float _ConeHeight = 10f;
17.     [Header("Seed Settings")]
18.     [SerializeField] [Range(1f, 1048576f)] int _SeedCount = 2048;
19.     [SerializeField] [Range(0f, 1f)] float _SeedSize = 0.2f;
20.     Material _Material;
21.     ComputeBuffer _Cone, _Seeds;
22.     Matrix4x4 _ModelViewProjection;
23.     RenderTexture _RenderTexture;
24.     GameObject _Plane;
25.  
26.     struct Seed
27.     {
28.         public Vector2 Location;
29.         public Vector3 Color;
30.     };
31.  
32.     List<Vector3> GenerateCone(int sides, float radius, float height)
33.     {
34.         List<Vector3> vertices = new List<Vector3>();
35.         List<Vector2> circle = new List<Vector2>();
36.         float radians = 0.01745329251f;
37.         float step = 360f / (float) sides * radians;
38.         for (int i = 0; i <= sides; i++)
39.         {
40.             float x = radius * Mathf.Cos(i * step);
41.             float y = radius * Mathf.Sin(i * step);
42.             circle.Add(new Vector2(x, y));
43.         }
44.         for (int i = 0; i < sides; i++)
45.         {
46.             vertices.Add(new Vector3(0f, 0f, -height));
47.             vertices.Add(new Vector3(circle[i + 1].x, circle[i + 1].y, 0f));
48.             vertices.Add(new Vector3(circle[i + 0].x, circle[i + 0].y, 0f));
49.         }
50.         return vertices;
51.     }
52.  
53.     void CreateCones()
54.     {
55.         if (_Shader == null) _Shader = Shader.Find("Hidden/VoronoiCones");
56.         _Material = new Material(_Shader);
57.         _RenderTexture = new RenderTexture(_Resolution, _Resolution, 16, RenderTextureFormat.ARGB32);
58.         _RenderTexture.Create();
59.         List<Vector3> vertices = GenerateCone(_ConeSides, _ConeRadius, _ConeHeight);
60.         _Cone = new ComputeBuffer(vertices.Count, Marshal.SizeOf(typeof(Vector3)), ComputeBufferType.Default);
61.         _Cone.SetData(vertices);
62.         _ModelViewProjection.SetRow(0, new Vector4(0.2f,  0.0f,        0.0f, 0.0f)); //orto size = 5, near = -15 and far = 0
63.         _ModelViewProjection.SetRow(1, new Vector4(0.0f, -0.2f,        0.0f, 0.0f));
64.         _ModelViewProjection.SetRow(2, new Vector4(0.0f,  0.0f, -0.0666667f, 0.0f));
65.         _ModelViewProjection.SetRow(3, new Vector4(0.0f,  0.0f,        0.0f, 1.0f));
66.         _Plane.GetComponent<Renderer>().sharedMaterial.mainTexture = _RenderTexture;
67.     }
68.  
69.     void CreateSeeds()
70.     {
71.         Seed[] seeds = new Seed[_SeedCount];
72.         for (int i = 0; i < seeds.Length; i++)
73.         {
74.             float x = UnityEngine.Random.Range(-5f, 5f);
75.             float y = UnityEngine.Random.Range(-5f, 5f);
76.             float r = UnityEngine.Random.Range( 0f, 1f);
77.             float g = UnityEngine.Random.Range( 0f, 1f);
78.             float b = UnityEngine.Random.Range( 0f, 1f);
79.             seeds[i] = new Seed{Location = new Vector2(x, y), Color = new Vector3(r, g, b)};
80.         }
81.         _Seeds = new ComputeBuffer(seeds.Length, Marshal.SizeOf(typeof(Seed)), ComputeBufferType.Default);
82.         _Seeds.SetData(seeds);
83.     }
84.  
85.     void DeleteCones()
86.     {
87.         if (_Cone != null) _Cone.Release();
88.         if (_Material != null) Destroy(_Material);
89.         if (_RenderTexture != null) _RenderTexture.Release();
90.     }
91.  
92.     void DeleteSeeds()
93.     {
94.         if (_Seeds != null) _Seeds.Release();
95.     }
96.  
97.     public void ApplyCones()
98.     {
99.         DeleteCones();
100.         CreateCones();
101.     }
102.  
103.     public void ApplySeeds()
104.     {
105.         DeleteSeeds();
106.         CreateSeeds();
107.     }
108.  
109.     void Start()
110.     {
111.         _Plane = GameObject.CreatePrimitive(PrimitiveType.Plane);
112.         CreateCones();
113.         CreateSeeds();
114.     }
115.  
116.     void OnRenderObject()
117.     {
118.         RenderTexture current = RenderTexture.active;
119.         RenderTexture.active = _RenderTexture;
120.         _Material.SetPass(0);
121.         _Material.SetBuffer("_Cone", _Cone);
122.         _Material.SetBuffer("_Seeds", _Seeds);
123.         _Material.SetMatrix("_ModelViewProjection", _ModelViewProjection);
124.         _Material.SetInt("_Animation", System.Convert.ToInt32(_Animation));
125.         _Material.SetFloat("_SeedSize", _SeedSize);
126.         _Material.SetFloat("_ConeHeight", _ConeHeight);
127.         GL.Clear(true, true, Color.clear);
128.         Graphics.DrawProceduralNow(MeshTopology.Triangles, _Cone.count, _Seeds.count);
129.         RenderTexture.active = current;
130.     }
131.  
132.     void OnDestroy()
133.     {
134.         DeleteCones();
135.         DeleteSeeds();
136.     }
137. }
138.  
139. #if UNITY_EDITOR
140. [CustomEditor(typeof(VoronoiCones))]
141. public class VoronoiConesEditor : Editor
142. {
143.     public override void OnInspectorGUI()
144.     {
145.         DrawDefaultInspector();  
146.         VoronoiCones vc = (VoronoiCones)target;
147.         if(GUILayout.Button("Apply Cone Settings")) vc.ApplyCones();
148.         if(GUILayout.Button("Apply Seed Settings")) vc.ApplySeeds();
149.     }
150. }
151. #endif

Code (CSharp):

1. Shader "Hidden/VoronoiCones"
2. {
3.     Subshader
4.     {
5.         Pass
6.         {
7.             CGPROGRAM
8.             #pragma vertex VSMain
9.             #pragma fragment PSMain
10.             #pragma target 5.0
11.  
12.             struct Seed
13.             {
14.                 float2 Location;
15.                 float3 Color;
16.             };
17.    
18.             StructuredBuffer<float3> _Cone;
19.             StructuredBuffer<Seed> _Seeds;
20.             float4x4 _ModelViewProjection;
21.             int _Animation;
22.             float _SeedSize, _ConeHeight;
23.  
24.             float4 VSMain (uint id : SV_VertexID, uint instance : SV_InstanceID, out float3 color : COLOR, out float4 worldPos : WORLDPOS) : SV_POSITION
25.             {
26.                 worldPos = float4(_Cone[id] + float3(_Seeds[instance].Location, 0.0), 1.0);
27.                 worldPos.xy += _Animation * float2(sin(_Time.g + instance), cos(_Time.g - instance)) * 0.5;
28.                 color = _Seeds[instance].Color;
29.                 return mul(_ModelViewProjection, worldPos);
30.             }
31.  
32.             float4 PSMain (float4 vertex : SV_POSITION, float3 color : COLOR, float4 worldPos : WORLDPOS) : SV_Target
33.             {
34.                 return worldPos.z <= (-_ConeHeight + _SeedSize) ? float4(0, 0, 0, 1) : float4(color, 1);
35.             }
36.             ENDCG
37.         }
38.     }
39. }

GPU-based 3D Voronoi diagram generation, technique Jump Flooding, scene with 50 000 seeds:

Code (CSharp):

1. using UnityEngine;
2. using UnityEngine.Rendering;
3. using System.Runtime.InteropServices;
4.  
5. public class JumpFlooding3D : MonoBehaviour
6. {
7.     [SerializeField] ComputeShader _ComputeShader;
8.     [SerializeField] Shader _PixelShader;
9.     [SerializeField] int _Resolution = 128;
10.     [SerializeField] bool _Animation = true;
11.     [SerializeField] [Range(1, 100000)] int _SeedCount = 2048;
12.     [SerializeField] [Range(-0.5f, 0.5f)] float _SliceXMin = -0.5f, _SliceXMax = 0.5f;
13.     [SerializeField] [Range(-0.5f, 0.5f)] float _SliceYMin = -0.5f, _SliceYMax = 0.5f;
14.     [SerializeField] [Range(-0.5f, 0.5f)] float _SliceZMin = -0.5f, _SliceZMax = 0.5f;
15.     [SerializeField] [Range( 0.0f, 1.0f)] float _Alpha = 1.0f;
16.     ComputeBuffer _Seeds, _Voxels;
17.     Material _Material;
18.     RenderTexture[] _RenderTextures = new RenderTexture[2];
19.     int _CVID, _BVID, _JFID;
20.     bool _Swap = true;
21.  
22.     struct Seed
23.     {
24.         public Vector3 Location;
25.         public Vector3 Color;
26.     };
27.  
28.     void Start()
29.     {
30.         GameObject cube = GameObject.CreatePrimitive(PrimitiveType.Cube);
31.         cube.transform.position = Vector3.zero;
32.         _Material = new Material(_PixelShader);
33.         cube.GetComponent<Renderer>().sharedMaterial = _Material;
34.         RenderTextureDescriptor descriptor = new RenderTextureDescriptor(_Resolution, _Resolution, RenderTextureFormat.ARGBFloat);
35.         descriptor.dimension = TextureDimension.Tex3D;
36.         descriptor.volumeDepth = _Resolution;
37.         for (int i = 0; i < 2; i++)
38.         {
39.             _RenderTextures[i] = new RenderTexture(descriptor);
40.             _RenderTextures[i].enableRandomWrite = true;
41.             _RenderTextures[i].Create();
42.             _RenderTextures[i].filterMode = FilterMode.Point;
43.         }
44.         _Voxels = new ComputeBuffer(_Resolution * _Resolution * _Resolution, sizeof(float) * 3, ComputeBufferType.Default);  
45.         Seed[] seeds = new Seed[_SeedCount];
46.         for (int i = 0; i < seeds.Length; i++)
47.         {
48.             int x = Random.Range(0, _Resolution);
49.             int y = Random.Range(0, _Resolution);
50.             int z = Random.Range(0, _Resolution);
51.             float r = Random.Range(0f, 1f);
52.             float g = Random.Range(0f, 1f);
53.             float b = Random.Range(0f, 1f);
54.             seeds[i] = new Seed{Location = new Vector3(x, y, z), Color = new Vector3(r, g, b)};
55.         }
56.         _Seeds = new ComputeBuffer(seeds.Length, Marshal.SizeOf(typeof(Seed)), ComputeBufferType.Default);
57.         _Seeds.SetData(seeds);
58.         _CVID = _ComputeShader.FindKernel("ClearVoxelsKernel");
59.         _BVID = _ComputeShader.FindKernel("BuildVoxelsKernel");
60.         _JFID = _ComputeShader.FindKernel("JumpFloodKernel");
61.     }
62.  
63.     void Update()
64.     {
65.         _Material.SetVector("_SliceMin", new Vector3(_SliceXMin, _SliceYMin, _SliceZMin));
66.         _Material.SetVector("_SliceMax", new Vector3(_SliceXMax, _SliceYMax, _SliceZMax));
67.         _Material.SetFloat("_Alpha", _Alpha);
68.         _ComputeShader.SetInt("_Resolution", _Resolution);
69.         _ComputeShader.SetInt("_Animation", System.Convert.ToInt32(_Animation));
70.         _ComputeShader.SetFloat("_MaxSteps", Mathf.Log((float)_Resolution, 2.0f));
71.         _ComputeShader.SetFloat("_Time", Time.time);
72.         _ComputeShader.SetBuffer(_CVID, "_Voxels", _Voxels);
73.         _ComputeShader.Dispatch(_CVID, _Resolution / 8, _Resolution / 8, _Resolution / 8);
74.         _ComputeShader.SetBuffer(_BVID, "_Seeds", _Seeds);
75.         _ComputeShader.SetBuffer(_BVID, "_Voxels", _Voxels);
76.         _ComputeShader.Dispatch(_BVID, (_Seeds.count + 8) / 8, 1, 1);
77.         int frameCount = 0;
78.         for (int i = 0; i < _Resolution; i++)
79.         {
80.             _ComputeShader.SetInt("_Frame", frameCount);
81.             int r = System.Convert.ToInt32(!_Swap);
82.             int w = System.Convert.ToInt32(_Swap);
83.             _ComputeShader.SetTexture(_JFID, "_Texture3D", _RenderTextures[r]);
84.             _ComputeShader.SetTexture(_JFID, "_RWTexture3D", _RenderTextures[w]);
85.             _ComputeShader.SetBuffer(_JFID, "_Voxels", _Voxels);
86.             _ComputeShader.Dispatch(_JFID, _Resolution / 8, _Resolution / 8, _Resolution / 8);
87.             _Material.SetTexture("_Volume", _RenderTextures[w]);
88.             _Swap = !_Swap;
89.             frameCount++;
90.         }
91.     }
92.  
93.     void OnDestroy()
94.     {
95.         Destroy(_Material);
96.         _Seeds.Release();
97.         _Voxels.Release();
98.         for (int i = 0; i < 2; i++) _RenderTextures[i].Release();
99.     }
100. }

Code (CSharp):

1. Shader "JumpFlooding3D"
2. {
3.     SubShader
4.     {
5.         Tags { "Queue" = "Transparent" "RenderType" = "Transparent" }
6.         Blend One OneMinusSrcAlpha
7.         Pass
8.         {
9.             CGPROGRAM
10.             #pragma vertex VSMain
11.             #pragma fragment PSMain
12.  
13.             sampler3D _Volume;
14.             float _Alpha;
15.             float3 _SliceMin, _SliceMax;
16.  
17.             float4 FloatToRGBA( float f )
18.             {
19.                 uint q = (uint)(f * 256.0 * 256.0 * 256.0 * 256.0);
20.                 uint r = (uint)(q / (256 * 256 * 256) % 256);
21.                 uint g = (uint)((q / (256 * 256)) % 256);
22.                 uint b = (uint)((q / (256)) % 256);
23.                 uint a = (uint)(q % 256);
24.                 return float4(r / 255.0, g / 255.0, b / 255.0, a / 255.0);
25.             }
26.  
27.             float4 VSMain (float4 vertex : POSITION, out float4 localPos : LOCALPOS, out float3 direction : DIRECTION) : SV_POSITION
28.             {
29.                 localPos = vertex;
30.                 direction = mul(unity_ObjectToWorld, vertex).xyz - _WorldSpaceCameraPos;
31.                 return UnityObjectToClipPos(vertex);
32.             }
33.  
34.             float4 PSMain (float4 vertex : SV_POSITION, float4 localPos : LOCALPOS, float3 direction : DIRECTION) : SV_Target
35.             {
36.                 float3 ro = localPos;
37.                 float3 rd = mul(unity_WorldToObject, float4(normalize(direction), 1));
38.                 float4 result = float4(0, 0, 0, 0);
39.                 int steps = 256;
40.                 float t = 2.0 / float(steps);
41.                 for (int i = 0; i < steps; i++)
42.                 {
43.                     if(max(abs(ro.x), max(abs(ro.y), abs(ro.z))) < 0.500001f)
44.                     {
45.                         float4 voxel = tex3D(_Volume, ro + float3(0.5f, 0.5f, 0.5f));
46.                         float4 color = float4(FloatToRGBA( voxel.w ).rgb, 1.0);
47.                         color.a *= _Alpha;
48.                         bool blend = all(ro > _SliceMin) && all(ro < _SliceMax);
49.                         result.rgb += blend ? (1.0 - result.a) * color.a * color.rgb : 0..xxx;
50.                         result.a += blend ? (1.0 - result.a) * color.a : 0.0;
51.                         ro += rd * t;
52.                     }
53.                 }
54.                 return result;
55.             }
56.             ENDCG
57.         }
58.     }
59. }

Code (CSharp):

1. #pragma kernel ClearVoxelsKernel
2. #pragma kernel BuildVoxelsKernel
3. #pragma kernel JumpFloodKernel
4.  
5. struct Seed
6. {
7.     float3 Location;
8.     float3 Color;
9. };
10.  
11. Texture3D<float4> _Texture3D;
12. RWTexture3D<float4> _RWTexture3D;
13. RWStructuredBuffer<float3> _Voxels;
14. RWStructuredBuffer<Seed> _Seeds;
15. uint _Frame, _Resolution, _Animation;
16. float _MaxSteps, _Time;
17.  
18. float RGBAToFloat( float4 rgba )
19. {
20.     uint r = (uint)(rgba.x * 255);
21.     uint g = (uint)(rgba.y * 255);
22.     uint b = (uint)(rgba.z * 255);
23.     uint a = (uint)(rgba.w * 255);
24.     uint q = (r << 24) + (g << 16) + (b << 8) + a;
25.     return q / (256.0 * 256.0 * 256.0 * 256.0);
26. }
27.  
28. float4 FloatToRGBA( float f )
29. {
30.     uint q = (uint)(f * 256.0 * 256.0 * 256.0 * 256.0);
31.     uint r = (uint)(q / (256 * 256 * 256) % 256);
32.     uint g = (uint)((q / (256 * 256)) % 256);
33.     uint b = (uint)((q / (256)) % 256);
34.     uint a = (uint)(q % 256);
35.     return float4(r / 255.0, g / 255.0, b / 255.0, a / 255.0);
36. }
37.  
38. float4 JFA3D (float3 fragCoord, float level)
39. {
40.     float range = clamp(level - 1.0, 0.0, _MaxSteps);
41.     float stepwidth = floor(exp2(_MaxSteps - range) + 0.5);
42.     float bestDistance = 9999.0;
43.     float3 bestCoord = float3(0.0, 0.0, 0.0);
44.     float3 bestColor = float3(0.0, 0.0, 0.0);
45.     for (int z = -1; z <= 1; ++z)
46.     {
47.         for (int y = -1; y <= 1; ++y)
48.         {
49.             for (int x = -1; x <= 1; ++x)
50.             {
51.                 float3 neighbour = fragCoord + float3(x,y,z) * stepwidth;
52.                 float4 source = _Texture3D.Load(int4(neighbour, 0));
53.                 float3 seedCoord = source.xyz;
54.                 float3 seedColor = FloatToRGBA( source.w ).xyz;
55.                 float magnitude = length(seedCoord - fragCoord);
56.                 if ((seedCoord.x != 0.0 || seedCoord.y != 0.0 || seedCoord.z != 0.0) && magnitude < bestDistance)
57.                 {
58.                     bestDistance = magnitude;
59.                     bestCoord = seedCoord;
60.                     bestColor = seedColor;
61.                 }
62.             }
63.         }
64.     }
65.     return float4(bestCoord, RGBAToFloat(float4(bestColor, 1.0)));
66. }
67.  
68. [numthreads(8,8,8)]
69. void ClearVoxelsKernel (uint3 id : SV_DispatchThreadID)
70. {
71.     uint instance = id.x * _Resolution * _Resolution + id.y * _Resolution + id.z;
72.     _Voxels[instance] = float3(-1.0, -1.0, -1.0);
73. }
74.  
75. [numthreads(8,1,1)]
76. void BuildVoxelsKernel (uint3 id : SV_DispatchThreadID)
77. {
78.     float factor = pow(_Resolution / 128.0, 4.0);
79.     float angle = _Time * 3.0 + id.x;
80.     _Seeds[id.x].Location += _Animation * float3(sin(angle), cos(angle), cos(1.0 - angle)) * factor;
81.     _Seeds[id.x].Location = clamp(_Seeds[id.x].Location, (float3)0.0, (float3)(_Resolution - 1));
82.     int3 location = int3(_Seeds[id.x].Location);
83.     int instance = location.x * _Resolution * _Resolution + location.y * _Resolution + location.z;
84.     _Voxels[instance] = _Seeds[id.x].Color;
85. }
86.  
87. [numthreads(8,8,8)]
88. void JumpFloodKernel (uint3 id : SV_DispatchThreadID)
89. {
90.     float3 fragCoord = float3(id.x, id.y, id.z);
91.     if (_Frame == 0u)
92.     {
93.         uint instance = id.x * _Resolution * _Resolution + id.y * _Resolution + id.z;
94.         float3 buffer = _Voxels[instance];
95.         _RWTexture3D[id] = (buffer.x < 0.0) ? float4(0,0,0,1) : float4(fragCoord, RGBAToFloat(float4(buffer, 1.0)));
96.         return;
97.     }
98.     _RWTexture3D[id] = JFA3D(fragCoord, floor(float(_Frame)));
99. }

Voronoi Dual Graph. Nothing more than Delaunay triangulation of Voronoi Diagram, in real-time.

Code (CSharp):

1. using System.Collections;
2. using System.Collections.Generic;
3. using UnityEngine;
4. using System.Runtime.InteropServices;
5.  
6. public class VoronoiDualGraph : MonoBehaviour
7. {
8.     [SerializeField] ComputeShader _ComputeShader;
9.     [SerializeField] Shader _VertexPixelShader;
10.     [SerializeField] int _SeedCount = 2048;
11.     [SerializeField] int _Resolution = 1024;
12.     [SerializeField] bool _Animation = true;
13.     ComputeBuffer _Seeds, _Triangles, _IndirectBuffer, _CounterBuffer;
14.     Material _Material;
15.     RenderTexture _RenderTexture;
16.     int _VK, _DK;
17.     Seed[] _SeedArray;
18.  
19.     struct Seed
20.     {
21.         public Vector2 Location;
22.         public Vector3 Color;
23.     };
24.  
25.     struct Triangle
26.     {
27.         public Vector2 A;
28.         public Vector2 B;
29.         public Vector2 C;
30.     }
31.  
32.     void Start()
33.     {
34.         _Material = new Material(_VertexPixelShader);
35.         _RenderTexture = new RenderTexture(_Resolution, _Resolution, 0, RenderTextureFormat.ARGBFloat);
36.         _RenderTexture.enableRandomWrite = true;
37.         _RenderTexture.Create();
38.         _RenderTexture.filterMode = FilterMode.Point;
39.         _RenderTexture.wrapMode = TextureWrapMode.Clamp;
40.         _SeedArray = new Seed[_SeedCount];
41.         for (int i = 0; i < _SeedArray.Length; i++)
42.         {
43.             float x = UnityEngine.Random.Range(16f, _Resolution - 16);
44.             float y = UnityEngine.Random.Range(16f, _Resolution - 16);
45.             float r = UnityEngine.Random.Range(0.1f, 0.9f);
46.             float g = UnityEngine.Random.Range(0.1f, 0.9f);
47.             float b = UnityEngine.Random.Range(0.1f, 0.9f);
48.             _SeedArray[i] = new Seed{Location = new Vector2(x, y), Color = new Vector3(r, g, b)};
49.         }
50.         _Seeds = new ComputeBuffer(_SeedArray.Length, Marshal.SizeOf(typeof(Seed)), ComputeBufferType.Default);
51.         _Seeds.SetData(_SeedArray);
52.         _Triangles = new ComputeBuffer(_SeedArray.Length * 16, Marshal.SizeOf(typeof(Triangle)), ComputeBufferType.Append);
53.         _IndirectBuffer = new ComputeBuffer (4, sizeof(int), ComputeBufferType.IndirectArguments);
54.         _IndirectBuffer.SetData(new int[] { 0, 1, 0, 0 });
55.         _CounterBuffer = new ComputeBuffer(1, 4, ComputeBufferType.Counter);
56.         GameObject plane = GameObject.CreatePrimitive(PrimitiveType.Plane);
57.         plane.transform.localScale = new Vector3(_Resolution / 10f, _Resolution / 10f, _Resolution / 10f);
58.         plane.transform.position = new Vector3(_Resolution / 2f, 0f, _Resolution / 2f);
59.         plane.transform.eulerAngles = new Vector3(0, 180f, 0f);
60.         plane.GetComponent<Renderer>().sharedMaterial = new Material(Shader.Find("Legacy Shaders/Diffuse"));
61.         plane.GetComponent<Renderer>().sharedMaterial.mainTexture = _RenderTexture;  
62.         _VK = _ComputeShader.FindKernel("VoronoiKernel");
63.         _DK = _ComputeShader.FindKernel("DelaunayKernel");
64.     }
65.  
66.     void OnRenderObject()
67.     {
68.         if (_Animation)
69.         {
70.             for (int i = 0; i < _SeedArray.Length; i++)
71.             {
72.                 _SeedArray[i].Location += new Vector2(Mathf.Cos(Time.time + i + 2), Mathf.Sin(Time.time + i + 2)) * 0.08f;
73.             }
74.             _Seeds.SetData(_SeedArray);
75.         }
76.         _ComputeShader.SetInt("_SeedsCount", _Seeds.count);
77.         _ComputeShader.SetInt("_Resolution", _Resolution);
78.         _ComputeShader.SetTexture(_VK,"_RWTexture2D", _RenderTexture);
79.         _ComputeShader.SetBuffer(_VK, "_Seeds", _Seeds);
80.         _ComputeShader.Dispatch(_VK, _Resolution / 8, _Resolution / 8, 1);
81.         _Triangles.SetCounterValue(0);
82.         _CounterBuffer.SetCounterValue(0);
83.         _ComputeShader.SetTexture(_DK,"_Texture2D", _RenderTexture);
84.         _ComputeShader.SetBuffer(_DK, "_Seeds", _Seeds);
85.         _ComputeShader.SetBuffer(_DK, "_Triangles", _Triangles);
86.         _ComputeShader.SetBuffer(_DK, "_CounterBuffer", _CounterBuffer);
87.         _ComputeShader.Dispatch(_DK, _Resolution / 8, _Resolution / 8, 1);
88.         int[] args = new int[] { 0, 1, 0, 0 };
89.         _IndirectBuffer.SetData(args);
90.         ComputeBuffer.CopyCount(_CounterBuffer, _IndirectBuffer, 0);
91.         _Material.SetPass(0);
92.         _Material.SetBuffer("_TriangleBuffer", _Triangles);
93.         Graphics.DrawProceduralIndirectNow(MeshTopology.Triangles, _IndirectBuffer);
94.     }
95.  
96.     void OnDestroy()
97.     {
98.         if (_Material != null) Destroy(_Material);
99.         if (_RenderTexture != null) _RenderTexture.Release();
100.         if (_Seeds != null) _Seeds.Release();
101.         if (_Triangles != null) _Triangles.Release();
102.         if (_IndirectBuffer != null) _IndirectBuffer.Release();
103.         if (_CounterBuffer != null) _CounterBuffer.Release();
104.     }
105. }

Code (CSharp):

1. Shader "Hidden/VoronoiDualGraph"
2. {
3.     SubShader
4.     {
5.         Cull Off
6.         Pass
7.         {
8.             CGPROGRAM
9.             #pragma vertex VSMain
10.             #pragma fragment PSMain
11.             #pragma target 5.0
12.  
13.             struct Triangle
14.             {
15.                 float2 Vertices[3];
16.             };
17.  
18.             uniform StructuredBuffer<Triangle> _TriangleBuffer;
19.  
20.             float4 VSMain (float4 vertex : POSITION, uint id : SV_VertexID, out float2 barycentric : BARYCENTRIC) : SV_Position
21.             {
22.                 uint index = id % 3u;
23.                 float3 worldPos = float3(_TriangleBuffer[id / 3u].Vertices[index], 0.02);
24.                 barycentric = float2(fmod(index, 2.0), step(2.0, index));
25.                 return UnityObjectToClipPos(float4(worldPos.xzy, 1.0));
26.             }
27.  
28.             float4 PSMain (float4 vertex : SV_POSITION, float2 barycentric : BARYCENTRIC) : SV_Target
29.             {
30.                 float3 coords = float3(barycentric, 1.0 - barycentric.x - barycentric.y);
31.                 float3 df = fwidth(coords);
32.                 float3 wireframe = smoothstep(df * 0.1, df * 0.1 + df, coords);
33.                 if ((1.0 - min(wireframe.x, min(wireframe.y, wireframe.z))) < 0.01) discard;
34.                 return float4(0.0, 0.0, 0.0, 1.0);
35.             }
36.             ENDCG
37.         }
38.     }
39. }

Code (CSharp):

1. #pragma kernel VoronoiKernel
2. #pragma kernel DelaunayKernel
3.  
4. struct Seed
5. {
6.     float2 Location;
7.     float3 Color;
8. };
9.  
10. struct Triangle
11. {
12.     float2 A;
13.     float2 B;
14.     float2 C;
15. };
16.  
17. Texture2D<float4> _Texture2D;
18. RWTexture2D<float4> _RWTexture2D;
19. StructuredBuffer<Seed> _Seeds;
20. AppendStructuredBuffer<Triangle> _Triangles;
21. RWStructuredBuffer<uint> _CounterBuffer;
22. uint _SeedsCount, _Resolution;
23.  
24. float RGBAToFloat( float4 rgba )
25. {
26.     uint r = (uint)(rgba.x * 255.0);
27.     uint g = (uint)(rgba.y * 255.0);
28.     uint b = (uint)(rgba.z * 255.0);
29.     uint a = (uint)(rgba.w * 255.0);
30.     uint q = (r << 24) + (g << 16) + (b << 8) + a;
31.     return float(q) / 4294967296.0;
32. }
33.  
34. float4 FloatToRGBA( float f )
35. {
36.     uint q = (uint)(f * 4294967296.0);
37.     uint r = (uint)((q / 16777216u) % 256u);
38.     uint g = (uint)((q / 65536u) % 256u);
39.     uint b = (uint)((q / 256u) % 256u);
40.     uint a = (uint)(q % 256u);
41.     return float4(r, g, b, a) / 255.0;
42. }
43.  
44. float Circle (float2 p, float2 c, float r)
45. {
46.     return step(length(p - c) - r, 0.0);
47. }
48.  
49. [numthreads(8,8,1)]
50. void VoronoiKernel (uint3 id : SV_DispatchThreadID)
51. {
52.     float2 fragCoord = float2(id.x, id.y);
53.     float4 result = float4(9999.0, 0.0, 0.0, 0.0);
54.     uint index = 0;
55.     for (uint i = 0; i < _SeedsCount; i++)
56.     {
57.         float3 seed = float3(_Seeds[i].Location, RGBAToFloat(float4(_Seeds[i].Color, 1.0)));
58.         float magnitude = distance(fragCoord.xy, seed.xy);
59.         if (magnitude < result.x)
60.         {
61.             result = float4(magnitude, seed);
62.             index = i;
63.         }
64.     }
65.     float3 circle = Circle(fragCoord, result.yz, 1.0).xxx;
66.     _RWTexture2D[id.xy] = float4(FloatToRGBA(result.w).rgb - circle, float(index));
67. }
68.  
69. [numthreads(8,8,1)]
70. void DelaunayKernel (uint3 id : SV_DispatchThreadID)
71. {
72.     float2 fragCoord = float2(id.x, id.y);
73.     float4 source = _Texture2D.Load(int3(fragCoord, 0));
74.     float4 neighbours[9];
75.     int cells[3] = {int(floor(source.a)), 0, 0};
76.     int count = 1;
77.     int index = 0;
78.     float2 border = float2(0.0, _Resolution - 1u);
79.     for (int y = -1; y <= 1; y++) // get all neighbour pixels
80.     {
81.         for (int x = -1; x <= 1; x++)
82.         {
83.             float2 coords = fragCoord + float2(x, y);
84.             bool off = (coords.x < border.x || coords.x > border.y || coords.y < border.x || coords.y > border.y);
85.             neighbours[index] = off ? source : _Texture2D.Load(int3(coords, 0));
86.             index++;
87.         }
88.     }
89.     for (int i = 1; i < 9; i++) // count distinct pixels in an array
90.     {
91.         int j = 0;
92.         for (j = 0; j < i; j++)
93.         {
94.             if (all(abs(neighbours[i].rgb - neighbours[j].rgb) < 0.001))
95.                 break;
96.         }
97.         if (i == j)
98.         {
99.             cells[count] = int(floor(neighbours[i].a));
100.             count += 1;
101.         }
102.     }
103.     if (count == 3) // if we found a contact point between three Voronoi cells, we can generate new triangle
104.     {
105.         Triangle polygon;
106.         polygon.A = _Seeds[cells[0]].Location;
107.         polygon.B = _Seeds[cells[1]].Location;
108.         polygon.C = _Seeds[cells[2]].Location;
109.         _Triangles.Append(polygon);
110.         _CounterBuffer.IncrementCounter();
111.         _CounterBuffer.IncrementCounter();
112.         _CounterBuffer.IncrementCounter();
113.     }
114. }