demo15 of Im2mesh package
demo15 - How to edit polygonal boundaries before meshing
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 = 350; height = 350;
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'))
Get boundries from image
We have demostrate how to do this in demo04.
im = imread('Shape.tif');
if size(im,3) == 3; im = rgb2gray( im ); end
% image to polygon boundary
boundsRaw = im2Bounds( im );
boundsCtrlP = getCtrlPnts( boundsRaw, false, size(im) );
% 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 );
% simplify polygon boundary
tolerance = 0.3; 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 );
Polyshape
It's tedious to edit boundaries directly. Using polyshape object will make things easier. We have demostrated polyshape in demo14. Please take a look at demo14 before you move to demo15.
Here, we use function bound2polyshape to convert a cell array of polygonal boundaries 'boundsClear' to a cell array of polyshape 'psCell'.
psCell = bound2polyshape(boundsClear);
Plot cell array of polyshape.
figure
hold on; axis equal;
for i = 1: length(psCell)
plot(psCell{i});
end
hold off
Coordinates of 4 corners
Before we edit polyshape, we want to know about the x y coordinates of 4 corners in in variable psCell.
The image size is
size(im)
The height and width of the input image is 80 and 92, respectively.
The x y coordinates of 4 corners in variable psCell will be
corners = [0.5, 0.5] + [ 0, 0;
size(im,2), 0;
size(im,2), size(im,1);
0, size(im,1) ];
corners
The lower left corner is [0.5, 0.5]. The upper right corner is [92.5, 80.5].
I forgot why there is 0.5 shifting in x y coordinates since some code in Im2mesh were written by me many years ago. Perhaps I want to be consistent with MATLAB built-in function bwboundaries.
Let's plot them.
figure
hold on; axis equal;
for i = 1: length(psCell)
plot(psCell{i});
end
plot(corners(:, 1), corners(:, 2), 'bo', 'MarkerFaceColor', 'b');
% Annotate each point with its coordinates
for i = 1:size(corners, 1)
text(corners(i, 1), corners(i, 2), ...
sprintf('(%.2f, %.2f)', corners(i, 1), corners(i, 2)), ...
'VerticalAlignment', 'bottom', ...
'HorizontalAlignment', 'center');
end
hold off
Therefore, based on 4 corners, we can know that the mid point of the upper boundary is [ 46.5, 80.5 ], and the mid point of the right boundary is [ 92.5, 40.5 ]. We will use them later.
% mid point of the upper boundary
xMidUp = 46.5;
yMidUp = 80.5;
% mid point of the right boundary
xMidRight = 92.5;
yMidRight = 40.5;
Add a block
Let's add a block to the left side of psCell.
% define a block
vertex = [ -20, 0.5; 0.5, 0.5; 0.5, 80.5; -20, 80.5 ];
psBlock = polyshape(vertex);
% add to the end of psCell
psCell{end+1} = psBlock;
figure
hold on; axis equal;
for i = 1: length(psCell)
plot(psCell{i});
end
hold off
Merge geometry
We have assumed the block is a new phase in the geometry.
What if the block and the background are in the same phase? We can merge the block into the background.
% psCell{1} is the background. psCell{5} is the block.
% merge the block into the background
psCell{1} = union( psCell{1}, psCell{5} );
psCell{5} = []; % delete the block
length( psCell )
ans = 5
The length of psCell is 5 although the 5-th element is [].
To avoid error in the following steps, we need to manually delete empty elements in the cell array.
% delete empty elements in cell array
psCell = psCell( ~cellfun('isempty', psCell) );
figure
hold on; axis equal;
for i = 1: length(psCell)
plot(psCell{i});
end
hold off
Add a notch
Let's add a oval-shape notch to the point [xMidUp, yMidUp].
Define a oval. The center of the oval is [xMidUp, yMidUp].
t = 0.05:0.03:2*pi;
x1 = xMidUp + 5*cos(t);
y1 = yMidUp + 40*sin(t);
psNotch = polyshape(x1,y1);
plot( psNotch ); axis equal
Add a oval-shape notch (by boolean operation - subtract).
for i = 1: length(psCell)
psCell{i} = subtract( psCell{i}, psNotch );
end
figure
hold on; axis equal;
for i = 1: length(psCell)
plot(psCell{i});
end
hold off
Add a rod
Let's add a rod to the mid point of the right boundary. The mid point is [xMidRight, yMidRight].
Define a rod. The center of the rod is [xMidRight, yMidRight].
% define a rod
vertex = [ 0 0; 1 0; 1 1; 0 1 ];
psUnitSq = polyshape(vertex); % unit square
% set length and width
psRod = scale( psUnitSq, [80, 5]);
% set the centroid to [xMidRight, yMidRight]
[xc,yc] = centroid(psRod);
vec = [ xMidRight - xc, yMidRight - yc];
psRod = translate( psRod, vec );
% rotate 30 degree with respect to centroid
[xc,yc] = centroid(psRod);
psRod = rotate( psRod, 30, [xc,yc] );
% plot together
figure
hold on; axis equal;
for i = 1: length(psCell)
plot(psCell{i});
end
plot( psRod );
hold off
We saw that they are overlapped.
Boolean operation
We use boolean operations of polyshape object to remove overlapped regions.
for i = 1: length(psCell)
psCell{i} = subtract( psCell{i}, psRod );
end
% add to the end of psCell
psCell{end+1} = psRod;
figure
hold on; axis equal;
for i = 1: length(psCell)
plot(psCell{i});
end
hold off
Add intersect points
Polyshape objects are not able to find intersect points automatically. Therefore, we convert the cell array of polyshape to a cell array of polygonal boundaries, and use function addIntersectPnts to add intersect points to boundaries.
Note that if we don't add intersect points, mesh generation will fail.
boundsNew = polyshape2bound(psCell);
tol_intersect = 1e-6; % distance tolerance for intersect
boundsNew = addIntersectPnts( boundsNew, tol_intersect );
% show all vertices
plotBounds( boundsNew, false, 'ko-' );
Zoom in
plotBounds( boundsNew, false, 'ko-' )
xlim([27.9 63.0])
ylim([35.9 57.8])
We saw that there are a lot of vertices on the notch.
Simplify boundaries
We want to reduce the number of vertices on the notch.
boundsNewCtrlP = getCtrlPnts( boundsNew, false );
tolerance = 0.02; % for Douglas-Peucker polyline simplification
boundsNewSimplified = simplifyBounds( boundsNewCtrlP, tolerance, 0 );
% show all vertices
plotBounds( boundsNewSimplified, false, 'ko-' );
Zoom in
plotBounds( boundsNewSimplified, false, 'ko-' );
xlim([27.9 63.0])
ylim([35.9 57.8])
Wonderful! We saw that the number of vertices on the notch reduced a lot.
Note that you need to adjust the value of tolerance to achieve your desired reduction.
Generate mesh
We use function poly2mesh to generate mesh.
[ poly_node, poly_edge ] = getPolyNodeEdge( boundsNewSimplified );
hmax = 500;
mesh_kind = 'delaunay';
grad_limit = 0.25;
[ vert,tria,tnum,vert2,tria2 ] = poly2mesh( poly_node, poly_edge, hmax, mesh_kind, grad_limit );
plotMeshes(vert,tria,tnum);
Zoom in
plotMeshes(vert,tria,tnum);
xlim([27.9 63.0])
ylim([35.9 57.8])
Wonderful!
We can also use MATLAB built-in function generateMesh or open-source software Gmsh to generate mesh. We have demonstrated how to do that in other examples of Im2mesh package.
We saw that the mesh in the rod is very coarse. It's possible to refine mesh near a specific polygon boundary. Please refer to demo04 and demo16.
% reset image size
set(groot, 'DefaultFigurePosition', 'factory')
% end of demo