## Posts Tagged ‘gptoolbox’

### Quad meshing in matlab using gptoolbox

Friday, February 17th, 2017

This morning I played around with a very simple way of quad meshing a triangle mesh. This is more or less the “ideal pipeline” for quad meshing via parameterization. It will likely only work for very simple meshes. I’m pleasantly surprised that it works at all.

% Load a mesh with boundary and no topological (donut) holes
V = V*axisangle2matrix([1 0 0],-pi/2);
% Construct the LSCM matrix
[~,Q] = lscm(V,F,[],[]);
M = massmatrix(V,F);
% Use Fiedler vector as parameterization
[EV,ED] = eigs(Q,repdiag(M,2),5,'sm');
U = reshape(EV(:,end-3),[],2);
% try to find "canonical" rotation
[sU,sS,sV] = svd(U'*U);
U = U*sU;
% Lay down a grid in 2D
h = ((max(U(:,1))-min(U(:,1)))/(160-2));
[X,Y] = meshgrid(min(U(:,1))-h:h:max(U(:,1))+h,min(U(:,2))-h:h:max(U(:,2))+h);
% Find all quads that are strictly inside the parameterized mesh
[QQ,QU] = surf2patch(X,Y,0*Y);QU = QU(:,1:2);
[sqrD,I,C] = signed_distance([QU QU(:,1)*0],[U U(:,1)*0],F);
QQ = QQ(all(abs(sqrD(QQ))<1e-10,2),:);
[QU,I] = remove_unreferenced(QU,QQ);
QQ = fliplr(I(QQ));
% Snap boundary vertices of quad mesh to boundary of mesh, and
% solve little Dirichlet problem to diffuse the displacements smoothly
QF = [QQ(:,[1 2 3]);QQ(:,[1 3 4])];
L = cotmatrix(QU,QF);
b = unique(outline(QQ));
[~,~,QU(b,:)] = point_mesh_squared_distance(QU(b,:),U,outline(F));
% Locate each quad mesh vertex in the parameterized mesh
[sqrD,I,C] = signed_distance([QU QU(:,1)*0],[U U(:,1)*0],F);
B = barycentric_coordinates(C(:,1:2),U(F(I,1),:),U(F(I,2),:),U(F(I,3),:));
% Interpolate 3D positions at quad mesh vertex locations
QV = V(F(I,1),:).*B(:,1) + V(F(I,2),:).*B(:,2) + V(F(I,3),:).*B(:,3);


To create the image above use:

clf;
hold on;
tq = tsurf(QQ,QV,'FaceVertexCData',repmat(blue,size(QV,1),1),'EdgeColor','none',fsoft,fphong);
tt = tsurf(F,V-[0.8 0 0],'FaceVertexCData',repmat(1-0.8*(1-blue),size(V,1),1),'EdgeAlpha',0.5,fsoft,fphong);
axis equal;
camproj('persp');
view(-38,17);
drawn;
camlight;
plot_edges(QV,[QQ(:,2:3);QQ(:,[4 1])],'Color',orange);
plot_edges(QV,[QQ(:,1:2);QQ(:,3:4)],'w');hold off;
l = light('Position',[-10 0 10],'Style','infinite');
set(gca,'Visible','off');
set(gcf,'Color','w');


There’s a lot to improve about the quad mesh: the boundaries are not aligned with the meshing direction, the meshing directions are not always aligned with principle curvature, etc. Not bad for 30 mins of coding using tools that are just lying around.

### Plot a bunch of edges in matlab with per-edge color data

Monday, December 19th, 2016

Suppose you have a list of vertices V (#V by dim) and a list of edges E (#E by 2) indexing V, then one hacky to draw all of these edges with per-vertex color data DV (#V by 1) is to use trisurf:

trisurf(E(:,[1 1 2]),V(:,1),V(:,2),V(:,3),'CData',DV,'EdgeColor','interp');


However, if you try to use this trick to draw per-edge color data DE (#E by 1), you’ll find that 'EdgeColor','flat' does not work. Instead you can explode you compact “mesh” of edges into individual edges and repeat you edge color data for each edge point:

trisurf(reshape([1:size(E,1),1:size(E,1)*2],[],3),V(E(:),1),V(E(:),2),V(E(:),3),'CData',repmat(DE,2,1),'EdgeColor','interp');


### Mesh Arrangements for Solid Geometry preprint

Thursday, April 21st, 2016

I’m very excited to publish my work with Qingnan Zhou, Eitan Grinspun and Denis Zorin on mesh booleans at this year’s SIGGRAPH. The paper is titled Mesh Arrangements for Solid Geometry. The code is (and has been) in libigl and gptoolbox and pymesh.

### Variadic mesh boolean operations in libigl and gptoolbox

Thursday, April 14th, 2016

I’ve just pushed some new changes to libigl and gptoolbox to expose the variadic implementation of our robust mesh boolean operations. By variadic, I mean that the boolean function can take one mesh as input or two meshes as input (usual binary case) or three or four and so on. This means you can easily take the union or intersection of n objects, but it also means that complex operations can be reduced to single call. For example, identifying regions inside at least k out of n objects. I first became aware of this variadic concept when reading “QuickCSG: Arbitrary and Faster Boolean Combinations of N Solids”.

In libigl, you can call igl::copyleft::cgal::mesh_boolean with an n-long list of mesh vertex arrays and an n-long list of mesh face arrays. Rather than the usual union, intersection, minus etc., you can also pass a C++ function handle. This function handle will filter the winding number vector at any point in space. For example, if you’d like to extract the region inside object 0 and not inside object 1 and either inside object 2 or not inside object 3. The filter would return w(0) > 0 && w(1)<=0 && (w(2)>0 || w(3)<= 0).

After a bit of pointer tweaking, I’ve also exposed this interface to the matlab wrapper in gptoolbox. You can pass a matlab function handle after the optional argument ‘WindingNumberFilter’. For example, assuming V and F are cell arrays containing vertex position arrays and face indices respectively, the following command will extract the region inside at least 3 of the inputs:

 [VV,FF,J] = mesh_boolean(V,F,'','WindingNumberFilter',@(w) sum(w>0)>=3);


You could also call it with

 [VV,FF,J] = mesh_boolean(V{1},F{1},V{2},F{2},V{3},F{3},V{4},F{4},'','WindingNumberFilter',@(w) sum(w>0)>=3);


Here’s an example of extracting min-k results from 4 spheres:

Since extraction from our mesh arrangement is cheap, each operation takes the same amount of time. And fortunately it seems that the mex function handle overhead is not so bad.

Update: Variadic operations are also useful for condensing entire binary CSG trees: (A union B) minus (C union D) could be a single extraction function on (A,B,C,D) return (w(0)>0 || w(1)>0) && !(w(2)>0 || w(3)>0)

### 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.

### Conservative voxelization in gptoolbox

Tuesday, November 3rd, 2015

A while ago I implemented a function to voxelize a given triangle mesh in matlab as part of gptoolbox. Rather than use 3D-Bresenham-style rasterization, I needed a conservative voxelization: where every point (not just vertex) on the mesh is contained inside the voxelization. In other words, output all voxels that intersect the triangle mesh.

Here’s a little snippet to produce the following meshes:

for k = [10 20 40 80]
[DV,I] = remove_unreferenced(DV,Q);
Q = I(Q);
d = render_in_cage(V,F,DV,Q,'ColorIntersections',false);
d.tc.LineWidth = 2;
end


### Mesh decimation (aka simplification) in matlab

Thursday, October 15th, 2015

I’d forgotten that I’d discovered that matlab has a built-in function for triangle mesh simplification: reduce_patch. Here’s a little comparison between this method and the ones I’ve wrapped up in gptoolbox.

### Input mesh

Given an input mesh with vertices V and faces F (this elephant has 14732 faces):

### Built-in matlab reduce patch

We can decimate the mesh to 1000 faces using pure built-in matlab:

[mF,mV] = reducepatch(F,V,1000);


Notice that it’s a rather adaptive remeshing. There aren’t (m)any options for controlling reducepatch, just the 'fast' flag, but in this case it will produce an identical output.

### Libigl’s decimate

Libigl has a framework for edge-collapse decimation on manifold meshes and a default setting for regular meshing:

[iV,iF] = decimate_libigl(V,F,1000);


### CGAL’s decimation

I have a wrapper for CGAL’s decimation algorithms. These can result in a regular remeshing with 1000 faces:

[c1V,c1F] = decimate_cgal(V,F,1000/size(F,1));


[c2V,c2F] = decimate_cgal(V,F,1000/size(F,1),'Adaptive',true);


CGAL has a fairly general interface for an edge-collapser, but these are the only metrics provided by default from CGAL.

### Timing

Suppose my original elephant had 3.7 million faces instead of just 14K, how do these methods measure up when reducing to 1000 faces:

Method Time
reducepatch 64.4s
libigl 74.56s
cgal (regular) 117.7s
qslim 108.0s + 507.5s + 59.4s

My qslim wrapper is currently very silly. It’s calling qslim as an executable and then reading in the entire collapse log and conducting the collapse within matlab (the second time is for reading the log, the third for conducting the collapse).

None of these methods guarantee that the result is self-intersection free. I’ve noticed that qslim and matlab will create non-manifold meshes from manifold inputs. The current libigl implementation should not creative non-manifold output, but assumes the input is a closed surface. For casual decimation, it seems like matlab’s reducepatch and libigl’s decimate_libigl are the best bets for speedy adaptive and regular meshing respectively. Conveniently they also require the least in terms of dependencies.

### Bounded biharmonic weights demo binary

Thursday, August 13th, 2015

I’ve finally gotten around to creating a stand-alone mac os x app for the prototype I use to give demos of our bounded biharmonic weights. The only catch is that you’ll have to grab a mosek (free academic trial) license and put it in the right place

Unfortunately do to license reasons I cannot release the source code for this app, but the code for bounded biharmonic weights is in my matlab gptoolbox and in our C++ libigl

### 3D printable stencil from binary image

Tuesday, July 21st, 2015

Using code lying around in gptoolbox it was easy to whip up a little function that converts a binary image into a 3d-printable stencil.

I’ve added a new function bwstencil.m to gptoolbox. Here’re the simple contents:

  imh = imfill(im,'holes');
im = imresize(im,2,'nearest');
[V,F] = bwmesh(~im,varargin{:});
[V,F] = extrude(V,F);
[V,F] = remesh_planar_patches(V,F);


Notice that I’m removing “islands” (regions that would be part of the stencil but not connected to the boundary). You can see that the eye of Hans Hass above has been removed in the stencil.

For this example, I’ve 3D printed the result as a proof of concept:

In retrospect, for this example I should have created the opposite stencil if I’m painting black on white:

and the physical proof:

### Wave equation on surfaces

Tuesday, July 21st, 2015

Inspired by a course Misha Kazhdan gave at SGP this year, I started playing around with the wave equation on surfaces. Below is a simple matlab demo:

The core formula here is to update the “height” H_i (extrusion along normal) at vertices for each time step i based on the last two steps:

((1 + b * dt) * M - a * dt^2 *L) H_i = (M * (( 2 - b * dt) * H_{i-1} - H_{i-2}));


b is the dampening factor, dt is the time step size, M is the mass matrix, 1/sqrt(a) is the speed of the wave, L is the cotangent Laplacian, H_i is the unknown heights, H_{i-1} and H_{i-2} are the heights at the previous two time steps. This is discretization with implicit time stepping of the wave equation:

∂²/∂t² u = a∆u - b ∂u/∂t


There’s surely a more canonical place to find this, but I liked the description in “Fast and Exact (Poisson) Solvers on Symmetric Geometries” [Kazhdan 2015].

Save in wave_demo.m

function wave_demo(V,F)
V = bsxfun(@minus,V,mean(V));
V = V./max(abs(V(:)));
V(:,end+1:3) = 0;
N = per_vertex_normals(V,F);
L = cotmatrix(V,F);
M = massmatrix(V,F);
M = M / max(abs(M(:)));

dt = 0.0075;
% 1/sqrt(a) is the speed of the wave
a = 1;
H = zeros(size(V,1),1);
% dampening
b = 0.05;
A = ((1+b*dt)*M - a*dt^2*L);
% Pre-factor
[cL,p,cS] = chol(A,'lower');
cU = cL';
cP = cS';
cQ = cS;
chol_solve = @(X) cQ * (cU \ (cL \ ( cP * X)));
G = H;
G_m1 = G;

[VV,FF,S] = loop([V+bsxfun(@times,G,N) H],F);
tsh = trisurf(FF,VV(:,1),VV(:,2),VV(:,3),  ...
'FaceColor','interp','FaceLighting','phong', ...
'EdgeColor','none','CData',VV(:,4));
tsh.SpecularStrength = 0.2;
tsh.SpecularExponent = 100;
tsh.ButtonDownFcn = @ondown;

CM = parula(30);colormap(CM(end-11:end,:));
set(gcf,'Color',[0.3 0.4 1.0]);
set(gca,'Visible','off');
l = light('Position',[-0.4 -0.4 1.0],'Style','infinite');
axis equal;
axis manual;
expand_axis(1.2);
view(2);
caxis([-max(abs(H)) max(abs(H))]);
T = get(gca,'tightinset');
set(gca,'position',[T(1) T(2) 1-T(1)-T(3) 1-T(2)-T(4)]);

t = 0;
draw_t = 0;
tic;
while true
G_m2 = G_m1;
G_m1 = G;
G = chol_solve(M*((2-b*dt)*G_m1 - G_m2 + dt*0));

t = t+dt;
draw_t = draw_t+toc;
if draw_t > 0.033
VV = S*[V+bsxfun(@times,G,N) G];
tsh.Vertices = VV(:,1:3);
tsh.CData = VV(:,4);
title(sprintf('t: %g, max(H): %g',t,max(abs(G))),'FontSize',15);
drawnow;
tic;
draw_t = 0;
end
end

function ondown(src,ev)
H = a/10*mean(diag(M))*(M\(S'*sparse( ...
knnsearch(VV(:,1:3),ev.IntersectionPoint,'k',2),1,[1;-1],size(VV,1),1)));
G = G+H;
end
end