demo04 of Im2mesh package

demo04 - What is inside function im2mesh.
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 poly2mesh 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'))

Shape

Let's start demo. Import image Shape.tif.
im = imread('Shape.tif');
if size(im,3) == 3; im = rgb2gray( im ); end
imshow( im,'InitialMagnification','fit' );

Extract boundaries

We use function im2Bounds to extract boundaries.
% image to polygon boundary
boundsRaw = im2Bounds( im );
% plot boundary using plotBounds
plotBounds(boundsRaw);
We can use function totalNumVertex to get total number of vertices in boundsRaw.
totalNumVertex(boundsRaw)
ans = 1320

Find control points

We use function getCtrlPnts to find and label control points.
% label control points
tf_avoid_sharp_corner = false;
boundsCtrlP = getCtrlPnts( boundsRaw, tf_avoid_sharp_corner, size(im) );
plotBounds(boundsCtrlP, true); % show starting and control points
We can use function totalNumCtrlPnt to get total number of control points.
totalNumCtrlPnt(boundsCtrlP)
ans = 22

Smooth boundary

We use function smoothBounds to smooth boundary.
lambda = 0.5;
mu = -0.5;
iters = 100;
threshold_num_turning = 0;
threshold_num_vert_Smo = 0;
boundsSmooth = smoothBounds( boundsCtrlP, lambda, mu, iters, ...
threshold_num_turning, threshold_num_vert_Smo );
plotBounds(boundsSmooth);

Simplify boundary

We use function simplifyBounds to simplify boundary. Other operation shown here is used to clear up redundant vertices.
% simplify polygon boundary
tolerance = 0.5;
threshold_num_vert_Sim = 0;
boundsSimplified = simplifyBounds( boundsSmooth, tolerance, ...
threshold_num_vert_Sim );
boundsSimplified = delZeroAreaPoly( boundsSimplified );
% clear up redundant vertices
% only control points and turning points will remain
boundsClear = getCtrlPnts( boundsSimplified, false );
boundsClear = simplifyBounds( boundsClear, 0 );
plotBounds(boundsClear);
We can use function totalNumVertex to get total number of vertices in boundsClear.
totalNumVertex(boundsClear)
ans = 104

Area percentage

We can use the following code to calculate area perccentage.
% column vector for intensity
intensity = unique( im );
% calculate the area perccentage of grayscale in image
percent_pixel = getPixelPercent( im );
% calculate the area perccentage in simplified polygonal boundary
percent_polyarea = getPolyShapePercent( boundsSimplified );
% create table
T = table( intensity, percent_pixel, percent_polyarea );
% show table
T
T = 4×3 table
 intensitypercent_pixelpercent_polyarea
1046.657646.8807
210017.989117.7978
318017.255417.2543
425518.097818.0672
Function getPixelPercent is used to calculate the area perccentage of grayscale in an image.
Function getPolyShapePercent is used to calculate the area perccentage in simplified polygonal boundary. Note that function getPolyShapePercent does not work in GNU Octave.
We saw minor changes in area perccentage. That's good.

Generate mesh

We can use function getPolyNodeEdge to get the nodes and edges of polygonal boundary.
Then, we can use function poly2mesh to generate mesh. Function poly2mesh is using MESH2D as mesh generator.
% get nodes and edges of polygonal boundary
[ poly_node, poly_edge ] = getPolyNodeEdge( boundsClear );
grad_limit = 0.25;
hmax = 500;
mesh_kind = 'delaunay';
% generate triangular mesh
[ vert,tria,tnum ] = poly2mesh( poly_node, poly_edge, ...
hmax, mesh_kind, grad_limit );
Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 5 81 96 10 95 205 10 95 205 Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 9 149 96 10 149 234 16 161 716 Smooth triangulation... ------------------------------------------------------- |ITER.| |MOVE(X)| |DTRI(X)| ------------------------------------------------------- 4 168 712 8 14 712
% plot mesh using function plotMeshes
plotMeshes(vert,tria,tnum)
Zoom in
plotMeshes(vert,tria,tnum)
xlim([32.5 62.4])
ylim([34.8 60.8])

Refine mesh near a specific polygon

We can refine the mesh near a specific polygon.
Suppose we are interested in polygon - boundsClear{1}{3}.
poly = boundsClear{1}{3};
plot( poly(:,1), poly(:,2), 'ko-' ); axis equal
We can add more vertices to the interested polygon by repeatedly inserting midpoints to edges. We do this by function insertMidPnt.
polyNew = insertMidPnt( poly );
polyNew = insertMidPnt( polyNew );
polyNew = insertMidPnt( polyNew );
plot( polyNew(:,1), polyNew(:,2), 'ko-' ); axis equal
We use function addPnt2Bound to insert polyNew into global boundaries. Please refer to demo16 about the usage of function addPnt2Bound.
tol_dist = 1e-2; % distance tolerance
boundsNew = addPnt2Bound( polyNew, boundsClear, tol_dist );
% show all vertices
plotBounds( boundsNew, false, 'ko-' );

Generate mesh

% get nodes and edges of polygonal boundary
[ poly_node, poly_edge ] = getPolyNodeEdge( boundsNew );
grad_limit = 0.25;
hmax = 500;
mesh_kind = 'delaunay';
% generate triangular mesh
[ vert,tria,tnum ] = poly2mesh( poly_node, poly_edge, ...
hmax, mesh_kind, grad_limit );
Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 5 102 138 10 118 281 14 121 323 Refine triangulation... ------------------------------------------------------- |ITER.| |CDT1(X)| |CDT2(X)| ------------------------------------------------------- 9 179 138 10 179 294 20 191 1086 20 191 1086 Smooth triangulation... ------------------------------------------------------- |ITER.| |MOVE(X)| |DTRI(X)| ------------------------------------------------------- 4 265 1076 8 49 1076 12 0 1076
plotMeshes(vert,tria,tnum);
Zoom in
plotMeshes(vert,tria,tnum)
xlim([32.5 62.4])
ylim([34.8 60.8])
Awesome!
% reset image size
set(groot, 'DefaultFigurePosition', 'factory')
% end of demo