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
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
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 );
% 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 );
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