# A Mathematical Analysis of Image Mesh Generation Using Delaunay Triangulation and Image Processing Techniques

José A Rodrigues, Beatriz Vieira

► To cite this version:

José A Rodrigues, Beatriz Vieira. A Mathematical Analysis of Image Mesh Generation Using Delaunay Triangulation and Image Processing Techniques. 4th ROME International Conference on Challenges in Engineering, Medical, Economics and Education: Research & Solutions (CEMEERS-24b), International Forum of Engineers & Practitioners, Dec 2024, Rome, Italy. hal-04822865

HAL Id: hal-04822865

https://hal.science/hal-04822865v1

Submitted on 6 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Mathematical Analysis of Image Mesh Generation Using Delaunay Triangulation and Image Processing Techniques

José A. Rodrigues and Beatriz Vieira

**Abstract**—In this paper, we explore the mathematical foundations and algorithms behind the generation of triangular meshes from image data, focusing on binary images that contain distinct regions, such as a white zone. The approach involves image processing techniques to binarize the input, followed by Delaunay triangulation to create a mesh fitted to the detected region. Key concepts include grid generation, spatial partitioning, and Delaunay triangulation. We explore the mathematical theories supporting these techniques and their computational realization in Python.

**Keywords**—Delaunay Triangulation, Mesh Generation, Image Processing, Binary Image Segmentation, Grid Discretization, Computational Geometry.

## I. INTRODUCTION

Mesh generation is a critical process in computational fields such as computer graphics, numerical simulations, and computer vision. A mesh subdivides a domain into smaller, simple geometric shapes such as triangles, allowing for efficient computation and visualization. In this paper, we generate triangular meshes and fit them to a specific region within an image using Delaunay triangulation.

We detail the step-by-step implementation, discussing how the mathematical concepts like grid discretization, point filtering, and spatial tessellation are used. The code is implemented in Python [1] with common libraries such as OpenCV, NumPy, SciPy, and Matplotlib.

Manuscript received October 14, 2024

José A. Rodrigues (Author), CIMA and Department of Mathematics of ISEL-Higher Institute of Engineering of Lisbon, jose.rodrigues@isel.pt

Beatriz Vieira (Author), Master in Applied Mathematics for Industry, Instituto Superior de Engenharia de Lisboa, Portugal

## A. Motivation

The triangulation and meshing of regions based on image data are essential in many applications like computer vision, image segmentation, finite element analysis (FEA), and physics simulations. This paper explores how Delaunay triangulation is employed to create a triangular mesh that conforms to a selected region in a binary image.

## II. MATHEMATICAL FOUNDATIONS

### A. Image Binarization

An image is represented as a 2D array of pixels, each acting as a point in a Cartesian plane. Image binarization transforms a grayscale image into a binary image, segmenting it into two regions based on pixel intensities. For any pixel intensity function $I(x, y)$, binarization is defined by a threshold function $T$:

$$I_{\text{binary}}(x, y) = \begin{cases} 255 & \text{if } I(x, y) > T \\ 0 & \text{otherwise} \end{cases}$$

In the code, the threshold $T = 240$ is applied to isolate the white zone of the image, resulting in a binary image.

### B. Structured Grid Generation

Grid generation involves discretizing the image into a lattice of points. Given image dimensions $(W, H)$, a grid of points is formed by spacing them evenly at intervals of 10 pixels, which is chosen for simplicity. The mathematical formulation for the grid is:

$$x_i = i \times \Delta x, \quad y_j = j \times \Delta y$$

where $\Delta x = \Delta y = 10$. This discretization divides the image into regularly spaced points, serving as candidate vertices for triangulation.
Once the grid is generated, we aim to select points that lie within the white zone. The binary image is used as a mask to filter out points outside the white region. Mathematically, for a point $p = (x_p, y_p)$, it is retained if:

$$I_{\text{binary}}(x_p, y_p) = 255$$

This forms a subset of the grid points that fit the shape of the white zone. Using set notation, we can express this filtering as:

$$\text{Filtered Points} = \{p \in \text{Grid Points} \mid I_{\text{binary}}(x_p, y_p) = 255\}$$

### C. Unstructured Grid Generation

The structured mesh has limitations, as it cannot be used to create meshes for all types of geometries. To give us greater flexibility, we developed an un-structured mesh approach. This approach allows us to define meshes for any limited 2D connex geometry and provides the freedom to choose the spacing between points. Additionally, it lets us control the number of points along the boundary and within the interior. This flexibility is especially beneficial, as the boundary typically requires a higher density of points compared to the interior.

To create a mesh for a segmented area in an image, a set of points was selected along the boundary of the area, followed by points within the interior.

Let $B = \{b_1, b_2, \dots, b_n\}$ be an ordered set of points representing the boundary of the white area in the image. The objective is to distribute $n$ points evenly along the boundary, ensuring that it is closed. For each segment along the boundary $[b_i, b_{i+1}]$, let $L_i = |b_{i+1} - b_i|$ represent the length of the segment. The position of a new point $p_i$ within the segment $[b_i, b_{i+1}]$ is given by:

$$\mathbf{p}_i = b_i + t(b_{i+1} - b_i)$$

where $t = \frac{d-\text{current\_length}}{L_i}$, with $d$ being the target spacing between consecutive points and current length representing the accumulated length along the boundary up to the start of this segment.

To uniformly distribute points within the region defined by the boundary, Algorithm 1 is applied. In this method, Poisson disk sampling is employed

to ensure a minimum distance between the interior points and the boundary.

### Algorithm 1 Constrained Poisson Disk Sampling for Interior Points

**Input:** Let $R$ be the region of interest defined by the previously defined boundary, $d_{\min}$ the minimum distance between points in the interior of the region, $d_{\text{boundary}}$ the minimum distance between points on the boundary and those in the interior, and $\text{num\_samples}$ the desired number of samples.

**Output:** A set of uniformly distributed points within $R$ that satisfy the distance constraints.

1. Compute the bounding box of $R$ to limit the sampling area.

2. Initialize an empty set of points, $\text{points} = \emptyset$: **while** size of $\text{points} < \text{num\_samples}$

Generate a candidate point $q = (x, y)$ randomly within the bounding box.

if $q \in R$ and $\min_{p \in \text{points}} \|q - p\| > d_{\min}$

if $\min_{b \in \text{boundary}} |q - b| > d_{\text{boundary}}$

Add $q$ to points

end if

end if

end while

return points

### D. Delaunay Triangulation

Delaunay triangulation [2] is the process of defining triangles such that no point lies inside the circumcircle of any triangle. For a set of points $P = \{\mathbf{p}_1, \mathbf{p}_2, \dots, \mathbf{p}_n\}$ on the segmented area, the Delaunay triangulation ensures that for any triangle $T$, its circumcircle does not enclose any point from $P$.

Given three points $\mathbf{p}_i, \mathbf{p}_j, \mathbf{p}_k$, the circumcircle condition is satisfied if:

$$|\mathbf{p}_i - \mathbf{p}_j| \cdot |\mathbf{p}_j - \mathbf{p}_k| \cdot |\mathbf{p}_k - \mathbf{p}_i|$$

This is the primary criterion that Delaunay triangulation algorithms use to form well-conditioned triangles, avoiding slivers.

### E. Laplacian Smoothing

Given a set of points $B$ that define the boundary of a segmented region in an image, the set
Polygon(B) is defined as follows: Supposing, as above, the finite set of points ordered such that each consecutive pair $(b_i, b_{i+1})$ defines an edge of the polygon. Each edge $e_i$ of the polygon is the line segment connecting points $b_i$ and $b_{i+1}$, that is,

$$e_i = \{(1-t)b_i + tb_{i+1} \mid t \in [0, 1]\}, \\ i = 1, 2, \dots, n,$$

where $b_{n+1} = b_1$ to close the loop. The polygon Polygon(B) is the region enclosed by the edges $e_1, e_2, \dots, e_n$.

Our goal is to define a set $P = \{\mathbf{p}_1, \mathbf{p}_2, \dots, \mathbf{p}_n\}$ such that $P \subseteq \text{Polygon}(B)$ subject to certain restrictions to be defined.

Using the points $P$, we aim to define a Delaunay triangulation, ensuring that all triangles lie within Polygon(B). This will be achieved by imposing the condition that the centroid of each triangle belongs to Polygon(B):

For each triangle $t = \{\mathbf{p}_1, \mathbf{p}_2, \mathbf{p}_3\}$, where $\mathbf{p}_i$ are the previously established points, the centroid is given by,

$$\text{centroid} = \frac{\mathbf{p}_1 + \mathbf{p}_2 + \mathbf{p}_3}{3}.$$

The triangle is accepted if

$$\text{centroid} \in \text{Polygon}(B).$$

With a valid mesh in hand, the quality and appearance are ensured by applying Laplacian Smoothing [3] to each triangle.

The set of points $P$, represents the mesh nodes and the set of triangles is denoted by

$$T = \{t_1, t_2, \dots, t_n\}.$$

Additionally, a subset of these points is designated as boundary points, which should remain fixed and not be smoothed. Given a number of smoothing iterations $k$, the objective of Laplacian smoothing is to update each interior point's position to the average position of its neighboring points over $k$ iterations.

For each iteration, consider each interior point $\mathbf{p}_i$ (i.e., points that are not boundary points) and identify the set of neighboring points $\mathbf{p}_j: j \in N(i)$,

where $N(i)$ denotes the set of indices of neighbors of point $i$. Then, update the position of $\mathbf{p}_i$ to the average of its neighboring points as follows:

$$\mathbf{p}_i^{\text{new}} = \frac{1}{|N(i)|} \sum_{j \in N(i)} \mathbf{p}_j$$

Repeat this process for $k$ iterations. At the end of each iteration, an updated set of points $\mathbf{p}_i^{\text{new}}$ is obtained, while boundary points remain unchanged.

To find the neighboring points of a given point $\mathbf{p}_i$ in the mesh, it is sufficient to examine the triangles that include as $\mathbf{p}_i$ a vertex. Any points that share a triangle with $\mathbf{p}_i$ are considered its neighbors, as they form edges with $\mathbf{p}_i$ within the triangle.

### III. COMPUTATIONAL IMPLEMENTATION

We now discuss how the mathematical concepts are implemented in the Python code. The process can be divided into four key functions:

#### A. Image Reading and Binarization

The function `read_image` reads the image, converts it to grayscale, and applies a threshold to generate a binary image. This step identifies the white region in the image:

```python
def read_image(image_path):
    image = cv2.imread(image_path)
    gray = cv2.cvtColor(image,
                         cv2.COLOR_BGR2GRAY)
   _, binary = cv2/threshold(grey,
                             240, 255, cv2.THRESH_BINARY)
    return image, binary
```

#### B. Structured Grid Generation

The grid is created by spacing points every 10 pixels. This regular grid provides potential points for triangulation.

```python
def generateMesh(imageShape):
    height, width = imageShape
    x = np.arange(0, width, 10)
    y = np.arange(0, height, 10)
    xv, yv = np.meshgrid(x, y)
    points = np.vstack([xv, ravel(),
                         yv, ravel()]).T
    return points
```

The following function `adjustMesh` adjusts the grid points by filtering out those that do not lie
within the white zone, using the binary image as a mask.

```python
def adjustMesh(points, binary_image):
    y, x = np wherever(binary_image == 255)
    mask = np zeros(binary_image.shape,
        dtype = bool)
    mask[y, x] = True
    adjusted_points =
        [point for point in points if
         mask[int(point[1]), int(point[0])]]
    return np array(adjusted_points)
```

### C. Unstructured Grid Generation

The first function presented here, `distribute_boundary_points`, is designed to evenly distribute a specified number of points along the perimeter of a closed boundary path. To achieve this, the function first calculates the total perimeter length, then determines the spacing between points, and finally distributes the points along the boundary. The point distribution is accomplished by iterating over each boundary segment to check if the next spaced point fits within the remaining segment. If it does, the function uses linear interpolation between the start and end points of the segment to place the new point. This process repeats until the desired number of points is placed. Finally, the function checks if the last point in the list is close enough to the first point based on the calculated spacing. If not, it adds the first point again to ensure the boundary is fully closed.

```python
def distribute_boundary_points
    (boundary, n_points):
        boundary = np vstack([boundary,
                                boundary[0]])
        total_length = np sum(np sqrt
            (np sum(np diff
                (boundary, axis = 0) ** 2,
                    axis = 1))))
        spacing = total_length / n_points
        boundary_points = []
        current_length = 0
        i = 0
    while len(boundary_points) < n_points:
        start = boundary[i]
        end = boundary[i + 1]
        segment_length = np linalg.
            norm(end - start)
        while current_length +
            segment_length >= spacing:
```

```python
t = (spacing -
    current_length) / segment_length
new_point = start + t * (end - start)
boundary_points.append(new_point)
current_length = 0
start = new_point
segment_length = np.
    linalg .norm(end - start)
    current_length += segment_length
i += 1
if i >= len(boundary) - 1:
    break
if np linalg .norm(
    boundary_points[0] -
    boundary_points[-1]) >
    spacing / 2:
    boundary_points.append
        (boundary_points[0])
    return np array
        (boundary_points)
```

The function `poisson_disk_sampled`, which takes as arguments the minimum distance between points in the interior of the region, the minimum distance between boundary points and interior points, and the desired number of samples, will apply Algorithm 1 to set the points within the interior region.

```python
def poisson_disk_sampled
    (boundary, boundary_points,
        min_dist, num_sampled,
        min_dist_to_boundary):
        polygon = Polygon(boundary)
        bbox = polygon bounds
        points = []
    while len points) < num_sampled:
        x = np random uniform
            (bbox 0, bbox 2])
        y = np random uniform
            (bbox 1, bbox 3])
        candidate = Point(x, y)
        if polygon contains
            (candidate):
                if len points) == 0
                or np min(cdist(
                    [candidate.
                        coordinates[0], points])) > min_dist:
                    if np min(cdist
                        ([candidate .coordinates[0]],
                            boundary_points))
                        > min_dist_to_boundary:
                            points.append([x, y])
    return np array points)
```
Given the image processed by the function `read_image`, the `findContours` function, from the library `cv2`, identifies the contours of the white region within the image. From these contours, the boundary points are extracted, followed by the interior points. The boundary and interior points are then combined to form the final set of points representing the region of interest.

```python
def put_together(image,
                  n_boundary_points,
                  n_interior_points,
                  min_distance,
                  min_dist_to_boundary):
    contours, _ = cv2.findContours
    (binary_image,
     cv2.RETR_EXTERNAL,
     cv2(CHAINciAPPROXtSIMPLE)
    for contour in contours:
        boundary = contour.
        reshape(-1, 2)
        boundary_points =
            distributeXTuniform_
            boundary_points
            (boundary,
             n_boundary_points)
        interior_points =
            poisson_disk王ampling
            (boundary,
             boundary_points,
             min_distance,
             n_interior_points,
             min_dist_to_boundary)
        all_points = np.vstack
        ([boundary_points,
         interior_points])
```

### D. Delaunay Triangulation

Delaunay triangulation is performed on the filtered points, and the resulting triangles are plotted on the original image:

```python
def createMesh(tri_points):
    return Delaunay(tri_points)
```

### E. Laplacian Smoothing

To ensure all triangles remain within the boundary of the mesh, `validTriangles` is applied.

```python
validTriangles = []
for simplex in tri.simplices:
    vertices =
        all_points[simplex]
```

```python
triangle_center =
    vertices.mean(axis=0)
if polygon contains
    (Point(triangle_center)):
    validTriangles.append
    (simplex)
```

The function `find_neighbors` identifies the neighboring points of each previously obtained point.

```python
def find_neighbors(point_index,
                   triangles):
    neighbors = set()
    for tri in triangles:
        if point_index in tri:
            neighbors.update(tri)
        if point_index in neighbors:
            neighbors.remove(point_index)
    return list(neighbors)
```

Considering the interior points of the mesh and the valid triangles, Laplacian Smoothing is applied to these points.

```python
def laplacian_smothing points,
triangles, boundary_points,
iterations=10):
    for _ in range(iterations):
        new_points = np.copy points)
        for i in range(len points)):
            if i < len(boundary_points):
                continue
            neighbors = find_neighbors
            (i, triangles)
            if neighbors:
                new_points[i] =
                    np.mean points
                    [neighbors],
                    axis=0)
            points[:] = new_points
        return points
```

## IV. RESULTS AND DISCUSSION

The result is a triangular mesh that conforms to the white zone of the binary image. Delaunay triangulation ensures well-shaped triangles that avoid sharp angles, which are advantageous for numerical simulations and graphical modeling.

For structured grids, the grid spacing and the threshold used in binarization both influence the quality of the mesh. Structured grids have limitations, as they cannot be used to create meshes for
all types of geometries. Furthermore, the placement
of points is largely predetermined, which can result
in a denser mesh. This yields finer triangulations
but comes at a higher computational cost.

For unstructured grids, we can manipulate the
distance between the points and the number of
points, providing better control. This allows us to
achieve a compromise between finer triangulations
and computational cost.

Figures 1 and 2 show the results of applying Delaunay triangulation to images with irregular shapes, represented in red. The meshes provide good approximations of the shapes, successfully delineating their boundaries. Delaunay triangulation succeeds in avoiding sharp angles, resulting in well-formed meshes for both images.

<figure><img src="image_2.png" /><figcaption>(a) Segmented image illustrating an irregular white shape.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Mesh generated for the irregular shape using Delaunay triangulation.</figcaption></figure>

Fig. 1: Results obtained from applying Delaunay triangulation to segmented images with irregular shapes.

<figure><img src="image_4.png" /><figcaption>(a) Segmented image of a second distinct irregular white shape.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Resulting mesh for the second shape, created with Delaunay triangulation.</figcaption></figure>

Fig. 2: Additional results showcasing Delaunay triangulation applied to different irregular shapes.

Figure 3 shows an example of applying mesh
generation to a practical scenario. A mesh was
created for a human femur, making it possible to

<figure><img src="image_6.png" /></figure>

Fig. 3: Generated mesh overlay for a segmented image of a femur bone, illustrating boundary and interior triangulation.

apply the finite element method (FEM) to calculate
deformations.

V. CONCLUSION

This paper presents a mathematical and compu-
tational approach to generating triangular meshes
based on image data, focusing on Delaunay tri-
angulation. The method effectively generates well-
shaped triangles and fits them to a region of interest
in an image. This work is applicable in fields such
as computer vision, image-based modeling, and
finite element analysis. There are other methods that
segment and generate a grid for an image, such as
FreeFem [4], where the grid is generated by dif-
ferential equations, which is more computationally
expensive. Other methods exist, but they are not
appropriate for the finite element method. Future
work could focus on adaptive meshing techniques
for more complex geometries.

ACKNOWLEDGMENT

This research was partially sponsored with na-
tional funds through the Fundação Nacional para a
Ciência e Tecnologia, Portugal-FCT, under projects
UIDB/04674/2020 (CIMA).
https://doi.org/10.54499/UIDB/04674/2020

REFERENCES

[1] Python Software Foundation. Python Language Reference, version 3.x. Available at: https://www.python.org/.

[2] B. Delaunay, *Sur la sphère vide. A la mémoire de Georges Voronoï*, Bulletin de l'Académie des Sciences de l'URSS. Classe des sciences mathématiques et na, no. 6, 793-800, 1934

[3] Botsch, Mario, et al. Polygon Mesh Processing. AK Pe-ters/CRC Press, 2010. ISBN: 978-1568814261

[4] FreeFem++ Team. FreeFem++: A software for solving par-
tial differential equations. Available at: https://www.freefem.
org/.