demo15 of Im2mesh package
demo15 - How to edit polygonal boundaries before meshing
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.
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.
Set default image size.
x = 250; y = 250; width = 350; height = 350;
set(groot, 'DefaultFigurePosition', [x,y,width,height])
% 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) );
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.
for i = 1: length(psCell)
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
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;
corners
0.5000 0.5000
92.5000 0.5000
92.5000 80.5000
0.5000 80.5000
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.
for i = 1: length(psCell)
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');
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
% mid point of the right boundary
Add a block
Let's add a block to the left side of psCell.
vertex = [ -20, 0.5; 0.5, 0.5; 0.5, 80.5; -20, 80.5 ];
psBlock = polyshape(vertex);
% add to the end of psCell
for i = 1: length(psCell)
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
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) );
for i = 1: length(psCell)
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].
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 );
for i = 1: length(psCell)
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].
vertex = [ 0 0; 1 0; 1 1; 0 1 ];
psUnitSq = polyshape(vertex); % unit square
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] );
for i = 1: length(psCell)
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 );
% add to the end of psCell
for i = 1: length(psCell)
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 );
plotBounds( boundsNew, false, 'ko-' );
Zoom in
plotBounds( boundsNew, false, 'ko-' )
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 );
plotBounds( boundsNewSimplified, false, 'ko-' );
Zoom in
plotBounds( boundsNewSimplified, false, 'ko-' );
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.
Function poly2mesh uses mesh generator MESH2D (developed by Darren Engwirda). [ poly_node, poly_edge ] = getPolyNodeEdge( boundsNewSimplified );
[ vert,tria,tnum,vert2,tria2 ] = poly2mesh( poly_node, poly_edge, hmax, mesh_kind, grad_limit );
Refine triangulation...
-------------------------------------------------------
|ITER.| |CDT1(X)| |CDT2(X)|
-------------------------------------------------------
5 158 188
10 191 433
15 194 519
Refine triangulation...
-------------------------------------------------------
|ITER.| |CDT1(X)| |CDT2(X)|
-------------------------------------------------------
9 299 188
10 299 442
20 308 1870
22 308 1886
Smooth triangulation...
-------------------------------------------------------
|ITER.| |MOVE(X)| |DTRI(X)|
-------------------------------------------------------
4 529 1860
8 106 1860
12 4 1860
16 0 1860
plotMeshes(vert,tria,tnum);
Zoom in
plotMeshes(vert,tria,tnum);
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.
set(groot, 'DefaultFigurePosition', 'factory')