## Posts Tagged ‘triangle mesh’

### Offset surface of triangle mesh in matlab

Wednesday, March 22nd, 2017

Here’s a little demonstration of how to use gptoolbox and MATLAB to generate an offset surfaces of a triangle mesh. This takes a mesh in V,F and creates a mesh SV,SF of the isosurface at signed distance iso:

% Extract offset at minus 3% of bounind box diagonal length
iso = -0.03;
% Resolution grid → resolution of extracted offset surface
side = 60;
% Amount of smoothing to apply to distance field
sigma = 1.4;
bbd = norm(max(V)-min(V));
% Pad generously for positive iso values
[BC,side,r] = voxel_grid([V;max(V)+iso*1;min(V)-iso*1],side);
D = signed_distance(BC,V,F);
D = reshape(D,side([2 1 3]));
% Smooth signed distance field
D = imfilter(D,fspecial('gaussian',9,sigma),'replicate');
BC3 = reshape(BC,[side([2 1 3]) 3]);
% Use matlab's built-in marching cubes iso-surface mesher (works most of the time)
surf = isosurface(BC3(:,:,:,1),BC3(:,:,:,2),BC3(:,:,:,3),D,iso*bbd);
SV = surf.vertices;
SF = surf.faces;

Here’s a blue bunny with a positive offset surface, an orange “cage”:

Here’s a blue bunny with a negative offset surface. This is useful for hollowing out objects to 3d print:

Because the iso-surface extraction will over tesselate low curvature patches of the output surface, it would make a lot of sense to remesh/decimate this mesh.

(to create these fancy renderings:)

clf;
hold on;
t  = tsurf(F,V,'EdgeColor','none',fsoft,  'FaceVertexCData',repmat(blue,size(V,1),1),'FaceAlpha',1+(iso<0)*(0.35-1),fphong);
ts = tsurf(SF,SV,'EdgeAlpha',0.2+(iso<0)*(0-0.2),fsoft,'FaceVertexCData',repmat(orange,size(SV,1),1),fphong,'FaceAlpha',1+(iso>0)*(0.2-1));
apply_ambient_occlusion(ts);
hold off;
axis equal;
view(-20,20)
camlight;
t.SpecularStrength = 0.04;
l = light('Position',[5 -5 10],'Style','local');
set(gca,'pos',[0 0 1 1])
set(gca,'Visible','off');
set(gcf,'Color','w');
drawnow;

### Boolean operations using generalized winding numbers

Tuesday, February 2nd, 2016

I posted a little tech report about how to use the generalized winding number for constructive-solid geometry (CSG) style boolean operations (union, intersection, difference etc.) on nasty meshes.

Code implementing this has existed for a little while now in the mesh_boolean_winding_number function of gptoolbox.

### How does Galerkin multigrid scale for irregular grids?

Tuesday, August 4th, 2015

There are different flavors of the multigrid method. Let’s consider the geometric subclass and within that the “Galerkin” or projection-style multigrid method. This method is favorable because it only requires discretization at the finest level and prolongation and restriction between levels. This means, given a problem A1 x1 = B1 on the finest mesh, we don’t need to know how to create an analogous problem A2 x2 = B2 on the second finest mesh and so on. Instead, we define A2 = R A1 P where R is the #coarse vertices by #fine vertices restriction operator taking fine mesh values to coarse mesh values and P is the #fine vertices by #coarse vertices prolongation operator. A Galerkin multigrid V-cycle looks something like this:

• relax current guess x1 according to A1 x1 = B1
• compute residual r1 = B1 - A1 x1
• restrict residual to coarse mesh r2 = R r1
• restrict system matrix A2 = R A1 P
• solve (recursively) for update A2 u2 = r2
• prolong update and add to guess u1 = P u2, x1 += u1
• relax current guess x1 according to A1 x1 = B1

Often (for good reasons), we take the restriction operator to be the transpose of the prolongation operator R = P'. When we’re working with nested meshes, it’s natural to use the linear interpolation operator as P.

This V-cycle if A1 x1 = B1 is coming from a quadratic energy minimization problem:

min  ½ x1' A1 x1 - x1' B1
x1

Suppose we freeze our current guess x1 and look for the best update u1, we’d change the problem to

min  ½ (x1 + u1)' A1 (x1 + u1) - (x1+u1)' B1
u1

Rearranging terms according to u1 this is

min  ½ u1' A1 u1 - u1' (B1 - A1 x1) + constants
u1

Now we recognize our residual r1 = B1 - A1 x1. Let’s substitute that in:

min  ½ u1' A1 u1 - u1' r1
u1

If we were to optimize this for u1 we’d get a system of equations:

A1 u1 = r1

Instead, let’s make use of our multires hierarchy. The coarse mesh spans a smaller space of functions that the finer mesh. We can take any function living on the coarse mesh u2 and interpolate it to give us a function on the fine mesh u1 = P u2, thus we could consider restricting the optimization search space of u1 functions to those we can create via P u2:

min  ½ (P u2)' A1 (P u2) - (P u2)' r1
u2

rearranging terms and substituting R = P' we get:

min  ½ u2' R A1 P u2 - u2' R r1
u2

finally optimizing this for u2 we get a system of equations:

R A1 P u2 = R r1

or

A2 u2 = r2

This Galerkin-style multigrid is great because we don’t need to worry about re-discretizing the system matrix (e.g. using FEM on the coarse mesh) or worry about redefining boundary conditions on the coarse mesh. The projection takes care of everything.

But there’s a catch.

For regular grids with carefully chosen decimation patterns (#vertices on a side is a power of 2 plus 1) then the projection matrix will be well behaved. Each row corresponding to a fine mesh vertex will at most 1 vertex if perfectly align atop a coarse mesh vertex and will involve the nearest 2 vertices on the coarse mesh if lying on a coarse mesh edge. Then examining a row of A2 = P' A1 P we can read off the change in the stencil. Let’s say A1 is a one-ring Laplacian on the fine mesh. Hitting it on the left with P' will form the matrix P' A1, a row of which corresponds to a coarse mesh vertex and the columns to the fine mesh vertices it “touches”. Since coarse vertices lie exactly at fine mesh vertices, this will simply connect each coarse mesh to its “half-radius-one-ring”, i.e. the fine mesh vertices lying on coarse mesh edges around it. Finally hitting this matrix with P on the right creates A2 = P' A1 P connecting coarse to coarse vertices. Since each coarse vertex was connected to its “half-radius-one-ring” via the fine mesh, now in P' A P` all coarse mesh vertices are connected to its true one-ring neighbors: just like if we built the Laplacian directly on the coarse mesh. The sparsity pattern is maintained.

For irregular grids this is unfortunately not the case. Each row in P will contain, in general, d+1 entries. When we examine P' A1, we see that we smear the stencil around the coarse mesh to include the “half-radius-one-rings” from d+1 nearest fine-mesh vertices. This effect is repeated when finally creating P' A1 P. We see the sparsity worsening with each level.

Here’s a plot of the number non-zeros per row of the projected system matrix for a multires hierarchy with 7 levels.

The orange lines shows the average number of non-zeros per row for the cotangent Laplacian projected down different sequences of random irregular meshes of the unit square. You can see a rough 2x increase in the number of non-zeros in the beginning and then things taper-off as the meshes get very coarse: meshes range from roughly 2562 to 42. The blue lines show the average number of non-zeros per row for the cotangent Laplacian rediscretized (aka rebuilt from scratch) on each mesh: this is and should be roughly 7 for any triangle mesh (the slight decrease is due to the increase effect of the boundary on the average).

This doesn’t look so bad. Even for the third level (562 vertices) the number of non-zeros per row is 25: still _very sparse).

Things get worse, much worse in 3D.

Already the average number of non-zeros per row of a the 3D tet-mesh Laplacian is not neatly bounded to 7, though it tends to be around 16 for a “reasonable mesh”.

In this plot we see the number of non-zeros in the projected system matrix increasing by around a factor of 3 each level, peaking around 112 on the 173 mesh. Notice on the 93=729 vertex grid the average number of non-zeros is 100. That’s a sparsity ratio of 13% or 729,000 non-zeros: not very sparse for such a small mesh. Especially considering that the rediscretization system matrix would have 12 non-zeros per row, a sparsity ratio of 1.6%, or only 8,700 total non-zeros.

I didn’t have the patience to let this run on a larger initial fine-mesh, but you can imagine that the situation gets drastically worse. Galerkin multigrid fails to improve performance and gets choked by memory consumption on irregular 3D tetrahedral meshes using piecewise linear prolongation operators.

$E = mc^2$

### Edge collapse and mesh decimation in libigl

Thursday, April 16th, 2015

Libigl purposefully does not build itself around complicated mesh data-structures. This is freeing for many reasons: it’s very easy to include our code in an arbitrary project, functions are not artificially limited to manifold meshes, array based data structures are easily serialized and fast, matrix operations on vertex positions are directly exposed, etc.

But our choice is also limiting. In particular, combinatorial changes to the mesh are potentially expensive and difficult to implement. Recently, I took a first stab at implementing an efficient edge-collapse routine for libigl. Indeed I’m seeing good performance: O(1) constant time for a single edge collapse and O(log m) time for cost-priority-queue-based collapses. The data structures are a little “messy” in the sense that when edges or faces disappear their rows are not deleted, but rather just redacted: entries are replaced with NULL values. This is because my approach is array-based and constant resizing would ruin performance. Luckily for edge-collapse, the number of elements only decreases. For edge-split I’ll have to think about amortized costs…

For now, check out the updated tutorial entry for libigl’s new mesh decimation and edge collapse features.

The code is “programmable” in the sense I expose function handles for computing edge costs and merged vertex placements. Though the default functions I currently provide are quite naive, this should support rather advanced “Progressive Meshes” style simplification with creases, sharp features, adaptivity, etc.

### Slice through a tet-mesh in matlab

Monday, February 2nd, 2015

Matlab has only rudimentary support for visualizing data defined over a tetrahedral mesh. I use usually external calls to medit. Today I needed to implement slicing through a tet mesh for a different application and noticed that it’s also useful for visualization. I add slice_tets.m to gptoolbox. The slicing works in an obvious, yet vectorized way. The way a plane passes through a tet (if at all) can be described by the number of vertices on either side of the plane: 0-4, 1-3, 2-2, 3-1, 4-0. The 0-4 and 4-0 cases don’t add anything to the slice and are quickly ignored. The 1-3 and 3-1 cases add a single triangle to the slice and if we categorize the vertices with respect to their signed distance to the plane, then we only need to implement the 1-3 case and reverse the signs to handle the 3-1 case. Finally the 2-2 case adds a quad (i.e. two triangles) to the slice, and we should be careful to flip the orientation depending on which vertices end up on either side.

My vectorized matlab identifies the case each tet belongs to and computes triangles or split quads for all of that case in unison.

The slicing is actually pretty fast. For a 600,000 tetrahedra mesh, slicing through the middle takes ~0.11 secs (~9fps). With a little more time spent, I may also keep track of the linear interpolation coordinates so that I can visualize a function defined over the tet-mesh onto the slicing plane.

Here I’m showing a slice through a Poisson equation inside this reclining knight:

### Robust mesh boolean operations in libigl, gptoolbox

Tuesday, November 4th, 2014

I’ve added robust mesh boolean operations to libigl and a mex wrappers for matlab in gptoolbox. For comparison and as an alternative, I also included new wrappers cork’s boolean operations.

Check out the boolean entry in the libigl tutorial.

### MATLAB2014b features anti-aliasing

Monday, October 13th, 2014

Finally. I’m pretty happy about the results:

AO = ambient_occlusion(V,F,V,per_vertex_normals(V,F),1000);
t = tsurf(F,V,fphong,'EdgeColor','none');
C = squeeze(ind2rgb(floor(matrixnormalize(t.FaceVertexCData)*size(colormap,1))+1,colormap));
t.FaceVertexCData = bsxfun(@times,C,1-AO);
t.SpecularStrength = 0.1;
t.DiffuseStrength = 0.1;
t.AmbientStrength = 0.8;
l = light('Position',[1 1 100],'Style','infinite');
l2 = light('Position',[1 -100 1],'Style','infinite');
set(gca,'XTickLabel',[],'YTickLabel',[],'ZTickLabel',[],'Color',[0.94 0.94 0.94]);
set(gcf,'Color','w');

And to spin the camera around:

axis equal
axis vis3d;
camproj('persp');
for f = 1:numel(T)
t = T(f);
view(-cos(t*2*pi)*45,sin(t*2*pi)*45+45);
drawnow;
frame = getframe(gcf);
[SIf,cm] = rgb2ind(frame.cdata,256);
if f == 1
imwrite(SIf,cm,filename,'Loop',Inf,'Delay',0);
else
imwrite(SIf,cm, filename,'WriteMode','append','Delay',0);
end
end

With the awesome but now obsolete myaa.m hacked anti-aliasing, creating this gif would have taken many minutes. This runs in real time.

### Using patcht to texture map a triangle mesh in matlab

Friday, August 9th, 2013

Recently I found the patcht script which lets you texture map a triangle mesh in matlab. It unfortunately does it in a brute force way: creating a textured surface for every triangle. But it’s at least something. Here’s how I use it for 2d meshes of images using the xy positions as texture coordinates:

patcht(F,V,F,[max(V(:,2))-V(:,2) V(:,1)],im);
axis equal

which produces:

We thank Scott Schaefer for providing the wooden gingerbread man image from “Image Deformation Using Moving Least Squares”.

### public alpha release of libigl – a c++ geometry processing library

Monday, April 22nd, 2013

We’ve finally released our in-house C++ library for geometry processing called libigl.

libigl is a simple c++ geometry processing library. We have a wide functionality including construction of sparse discrete differential geometry operators and finite-elements matrices such as the contangent Laplacian and diagonalized mass matrix, simple facet and edge-based topology data structures, mesh-viewing utilities for opengl and glsl, and many core functions for matrix manipulation which make Eigen feel a lot more like MATLAB.

It is first and foremost a header library. Each header file contains a single function. Most are tailored to operate on a generic triangle mesh stored in an n-by-3 matrix of vertex positions V and an m-by-3 matrix of triangle indices F.

The library may also be compiled into a statically linked library, for faster compile times with your projects.

We use the Eigen library heavily in our code. Our group prototypes a lot in MATLAB, and we have a useful conversion table from MATLAB to libigl/Eigen.

### Isointerval contour plots on triangle meshes in Matlab

Wednesday, February 29th, 2012

Matlab let’s you plot contours of scalar functions defined on grid-surfaces. But this is not easily achieved if you’re working with a triangle mesh. I had previously attempted to achieve this by first resampling the function onto a grid surface and then using Matlab’s built-in contour. This eventually (very slowly) gives you nice contour intervals and isolines separating them, but it makes the boundary of the domain ugly and doesn’t easily extend to 3d surfaces.

But now I’ve found that using the right colormap and the right rendering flags you can get pretty close to that sort of contour plot without the slow resampling step. The trick seems to be using a relatively small colormap and setting the ‘FaceLighting’ flag to ‘phong’ and the ‘FaceColor’ flag to ‘interp’. If your mesh is in V and F and you have some scalar functions in the columns of W then you can render 10 isointervals of alternating green and yellow using:

t = trisurf(F,V(:,1),V(:,2),V(:,3),W(:,5),'EdgeColor','none','FaceColor','interp','FaceLighting','phong');
colormap(repmat([1.0 0.95 0;0.1 0.5 0.2],10/2,1));
light('Position',[-1.0,-1.0,100.0],'Style','infinite');
axis equal;
view(2);

Note: the light command is not necessary to achieve sharp isointervals. But the setting the face “lighting” to ‘phong’ is necessary

And with a little effort, you can achieve something close to isointervals separated by isolines:

% number of isointervals
nin = 8;
M = jet(nin);
% thickness ratio of each isointerval compared to "isoline"
thickness = 10;
MM = reshape(permute(cat(3,repmat(M,[1 1 thickness]),zeros([nin 3 1])),[3 1 2]),nin*thickness + nin,3);
colormap(MM)
colorbar

Or on a 2D mesh without the light call: