MATLAB is a really great tool for quickly prototyping an idea and there are lots of nice meshing data structures built in which many people don’t know about. I’m currently working on a problem where I’m trying to fit functions to a surface mesh. My mesh is defined by a set of nodes in the matrix X, where each row gives you (x,y,z) for a node, and the connectivity matrix tri in which each row gives you the the node numbers for each triangular element. In MATLAB I can plot this mesh like so…

`trisurf(tri,X(:,1),X(:,2),X(:,3));`

daspect([1 1 1]);

The daspect command is an important one a lot of people forget about when plotting meshes. It simply corrects the aspect ratio of the plot. As you can see that mesh is very crinkly and I want to smooth it out. My idea is to try

Laplace smoothing which, in its simplest form, works by taking each node and moving it towards the centroid of the nodes connected to it. It might not work with a surface mesh, I’ve read a few articles where they have used it but not tried it myself. I’m writing this blog post literally as I’m coding and I’m going to post the results whatever they turn out to be!

To do this I’m going to write a function

`[ X ] = surfLaplace( tri,X, alpha )`

which takes in the current connectivity and node positions and outputs the new ones. The value alpha will control how far to move the nodes towards the centroid of it’s neighbours, 0 means it won’t move, 0.5 means half way towards the centroid and 1 means all the way.

The first thing I need to do is identify the boundary nodes since we don’t want these to move. This is where MATLAB’s

TriRep class comes in handy. Once we construct it the TriRep representation of our mesh will contain all sorts of handy data structures we can exploit. Construct a TriRep simply with the command

`TR = TriRep(tri, X);`

.

In a triangular mesh the boundary is the set of edges which only belong to a single triangle. MATLAB has a function to find this set using a TriRep called freeBoundary. If we run the code…

`FF = freeBoundary(TR);`

trisurf(FF, X(:,1),X(:,2),X(:,3));

…we get the following plot…

…because the array FF is a list of the boundary edges of our mesh. For the Laplace algorithm we don’t really care about the edge information we just want the node indices. We can get a list of those with some MATLAB aerobics…

`BN = [FF(:,1);FF(:,2)];`

BN = unique(BN);

NN = length(X);

isBound = zeros(1,NN);

isBound(BN) = 1;

%Just to check

hold

trisurf(tri,X(:,1),X(:,2),X(:,3));

scatter3(X(BN,1),X(BN,2),X(BN,3));

daspect([1 1 1]);

Now, if we want to know if index i is on the boundary we just inspect the value of isBound(i). In my case my array of nodes X actually has more than one surface in it, so I’ve also got to generate an array isSurf which tells me if a node is in the surface or not.

With this information we can start traversing the mesh and smoothing. For each node we will need to know the neighboring nodes. I can’t find a function that does this directly so I use the function

`SI = vertexAttachments(TR)`

. Here SI is a cell array with an entry for each node, which are lists of all the triangular elements which are connected to those nodes. With that list we can then collect all the nodes which are in those elements and find the unique values in that list like so for node ip…

`attachElem = SI{ip};`

numElem = length(attachElem);

neigh = tri(attachElem(1),:);

for ie = 2:numElem

neigh = [neigh tri(attachElem(ie),:)];

end

neigh = unique(neigh);

To prove it works here is a plot of the nodes neigh for a particular ip.

We need to remember that the node ip is itself in the list neigh. Now we have the list it is just a simple case of looping round each node and moving it towards the centroid of its neighbours. There might be a more efficient way of writing this but since this is a quick prototype…

`numNeigh = length(neigh);`

xnode = -1*X(ip,1);

ynode = -1*X(ip,2);

znode = -1*X(ip,3);

for ineigh = 1:numNeigh

xnode = xnode + X(neigh(ineigh),1);

ynode = ynode + X(neigh(ineigh),2);

znode = znode + X(neigh(ineigh),3);

end

Xnew(ip,1) = X(ip,1) + alpha*((xnode/numNeigh)-X(ip,1));

Xnew(ip,2) = X(ip,2) + alpha*((ynode/numNeigh)-X(ip,2));

Xnew(ip,3) = X(ip,3) + alpha*((znode/numNeigh)-X(ip,3));

Running the smoothing function with alpha = 0.1 ten times on the surface mesh seems to work pretty well…

The mesh on the left is where we started and the one on the right is where we finished. I’m not entirely sure it is going to work well for all cases, but it might be a good first step. Check the MATLAB documentation for more useful mesh related data structures and any comments post here or talk to me on twitter

@DrSeanWalton
Slight caveat to this post is that I’m using MATLAB version 7.8.0 which uses the

TriRep class which is going to be removed and replaced with the

triangulation class in a future version. It seems to have all the same functionality but with a different name.