Delaunay triangulation
This page provides basic information about "delaunay_tri_struct". The "delaunay_tri_struct" contains a subroutine that calculates the Delaunay triangulation from given set of vertices(2D/3D).
Background
Delaunay triangulation (Feb. 5, 2016, 06:08 UTC). In Wikipedia: The Free Encyclopedia. Retrieved from https://en.wikipedia.org/wiki/Delaunay_triangulation
In mathematics and computational geometry, a Delaunay triangulation for a set P of points in a plane is a triangulation DT(P) such that no point in P is inside the circumcircle of any triangle in DT(P). Delaunay triangulations maximize the minimum angle of all the angles of the triangles in the triangulation; they tend to avoid skinny triangles. The triangulation is named after Boris Delaunay for his work on this topic from 1934.
For a set of points on the same line there is no Delaunay triangulation (the notion of triangulation is degenerate for this case). For four or more points on the same circle (e.g., the vertices of a rectangle) the Delaunay triangulation is not unique: each of the two possible triangulations that split the quadrangle into two triangles satisfies the "Delaunay condition", i.e., the requirement that the circumcircles of all triangles have empty interiors.
...
Module structure
The delaunay_tri_module is actually an interface for Fortran programs to call C++ functions.
Interface in Fortran:
<source lang="fortran">
!!! Interface for c++ void function.
INTERFACE
!!! This function calculates 2d Delaunay triangulation.
SUBROUTINE get2dtricpp(cp_vertices,vsize,cp_table) BIND(C,name='get2dtricpp')
USE,INTRINSIC :: iso_c_binding
TYPE(c_ptr) :: cp_vertices
INTEGER(c_int) :: vsize
TYPE(c_ptr) :: cp_table
END SUBROUTINE get2dtricpp
!!! This function calculates 3d Delaunay triangulation.
SUBROUTINE get3dtricpp(cp_vertices,vsize,cp_table) BIND(C,name='get3dtricpp')
USE,INTRINSIC :: iso_c_binding
TYPE(c_ptr) :: cp_vertices
INTEGER(c_int) :: vsize
TYPE(c_ptr) :: cp_table
END SUBROUTINE get3dtricpp
!!! This function releases the allocated memory in c++.
SUBROUTINE clear_table(cp_table) BIND(C,name='clear_table')
USE,INTRINSIC :: iso_c_binding
TYPE(c_ptr) :: cp_table
END SUBROUTINE
END INTERFACE
</source>
The actual struct:
<source lang="fortran">
!!! Struct for delaunay tessellation. TYPE delaunay_tri_struct CONTAINS PROCEDURE :: calctri !! (vertices(1:2,1:nsite),table(1:nsite,1:nsite)) END TYPE delaunay_tri_struct
CONTAINS
</source>
C++ functions:
<source lang="cpp">
extern "C"{
/*Void function that takes in flat 2d-array to form 2d-vertices set, calculates delaunay triangulation,
then returns flat T/F table*/
void get2dtricpp(double** verticesflat, int* vsize, bool** table1d){}
/*Void function that takes in flat 3d-array to form 3d-vertices set, calculates delaunay triangulation, then returns flat T/F table*/
void get3dtricpp(double** verticesflat, int* vsize, bool** table1d){}
/*Release memory allocated for sstable1d.*/
void clear_table(bool** table1d){}
}
</source>
Algorithm
The algorithm being used to calculate Delaunay triangulation is Bowyer–Watson algorithm. The c++ code was based on the code provided by : TERCEL::DIARY Here is the pseudocode for the algorithm:
<source lang="fortran">
!!Subroutine for triangulation. SUBROUTINE delaunay_tri(IN[vertex_dict],OUT[simplex_dict])
Create huge d-simplex that contains all vertices. insert "huge d-simplex" into "simplex_dict" FOR (vertex in "vertex_dict") DO
pick next vertex from "vertex_dict"
FOR (d-simplex in "simplex_dict") DO
pick next d-simplex from "simplex_dict"
CALL calculate_circumsphere(IN[d-simplex],OUT[circle])
IF (vertex inside circle) THEN
delete d-simplex from "simplex_dict"
divide d-simplex into (d+1) new and insert into "temp_simplex_dict"
-->{if already exist in "temp_simplex_dict",mark it}
END IF
END DO
FOR (d-simplex in "temp_simplex_dict") DO
pick next d-simplex from "temp_simplex_dict"
IF [(d-simplex not in "simplex_dict").AND.(d-simplex .NOT.marked)] THEN
insert d-simplex inside "temp_simplex_dict" into "simplex_dict"
END IF
END DO
END DO
delete d-simplex that shares vertex with huge d-simplex from "simplex_dict"
END SUBROUTINE delaunay_tri
</source>
Incomplete triangulation problem
When any of the 3 precomputed super triangles' point lies inside a circumcircle determined by any set of points in the given set, the Bowyer–Watson algorithm will suffer from incomplete triangulation problem.
For simplicity, we illustrate this problem in 2d situation, but be aware that this problem applies to arbitrary dimension for Bowyer–Watson algorithm.
The detailed description of the problem is shown here on the Stack Overflow question:Bowyer-Watson algorithm: how to fill “holes” left by removing triangles with super triangle vertices
No perfect solution exists, while a realistic solution that significantly reduces the chance of having the problem is simply making the super triangle as big as possible.