Description
A short example of using TetGen to find the Voronoi complex of a point set. Of course, one can apply this any time a Delaunay triangulation is computed, but for simplicity of the examples, we will look at only a point set. One should realize that actually calling TetGen is quite simple; it's the data extraction from the output structure that can be a bit tricky depending on the target application.
To generate voronoi output, use the "v" flag. This will fill in these elements in the output tetgenio:
- numberofvpoints, numberofvedges, numberofvfacets, numberofvcells
- vpointlist, vedgelist, vfacetlist, vcelllist
Potential Gotchas (things I did wrong)
The dual voronoi element should be located at the same index in its list as the primal Delaunay element. However, if redundant points are jettison'd, then you cannot reliably count on the input vertex indices matching with any output you get. Therefore, only use out.pointlist for these purposes. Remember to update the total point count to out.numberofpoints, too.
If you want to compute a Hodge star on edges and dual facets, then you will want to generate edge output from TetGen. This requires TWO 'e' switches: "ee" to get all the edges.
Do not be tempted to create a tetgenbehavior structure and fill in the fields you want; use the character string to specify the behavior. This leads to problems since there are internal switches that need to be set depending on the combination of switches you specify.
Example code
/* This example assumes that PointSet has members NumPoints() which will * return the number of points in the set, and GetPoint(i) which will return * the i-th point. * * Other data structures used are Pt3 (representing a point in R^3), and * Vec3 (representing a displacement in R^3). Subtraction of two points * results in a vector. Vec3Dot computes a dot product between two vectors. */ void PointSetVoronoi(const PointSet &PS){ tetgenio in, out; in.numberofpoints = PS.NumPoints(); in.pointlist = new REAL[3*N]; for(unsigned int i = 0; i < in.numberofpoints; i++){ Pt3 p = PS.GetPoint(i); in.pointlist[3*i+0] = p.e[0]; in.pointlist[3*i+1] = p.e[1]; in.pointlist[3*i+2] = p.e[2]; } // Define the what you want TetGen to do // It is much safer to use the character string list than to explicitly set the // members of tetgenbehavior since there are several dependent members that need // to be set sometimes, as is the case with Voronoi output: //tetgenbehavior behav; //behav.voroout = 1; //behav.useshelles = 1; //tetrahedralize(&behav, &in, &out); // We will not use this approach since it breaks the interface conventions and // does not gain us much in speed if we only need to do the computation once. // We will specify: // Q - quiet, no console output // v - output Voronoi stuff // If you want all edges of the primal complex, then also specify "ee" // to have all edges output. tetrahedralize("Qv", &in, &out); // At this point, out contains all the information about the voronoi complex // We will demonstrate several things you can do with it // Get all voronoi points first, for convenience std::vector<Pt3> vpoints; vpoints.reserve(out.numberofvpoints); for(unsigned int i = 0; i < out.numberofvpoints; i++){ vpoints.push_back(Pt3( out.vpointlist[3*i+0], out.vpointlist[3*i+1], out.vpointlist[3*i+2])); } // Example 1: Iterate over all edges for(unsigned int i = 0; i < out.numberofvedges; i++){ tetgenio::voroedge *cur_edge = &(out.vedgelist[i]); Vec3 normal(cur_edge->vnormal[0], cur_edge->vnormal[1], cur_edge->vnormal[2]); // We can easily recover the primal Delaunay triangle: // It should be at the same index as the Voronoi edge in the list Pt3 primal_triangle[3]; primal_triangle[0] = out.trifacelist[3*i+0]; primal_triangle[1] = out.trifacelist[3*i+1]; primal_triangle[2] = out.trifacelist[3*i+2]; // The dual triangle's plane should be normal to the edge direction. assert(Vec3Dot(normal, primal_triangle[1] - primal_triangle[0]) < EPSILON); assert(Vec3Dot(normal, primal_triangle[2] - primal_triangle[0]) < EPSILON); // There are now 4 cases: // The edge is finite in length, // The edge is a line (shouldn't be possible with more than 3 non-degenerate points) // The edge is a ray emanating from v1 // The edge is a ray emanating from v2 if(cur_edge->v1 == -1){ if(cur_edge->v2 != -1){ // Edge is a ray emanating from v2 Pt3 end2(vpoints[cur_edge->v2]); // Do something with end2 and normal }else{ // Edge is a line } }else if(cur_edge->v2 == -1){ if(cur_edge->v1 != -1){ // Edge is a ray emanating from v1 Pt3 end1(vpoints[cur_edge->v1]); // Do something with end2 and normal }else{ // Edge is a line } }else{ // Edge is finite length Pt3 end1(vpoints[cur_edge->v1]); Pt3 end2(vpoints[cur_edge->v2]); // Do something with the endpoints } } // Example 2: Recover each facet // We will run through all facets and extract the polygon. for(unsigned int i = 0; i < out.numberofvfacets; i++){ tetgenio::vorofacet *cur_facet = &(out.vfacetlist[i]); // First element in the facet list is the number of elements remaining, // which are actual indices into out.vedgelist unsigned int size = cur_facet->elist[0]; // The primal edge should be at the same index in the edge list as this facet // if we used the "ee" flag when calling TetGen. Pt3 primal_edge[2]; primal_edge[0] = out.edgelist[2*i+0]; primal_edge[1] = out.edgelist[2*i+1]; // The adjacent cells are given by cur_facet->c1 and cur_facet->c2, // which are indices into out.vcelllist[] std::vector<unsigned int> vertex_indices; // the final list of indices in circulation order std::set<unsigned int> vertex_seen; // used to prevent duplicating voronoi vertex indices bool infinite_facet = false; // record whether the facet is finite-area or not for(unsigned int j = 0; j < size; j++){ tetgenio::voroedge *cur_edge = &(out.vedgelist[cur_facet->elist[j+1]]); // Technically, the vedgelist is given in circulation order such that if there // are vedge rays, they will be the first and last in the list. However, it is // not necessary to depend on this. if(cur_edge->v1 == -1 || cur_edge->v2 == -1){ infinite_facet = true; if(cur_edge->v1 != -1){ if(vertex_seen.find(cur_edge->v1) == vertex_seen.end()){ vertex_indices.push_back(cur_edge->v1); vertex_seen.insert(cur_edge->v1); } }else if(cur_edge->v2 != -1){ if(vertex_seen.find(cur_edge->v2) == vertex_seen.end()){ vertex_indices.push_back(cur_edge->v2); vertex_seen.insert(cur_edge->v2); } } }else{ if(vertex_seen.find(cur_edge->v1) == vertex_seen.end()){ vertex_indices.push_back(cur_edge->v1); vertex_seen.insert(cur_edge->v1); } if(vertex_seen.find(cur_edge->v2) == vertex_seen.end()){ vertex_indices.push_back(cur_edge->v2); vertex_seen.insert(cur_edge->v2); } } } if(!infinite_facet){ // For example, render the finite-area facet glBegin(GL_POLYGON); for(unsigned int j = 0; j < vertex_indices.size(); j++){ glVertex3fv(vpoints[vertex_indices[j]]); } glEnd(); }else{ // vertex_indices contains the indices into vpoints for the finite // points of the facet. } } // Example 3: Recover each cell // Here we will recover the entire polyhedron defining a Voronoi cell // Note, however, that the facets will be unoriented. To recover // oriented facets, one can exploit the fact that Voronoi cells are // convex, and test the normal of the facet against the vector from a // facet vertex to any point within the cell. If the dot product is // negative, the facet indices need reversal. A convenient interior // point to use is any average of all vertices used in the cell; any // convex combination will do. for(unsigned int i = 0; i < out.numberofvcells; i++){ // The primal vertex is at the same index: Pt3 primal_vertex( out.pointlist[3*j+0], out.pointlist[3*j+1], out.pointlist[3*j+2] ); std::set<unsigned int> cell_indices; // the set of unique Voronoi vertex indices used in the cell for(unsigned int j = 0; j < out.vcelllist[i][0]; j++){ tetgenio::vorofacet *cur_facet = &(out.vfacetlist[out.vcelllist[i][j+1]]); // Perform the iteration over facets like in Example 2, except also // add to cell_indices each time a Voronoi vertex is stored in vertex_indices. } } }