## Posts Tagged ‘gptoolbox’

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


### Physical simulation “skinning” with biharmonic coordinates

Thursday, June 25th, 2015

This is a cute little example we didn’t have time to put into our upcoming SIGGRAPH paper Linear Subspace Design for Real-Time Shape Deformation.

In this example I first decimate the high resolution (blue) octopus to create the low-res orange octopus. I’m careful to use a decimate that keeps vertices at (or at least near) their original positions. Then I compute “biharmonic coordinates” for vertices of the blue mesh with respect to the vertices of the coarse orange mesh. Notice that the orange mesh does not need to be an enclosing cage. Rather its points just need to lie in or on the blue octopus (there’re even ways to deal with points outside the blue octopus, but right now I’m skipping that). These weights are computed in gptoolbox with:

b = knnsearch(orange_V,blue_tet_V);
W = biharmonic_coordinates(blue_tet_V,blue_tet_T,b);


Then I run a sort of co-rotational linear-elasticity simulation on the orange mesh with collisions activated for the floor. The blue mesh tracks the orange mesh by simply interpolating the orange vertex positions via the coordinates:

blue_V = W * orange_V;


Simple as that. Smooth physics “skinning” with a matrix-matrix multiplication. This is due to the affine precision of our coordinates and their squared Laplacian energy minimization. The “standard” approach to physics “skinning” would be to create the orange tet-mesh so that it envelops the blue mesh and then use piecewise-linear (i.e. non-smooth) interpolation to move the blue mesh vertices.

Update: I should make it super clear that this is not a subspace (aka model reduction) simulation. The simulation of the orange mesh is not aware of the blue mesh at all. Sometimes subspace physics would be more appropriate, though other times pushing the subspace through all parts of the simulation pipeline is impractical/impossible. Hence this prototypical result.

For posterity here’s the entire example code:

% load hi-res mesh
[~,CV,CF] = qslim(uV,uF,900,'QSlimFlags','-O 0');
[CV,CF] = meshfix(CV,CF);
% Find closest points
b = knnsearch(uV,CV);
CV = uV(b,:);

% Tet-mesh fine resolution mesh and compute weights
[dV,dT,dF] = tetgen(uV,uF,'Flags','-q2');
W = biharmonic_coordinates(dV,dT,b);
W = W(1:size(uV,1),:);

[VV,VT,VF] = tetgen(CV,CF,'Flags','-q100 -Y');
off = mean(VV);
VV=bsxfun(@minus,VV,off);
UU = bsxfun(@plus,VV,[0 0 1]);
vel = zeros(size(UU));
data = [];
M = massmatrix(VV,VT);
% This normalization is very important
M = M/max(M(:));
dt = 1e-2;

t = tsurf(VF,UU, ...
'FaceColor',[0.8 0.5 0.2], ...
'SpecularStrength',0,'DiffuseStrength',0.3,'AmbientStrength',0.7);
%t.Visible = 'off';
l = light('Style','local','Position',[5 0 20]);
l2 = light('Style','infinite','Position',[-1 -1 2]);
hold on;
up = @(UU) bsxfun(@plus,[1 0 0],W*UU(1:numel(b),:));
u = tsurf(uF,up(UU),'EdgeColor','none', ...
'FaceVertexCData',repmat([0.3 0.4 1.0],size(uV,1),1), ...
fphong,'SpecularStrength',0,'DiffuseStrength',0.3,'AmbientStrength',0.7);
e = -1e-4;
QV = [-0.8,-0.5 e; -0.8  0.5 e; 1.5  0.5 e;1.5,-0.5 e];
QQ = [1 2 3 4];
trisurf(QQ,QV(:,1),QV(:,2),QV(:,3),'EdgeColor','none','FaceColor',[0.7 0.7 0.7]);
st.FaceColor = [0.5 0.5 0.5];
su.FaceColor = [0.5 0.5 0.5];
set(gca,'Visible','off');
set(gcf,'Color','w');
end
hold off;
axis equal;
axis([-0.8 1.8 -0.5 0.5 0 1.5]);
axis manual;
camproj('persp');
T = get(gca,'tightinset');
set(gca,'position',[T(1) T(2) 1-T(1)-T(3) 1-T(2)-T(4)]);

while true
dur = 0;
UU0 = UU;
[UU,data] = arap(VV,VT,[],[], ...
'Energy','elements', ...
'V0',UU0, ...
'Data',data, ...
'Dynamic',M*repmat([0 0 -9.8],size(VV,1),1), ...
'MaxIter',100, ...
'TimeStep',dt, ...
'Vm1',UU-vel);
vel = UU-UU0;
C = UU(:,3)<0;
UU(C,3) = -UU(C,3);
vel(C,3) = -vel(C,3);
dur = dur+dt;
t.Vertices = UU;
u.Vertices = up(UU);