ChatGPT:
Yes, that’s correct. Once you have a Delaunay triangulation, generating the corresponding Voronoi diagram is relatively straightforward.
In a Voronoi diagram, each cell consists of all the points that are closer to a particular site (i.e., an input point) than to any other. This makes the Voronoi cell the dual of the Delaunay triangle, because each Delaunay triangle is formed by connecting three sites, and the circumcenter of this triangle (the center of the circumscribed circle that passes through all three vertices of the triangle) is the point that is equidistant from these three sites. Therefore, this circumcenter is the vertex of the three adjacent Voronoi cells corresponding to these sites.
To transform a Delaunay triangulation into a Voronoi diagram, you need to:
- Compute the circumcenter of each Delaunay triangle.
- For each site, find all the circumcenters of the triangles that contain this site as a vertex. These circumcenters are the vertices of the Voronoi cell for this site.
So, you will need to adjust the code of the Delaunay triangulation to also compute and store the circumcenters of the triangles, and then create the Voronoi cells in a second pass.
Nick:
to start, I should compute the circumcenter of each triangle, what does that look like?
ChatGPT:
The circumcenter of a triangle is the center of the circle that passes through all three of its vertices. In the context of a Delaunay triangulation, it’s particularly interesting because the circumcenter of a Delaunay triangle is equidistant from all three vertices of the triangle. This property is used to generate Voronoi diagrams, as mentioned earlier.
Here’s the computation for the circumcenter in C++ code:
struct Circumcenter {
float x, y, r;
};
Circumcenter computeCircumcenter(Vertex* a, Vertex* b, Vertex* c) {
float D = 2.0 * (a->x * (b->y - c->y) + b->x * (c->y - a->y) + c->x * (a->y - b->y));
float Ux = ((a->x * a->x + a->y * a->y) * (b->y - c->y)
+ (b->x * b->x + b->y * b->y) * (c->y - a->y)
+ (c->x * c->x + c->y * c->y) * (a->y - b->y)) / D;
float Uy = ((a->x * a->x + a->y * a->y) * (c->x - b->x)
+ (b->x * b->x + b->y * b->y) * (a->x - c->x)
+ (c->x * c->x + c->y * c->y) * (b->x - a->x)) / D;
float x = Ux;
float y = Uy;
float r = sqrt((x - a->x) * (x - a->x) + (y - a->y) * (y - a->y));
return { x, y, r };
}
You can add this to your Triangle class, where a, b, and c are the vertices of the triangle. The x and y values in the result are the coordinates of the circumcenter, and r is the radius of the circumscribing circle.
Nick:
So, I put the existing code into a class, DelaunayMesh, and we will create a second class, VoronoiMesh; currently we have:
#include <iostream>
#include <random>
#include <algorithm>
#include <vector>
#include <array>
#include <format>
#include <span>
#include <execution>
#include <mutex>
struct Vertex {
float x{ 0 }, y{ 0 };
Vertex() = default;
Vertex(float x, float y)
:x(x), y(y)
{}
};
struct Triangle {
using Edge = std::pair<Vertex*, Vertex*>;
using Edges = std::array<Edge, 3 >;
Vertex* a{ nullptr };
Vertex* b{ nullptr };
Vertex* c{ nullptr };
Triangle(Vertex* a, Vertex* b, Vertex* c) : a(a), b(b), c(c) {
makeClockwise();
}
bool containsVertex(Vertex* v) const {
return a == v || b == v || c == v;
}
bool operator==(const Triangle& other) const {
return containsVertex(other.a) && containsVertex(other.b)
&& containsVertex(other.c);
}
bool circumCircleContains(Vertex* p) const {
//must be clockwise
float ax = a->x - p->x;
float ay = a->y - p->y;
float bx = b->x - p->x;
float by = b->y - p->y;
float cx = c->x - p->x;
float cy = c->y - p->y;
float detAB = ax * ax + ay * ay;
float detBC = bx * bx + by * by;
float detCA = cx * cx + cy * cy;
float det = ax * (by * detCA - detBC * cy)
- ay * (bx * detCA - detBC * cx)
+ detAB * (bx * cy - by * cx);
return det < 0;
}
bool isClockwise(Vertex* a, Vertex* b, Vertex* c) const {
return (b->y - a->y) * (c->x - b->x) > (b->x - a->x) * (c->y - b->y);
}
void makeClockwise() {
if (!isClockwise(a, b, c))
std::swap(b, c);
}
Edges getEdges() const {
return {
//we must ensure that the edge points are in a
//non-reflectable order which would cause sort issues
a < b ? Edge{a, b} : Edge{b, a},
b < c ? Edge{b, c} : Edge{c, b},
c < a ? Edge{c, a} : Edge{a, c}
};
}
struct Circumcenter {
float x{ 0 }, y{ 0 }, r{ 0 };
};
Circumcenter circumcenter;
void computeCircumcenter() {
float D = 2.0 * (a->x * (b->y - c->y) + b->x * (c->y - a->y) + c->x * (a->y - b->y));
float Ux = ((a->x * a->x + a->y * a->y) * (b->y - c->y)
+ (b->x * b->x + b->y * b->y) * (c->y - a->y)
+ (c->x * c->x + c->y * c->y) * (a->y - b->y)) / D;
float Uy = ((a->x * a->x + a->y * a->y) * (c->x - b->x)
+ (b->x * b->x + b->y * b->y) * (a->x - c->x)
+ (c->x * c->x + c->y * c->y) * (b->x - a->x)) / D;
float x = Ux;
float y = Uy;
float r = sqrt((x - a->x) * (x - a->x) + (y - a->y) * (y - a->y));
circumcenter = { x, y, r };
}
};
class DelaunayMesh {
public:
std::vector< Triangle > triangles;
std::vector< Vertex > vertices;
void generate(std::mt19937& random, std::size_t vertexNum = 40, float min=-1.0, float max=1.0 ) {
vertices.resize(vertexNum);
std::uniform_real_distribution<float> reals(min, max);
std::generate(vertices.begin(), vertices.end(), [&]()->Vertex { return { reals(random), reals(random) }; });
std::array<Vertex, 3> superTriangleVertices = { Vertex{0,0}, {0,0}, {0,0} }; //super triangle verticies
auto createSuperTriangle = [&]()->Triangle {
// Find super triangle
float minX = vertices.front().x, minY = vertices.front().y;
float maxX = minX, maxY = minY;
for (const auto& vertex : vertices) {
minX = std::min(minX, vertex.x);
minY = std::min(minY, vertex.y);
maxX = std::max(maxX, vertex.x);
maxY = std::max(maxY, vertex.y);
}
float dx = maxX - minX;
float dy = maxY - minY;
float deltaMax = std::max(dx, dy);
float midX = (minX + maxX) / 2;
float midY = (minY + maxY) / 2;
superTriangleVertices = { Vertex{midX - 20 * deltaMax, midY - deltaMax}
, { midX, midY + 20 * deltaMax}
, { midX + 20 * deltaMax, midY - deltaMax} };
return { &superTriangleVertices[0], &superTriangleVertices[1], &superTriangleVertices[2] };
};
auto superTriangle = createSuperTriangle();
triangles.push_back(superTriangle);
for (auto& v : vertices) {
auto findBadTriangles = [&]() {
std::vector<Triangle*> badTriangles;
std::mutex badMutex;
std::for_each(std::execution::par, triangles.begin(), triangles.end(), [&](auto& t) {
if (t.circumCircleContains(&v)) {
std::scoped_lock lock(badMutex);
badTriangles.push_back(&t);
}
});
return badTriangles;
};
auto createPolygonFromBadTriangles = [&](auto& badTriangles) {
//bad triangles tends to be a small list ( < 10)
std::vector<Triangle::Edge> edges;
for (const auto& triangle : badTriangles) {
auto tEdges = triangle->getEdges();
edges.insert(edges.end(), tEdges.begin(), tEdges.end());
}
std::sort(edges.begin(), edges.end());
std::vector<Triangle::Edge> polygon;
auto it = edges.begin();
do {
auto next = std::find_if_not(it, edges.end(), [&](auto& pair) {return pair == *it; });
auto count = std::distance(it, next);
if (count == 1) polygon.push_back(*it);
it = next;
} while (it != edges.end());
return polygon;
};
auto removeBadTriangles = [&](auto& badTriangles) {
triangles.erase(std::remove_if(std::execution::par, triangles.begin(), triangles.end(),
[&](const Triangle& t) {
auto isBad = std::find(badTriangles.begin(), badTriangles.end(), &t);
return isBad != badTriangles.end();
}), triangles.end());
};
auto createNewTriangles = [&](auto& edges) {
for (const auto& edge : edges) {
triangles.push_back({ &v, edge.first, edge.second });
}
};
auto badTriangles = findBadTriangles();
auto edges = createPolygonFromBadTriangles(badTriangles);
removeBadTriangles(badTriangles);
createNewTriangles(edges);
}
auto removeSuperTriangleConnections = [&]() {
triangles.erase(std::remove_if(triangles.begin(), triangles.end(),
[&](const Triangle& t) {
return t.containsVertex(superTriangle.a)
|| t.containsVertex(superTriangle.b)
|| t.containsVertex(superTriangle.c);
}), triangles.end());
};
removeSuperTriangleConnections();
std::cout << std::format("\nNumber of voronoi triangles: {}\n", triangles.size());
}
void calculateCircumcenters() {
std::for_each(std::execution::par, triangles.begin(), triangles.end(), [&](auto& triangle) { triangle.computeCircumcenter(); });
}
};
class VoronoiMesh {
public:
void generate(DelaunayMesh& delaunay) {
}
};
int main(){
std::mt19937 random;
DelaunayMesh delaunay;
delaunay.generate(random);
delaunay.calculateCircumcenters();
VoronoiMesh voronoi;
voronoi.generate(delaunay);
std::cout << "\nProgram finished\n";
}
This was an excellent starting point for Voronoi and in retrospect it is always a reminder that parallel code is easy to visualize in C++, but difficult to efficiently achieve, despite the deceptive ease at which it can be applied. Sequential execution these days is lightning fast, so all of that ended up being: std::execution::seq and single threaded, because I was going to generate many of these in parallel instead.
BadTriangles turns out to be a poor data structure choice, it would have been better to use a searchable structure such as an unordered_set, but maybe not, badTriangles tends to be a small list. The biggest and easiest optimization that could be made right now is to move the vectors outside of the vertex loop and call clear each iteration.
getEdges, it turns out, can be made much simpler as you will see…
ChatGPT:
Great! You have made good progress in segregating the functionality into classes. Now let’s move to VoronoiMesh generation.
The Voronoi diagram can be generated from a Delaunay triangulation quite simply, because the circumcenter of each triangle in the Delaunay triangulation is a vertex in the Voronoi diagram. For each site (the points you triangulated), the corresponding region in the Voronoi diagram is the polygon whose vertices are the circumcenters of the triangles that have that site as a vertex.
Here’s how to generate the Voronoi diagram:
class VoronoiMesh {
public:
std::vector<std::vector<Vertex>> cells;
std::vector<Vertex> vertices;
void generate(DelaunayMesh& delaunay) {
cells.resize(delaunay.vertices.size());
for (auto& triangle : delaunay.triangles) {
for (auto& vertex : { triangle.a, triangle.b, triangle.c }) {
auto it = std::find(delaunay.vertices.begin(), delaunay.vertices.end(), *vertex);
if (it != delaunay.vertices.end()) {
int index = std::distance(delaunay.vertices.begin(), it);
cells[index].push_back({triangle.circumcenter.x, triangle.circumcenter.y});
}
}
}
for (auto& cell : cells) {
std::sort(cell.begin(), cell.end(), [](Vertex a, Vertex b) {
return std::atan2(a.y, a.x) < std::atan2(b.y, b.x);
});
vertices.insert(vertices.end(), cell.begin(), cell.end());
}
}
};
In this code:
- We are iterating through each triangle in the Delaunay triangulation, and for each vertex in the triangle, we are finding the corresponding Voronoi cell and adding the circumcenter of the triangle to it.
- We are sorting the vertices in each cell in clockwise order to ensure they form a valid polygon. This is done by using the
std::atan2function to compute the angle from the x-axis to the point and sorting based on that. - We then add all the vertices in each cell to the overall vertices list.
The result is a Voronoi diagram where each site in the Delaunay triangulation corresponds to a cell in the Voronoi diagram. Note that this simple implementation does not handle situations where a Voronoi cell should be unbounded. For a complete implementation, additional handling would be required for cells that are connected to the super triangle in the Delaunay triangulation.
This code also doesn’t handle visualization of the Voronoi diagram. Depending on your use case, you might want to generate a graphical representation or output to a specific format.
Nick:
this does not seem correct. The cell shapes are not correct.
ChatGPT:
I apologize for the confusion. Indeed, the previous code was an oversimplification.
To create a Voronoi diagram, you first need to create the Delaunay triangulation, and then for each triangle in the Delaunay triangulation, you can calculate the circumcenter which will serve as the vertex of Voronoi polygons. However, creating the Voronoi cells from these vertices requires more steps:
- Identify all the triangles around each site (vertex from the initial input) in the Delaunay Triangulation.
- Order these triangles in the sequence they appear around the site (clockwise or counter-clockwise).
- Create a polygon by connecting the circumcenters of these triangles in the order determined in step 2. This polygon is the Voronoi cell for that site.
Here is a simple implementation of these steps, please note that this might need to be optimized based on your actual usage:
class VoronoiMesh {
public:
struct Cell {
Vertex site;
std::vector<Vertex> vertices;
};
std::vector<Cell> cells;
void generate(DelaunayMesh& delaunay) {
cells.resize(delaunay.vertices.size());
for (size_t i = 0; i < delaunay.vertices.size(); ++i) {
cells[i].site = delaunay.vertices[i];
std::vector<Triangle*> trianglesAroundSite;
for (auto& triangle : delaunay.triangles) {
if (triangle.containsVertex(&cells[i].site)) {
trianglesAroundSite.push_back(&triangle);
}
}
std::sort(trianglesAroundSite.begin(), trianglesAroundSite.end(),
[&](Triangle* a, Triangle* b) {
double t1 = atan2(a->circumcenter.y - cells[i].site.y, a->circumcenter.x - cells[i].site.x);
double t2 = atan2(b->circumcenter.y - cells[i].site.y, b->circumcenter.x - cells[i].site.x);
return t1 < t2;
});
for (auto& triangle : trianglesAroundSite) {
cells[i].vertices.push_back({triangle->circumcenter.x, triangle->circumcenter.y});
}
}
}
};
Please note that this code assumes that all points are unique and all circumcenters have been computed. It is also a basic implementation and doesn’t handle edge cases such as sites with only 2 triangles around them (these would be along the border of the point set and the corresponding Voronoi cell would be unbounded) or numerical precision issues.
Again, depending on your use case, you might want to generate a graphical representation or output to a specific format.
Nick:
it seems to work except there appear to be some rather very out of bounds verticies pulling the diagram oddly
ChatGPT:
The generation of out-of-bound vertices might be due to the inclusion of circumcenters from the supertriangle in your Voronoi diagram. You should exclude triangles that are connected to the supertriangle when generating the Voronoi diagram because the supertriangle is not a part of your original input, but only an auxiliary construct to simplify the implementation of the Delaunay triangulation.
In other words, the problem arises because we are treating the circumcenters of the supertriangle’s edges as valid vertices in the Voronoi diagram, which we should not.
A solution is to add an isSuperTriangle boolean member to the Triangle structure. This variable should be true for the supertriangle and false for all other triangles. Then, you can add a check to ignore the supertriangle when creating the Voronoi cells:
struct Triangle {
...
bool isSuperTriangle = false;
...
};
void generate(DelaunayMesh& delaunay) {
...
for (auto& triangle : delaunay.triangles) {
if (triangle.containsVertex(&cells[i].site) && !triangle.isSuperTriangle) {
trianglesAroundSite.push_back(&triangle);
}
}
...
}
Remember to set isSuperTriangle = true for the supertriangle during its creation.
This should eliminate the strange distortion of the Voronoi diagram. However, this is only a partial solution. To properly handle the cells that would share an edge with the supertriangle (i.e., cells along the border of your point set), you’ll need to determine a strategy for closing those cells, such as introducing bounding box or handling infinite vertices. Handling these edge cases would be more involved and depend on your particular requirements.
Nick:
but we already removed the super triangle connections earlier?
ChatGPT:
You’re correct, and I apologize for the confusion. The supertriangle should have already been removed from the list of triangles.
When the Voronoi cells are being generated, they may extend beyond the bounds of the original Delaunay triangulation, as they represent regions of the plane closest to the input points. For points on the edge of the point set, their Voronoi regions can theoretically extend to infinity. In practice, these are typically clipped at the bounds of the image or some other predetermined boundary.
That being said, if you’re still seeing oddly-shaped Voronoi cells, the issue could be in the implementation of the Voronoi diagram generation. Here’s a simplified version of what this process might look like:
- Iterate over each triangle in the Delaunay triangulation.
- The circumcenter of each Delaunay triangle becomes a vertex in the Voronoi diagram.
- Each triangle’s circumcenter is connected to the circumcenters of its neighboring triangles in the Delaunay triangulation to form the edges of the Voronoi cells.
If some vertices are very far from the original point set, one possible cause could be incorrect calculation of the circumcenters. You could verify this by visually inspecting the Delaunay triangulation and the corresponding circumcenters. If some circumcenters are far off, it may indicate that their parent triangles in the Delaunay triangulation are large or thin, which can occur if the triangulation is not quite right.
If everything appears correct with the Delaunay triangulation and circumcenters, another area to examine would be the code that connects the Voronoi vertices. There could be an issue with how neighboring triangles are identified or how edges are being drawn between their circumcenters.
Nick:
you are correct, the issue has to do with particularly long triangles about the edges of the delaunay mesh, these generate far off circumcenters.
ChatGPT:
That makes sense. As I mentioned, Voronoi regions for points on the edge of the point set can theoretically extend to infinity. In practice, these cells should be clipped at some reasonable boundary.
To deal with the issue, one approach is to introduce a bounding box or boundary around your Voronoi diagram and clip the Voronoi cells at this boundary.
Here’s a basic outline of the process:
- After you’ve calculated the circumcenters of the Delaunay triangles and created the Voronoi vertices, iterate through each Voronoi cell.
- For each cell, check if any of the cell’s edges extend beyond the bounding box.
- If an edge extends beyond the bounding box, replace that edge with two new edges that intersect the corners of the bounding box. This effectively clips the cell to the bounding box.
Please note that this might result in Voronoi cells that are not perfectly polygonal (i.e., some cells might be irregular polygons), but it will keep all cells within a reasonable boundary. This approach should be sufficient for most graphical or modeling applications.
For more precision, you could look into more complex algorithms that handle unbounded cells in Voronoi diagrams, but these are typically more involved and might be overkill for your purposes.
I rendered the meshes in my Ogre based game renderer using lines and the meshes looked typical. Delaunay and Voronoi geometry are now available to use in Golden Ticket Algorithms. We are going to get into drawing next, but first an intermission of some other content of mine: this is m101 pinwheel galaxy roygbiv 50. When I need a break, sometimes I zoom in on these very large images and mosey around for long adventures…
until next time, when we learn to draw in float space!
Spot C++ 