demo01 of Im2mesh package

demo01 - Demonstrate function im2mesh, which use MESH2D as mesh generator.
Cite as: Ma, J., & Li, Y. (2025). Im2mesh: A MATLAB/Octave package for generating finite element mesh based on 2D multi-phase image (2.1.5). Zenodo. https://doi.org/10.5281/zenodo.14847059
Jiexian Ma, mjx0799@gmail.com. Project website.
Table of Contents

Note

I suggest familiarizing yourself with Im2mesh_GUI before learning Im2mesh package. With graphical user interface, Im2mesh_GUI will help you better understand the workflow and parameters of Im2mesh package.
Im2mesh_GUI: https://www.mathworks.com/matlabcentral/fileexchange/179684-im2mesh_gui-2d-image-to-finite-element-mesh
If you are using Im2mesh package in MATLAB, you need to install MATLAB Image Processing Toolbox and Mapping Toolbox because Im2mesh package use a few functions in these toolboxes.

Setup

Before we start, please set folder "Im2mesh_Matlab" as your current folder of MATLAB.
clearvars
Set default image size.
x = 250; y = 250; width = 250; height = 250;
set(groot, 'DefaultFigurePosition', [x,y,width,height])
% To reset:
% set(groot, 'DefaultFigurePosition', 'factory')
Function im2mesh use a mesh generator called MESH2D (developed by Darren Engwirda). We can use the following command to add the folder 'mesh2d-master' to the path of MATLAB.
addpath(genpath('mesh2d-master'))

im2mesh

Let's start demo01. We'll demostate the usage of function im2mesh.
Import image kumamon.tif.
im = imread('kumamon.tif');
if size(im,3) == 3; im = rgb2gray( im ); end
imshow( im,'InitialMagnification','fit' );
Show the grayscale levels in the image
intensity = unique( im );
intensity
intensity = 4×1 uint8 column vector
0 100 180 255
Let's use function im2mesh to generate meshes based on the kumamon image. We will use default settings. Note that function im2mesh has incorporated the workflow that you saw in im2mesh_GUI: extract polygonal boundaries from image, search & label control points, smooth boundary, simplify boundary, select phase, and generate mesh.
[ vert, tria, tnum ] = im2mesh( im );
Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 4 147 227 10 163 504 12 163 508 Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 6 248 227 10 264 1400 18 270 1799 Smooth triangulation... ------------------------------------------------------- |ITER.| |MOVE(X)| |DTRI(X)| ------------------------------------------------------- 10 19 1783
Function im2mesh has 3 ouptut variables: vert, tria, tnum. They defines linear elements of generated triangular mesh.
vert - Node data. N-by-2 array.
vert(i,1:2) = [x_coordinate, y_coordinate] of the i-th node
tria - Node numbering for each triangle. M-by-3 array.
tria(j,1:3) = [node_numbering_of_3_nodes] of the j-th element
tnum - Label of material phase. P-by-1 array.
tnum(j,1) = k; means the j-th element is belong to the k-th phase
We can plot the mesh using function plotMeshes.
plotMeshes( vert, tria, tnum )

Mesh quality

We can use function tricost to check the mesh quality. Function tricost is written by by Darren Engwirda.
tricost( vert, tria, tnum );
We can get the mean value of Q.
mean(triscr2( vert, tria ))
ans = 0.9550
Wonderful! Let's check the total number of triangles in the mesh.
size( tria, 1 )
ans = 1783
That's a lot of triangles.

Change parameter

We can change the parameters of function im2mesh to see whether it help to reduce the number of triangles.
opt = []; % initialize. opt is a structure array.
opt.tolerance = 1; % default value for opt.tolerance is 0.3
[ vert, tria, tnum ] = im2mesh( im, opt );
Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 4 96 125 8 107 211 Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 4 153 125 10 178 690 15 180 714 Smooth triangulation... ------------------------------------------------------- |ITER.| |MOVE(X)| |DTRI(X)| -------------------------------------------------------
plotMeshes( vert, tria, tnum )
size( tria, 1 )
ans = 708
Now, the number of triangles decreases a lot. That's great.

Parameters of function im2mesh

As I showed before, there're two ways to call function im2mesh.
[ vert, tria, tnum ] = im2mesh( im ); % this use default opt setting
[ vert, tria, tnum ] = im2mesh( im, opt );
Variable opt means options. It's a structure array used to stored input parameters for function im2mesh. We can take a look at the default value of opt.
% opt is a structure array for parameters
opt.tf_avoid_sharp_corner = false;
opt.lambda = 0.5;
opt.mu = -0.5;
opt.iters = 100;
opt.threshold_num_turning = 0;
opt.threshold_num_vert_Smo = 0;
opt.tolerance = 0.3;
opt.threshold_num_vert_Sim = 0;
opt.select_phase = [];
opt.grad_limit = 0.25;
opt.hmax = 500;
opt.mesh_kind = 'delaunay';
opt.tf_mesh = true;
You don't need to memorize the names of these parameters. The meaning of these parameters are listed the file "Im2mesh functions and parameters.pdf". Please refer to "Im2mesh_GUI Tutorial.pdf" about how to choose the value of these parameters. You may change the default values by editting function setOption in "im2mesh.m".

Quadratic element

Function im2mesh can also generate quadratic element via the following syntax.
[ vert, tria, tnum, vert2, tria2 ] = im2mesh( im ); % this use default opt setting
[ vert, tria, tnum, vert2, tria2 ] = im2mesh( im, opt );
Ouput variables - vert, tria, tnum define linear elements.
Ouput variables - vert2, tria2, tnum define quadratic elements.
Their meanings are as follows.
vert - Node data. N-by-2 array.
vert(i,1:2) = [x_coordinate, y_coordinate] of the i-th node
tria - Node numbering for each triangle. M-by-3 array.
tria(j,1:3) = [node_numbering_of_3_nodes] of the j-th element
tnum - Label of material phase. P-by-1 array.
tnum(j,1) = k; means the j-th element is belong to the k-th phase
vert2 - Node data (2nd order element). P-by-2 array. Due to new vertices, the length of vert2 is much longer than vert.
vert2(i,1:2) = [x_coordinate, y_coordinate] of the i-th node
tria2 - Node numbering for each triangle (2nd order element). M-by-6 array.
tria2(j,1:6) = [node_numbering_of_6_nodes] of the j-th element
You can also plot mesh using function plotMeshes, or eveluate mesh quality using function tricost.
plotMeshes( vert2, tria2, tnum );
tricost( vert2, tria2, tnum );

Example 1

Let's show some examples with different opt.
% example 1
opt = []; % reset opt
opt.tolerance = 1e-10; % no polyline simplification
[ vert, tria, tnum ] = im2mesh( im, opt );
Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 6 714 1331 10 739 1681 20 766 3673 20 766 3673 Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 10 883 1331 16 900 1331 20 903 6012 30 904 10841 30 904 10841 Smooth triangulation... ------------------------------------------------------- |ITER.| |MOVE(X)| |DTRI(X)| ------------------------------------------------------- 10 332 10749
plotMeshes( vert, tria, tnum );

Example 2

% example 2
opt = []; % reset opt
opt.tolerance = 1;
opt.select_phase = [1 2 4]; % select phase for meshing
[ vert, tria, tnum ] = im2mesh( im, opt );
Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 3 84 117 7 91 183 Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 4 134 117 10 154 618 15 156 644 Smooth triangulation... ------------------------------------------------------- |ITER.| |MOVE(X)| |DTRI(X)| ------------------------------------------------------- 10 0 636
plotMeshes( vert, tria, tnum );

Example 3

What if we don't need to generate mesh and we want to check the simplified polygonal boundary?
In this special case, we can set "opt.tf_mesh = false;". It means we do not generate mesh. We will use this method in demo13, 15 and 16.
% example 3
opt = []; % reset opt
opt.tolerance = 1;
opt.select_phase = [1 2 4]; % select phase for meshing
opt.tf_mesh = false; % do not generate mesh
bounds = im2mesh( im, opt );
bounds is the simplified polygonal boundaries. We use function plotBounds to plot it.
plotBounds(bounds);
We can plot the polygons with a face color. We will show more about this in demo15. Please refer to the following link for how to set face color: https://www.mathworks.com/help/matlab/ref/polyshape.plot.html
% convert to a cell array of polyshapes
psCell = bound2polyshape(bounds);
figure
hold on; axis equal;
for i = 1: length(psCell)
plot( psCell{i} );
end
hold off
% reset image size
set(groot, 'DefaultFigurePosition', 'factory')
% end of demo