Let's remind briefly the main points of the Delaunay triangulation method^{[1]} together with their numerical implementation using Octave and Matlab software^{[2]}. Let P = {p_{i}, i = 1, 2, … ,N} be a set of points in two – dimensional Euclidean plane ℜ^{2}. They are called forming points of mesh^{[1]}. Let us define the triangle T as a set of three mesh points
T = {t_{j} ∈ P, j = 1,2,3}. 
Using the Delaunay criterion one can generate triangulation where no four points from the set of forming points P are co – circular:
∀ p_{i} ∈ P ∧ p_{i} ∉ T  p_{i}  u > ρ^{2} 
where u is the center of the T triangle and ρ is its radius. The proposed algorithm consists of the following steps:
A = t_{12}  t_{13} 
A_{n} = n ⋅ A 
The determinant of the square matrix D(T) is built on the basis of the set of triangle's nodes given by the equation
D(T) =  det  ( 

) 
Figure shows the main idea of the Delaunay criterion. a) Two triangles (with nodes abc and acd, respectively) are not Delaunay triangles, b) after exchange of the edge ac to the edge bd two new triangles abd and bcd replace the old ones. They are both of the Delaunay type. Circles represent the Delaunay zones.
next the following determinant is calculated in order to find out whether a mesh point p_{i} is outside or inside the Delaunay zone (see Fig. 13)
D(T)_{i} =  det  ( 

) 
The algorithm ends up with the new triangular mesh Ω_{new}.
If you wish you can have an insight into the program delaunay.m that implements the above – presented algorithm using the Octave and Matlab software (it is a part of non – free Octave and Matlab course). One can use also the appropriate Octave and Matlab library function.
Help. In order to find an orientation of a triangle T one can check the sign of A_{n} (see equation). If it is greater than 0 the triangle orientation is clockwise unless counterclockwise. In the latter case, to ensure the clockwise orientation one can once flip up and down matrix in equation then the triangle orientation turns into the opposite one. Obviously, this flipping results in the change of the sign of the matrix determinant D(T) → D(T).