## Archive for December, 2010

### Mesh animation with MATLAB

Thursday, December 30th, 2010

Here’s a small piece of code demonstration how to animate a triangle mesh using matlab. Assuming your mesh vertices are stored in V and your face indices in F then issue:


t = trisurf(F,V(:,1),V(:,2),V(:,3));
view(2)
axis equal
axis tight
ax = axis;
axis([ax(:,1)-10, ax(:,2)+10, ax(:,3), ax(:,4)+10]);
then = now;
while(true)
x = 10*sin(100000*(now-then));
y = abs(10*cos(100000*(now-then)));
set(t,'Vertices',[V(:,1)+x V(:,2)+y V(:,3)]);
drawnow;
pause(0.1);
end


It’s a pretty slow animation, but any faster on a 800-vertices mesh MATLAB drops frames. The pause helps but of course makes it slower, as a side effect though the my CPU is only at 20%.

### Mouse interaction with meshes in matlab

Thursday, December 30th, 2010

Here’s a small piece of code illustrating how to receive mouse clicks on a mesh plotted with trisurf in matlab. Assuming your mesh is stored with vertices in V and face indices in F, then have onmeshdown be the function you want to be called when the mesh is clicked and onvoiddown be the function you want to be called when the void in the plot outside the mesh is clicked (the white and the grey parts).


tsh = trisurf(F,V(:,1),V(:,2),V(:,3),'tag','trisurf','ButtonDownFcn',@onmeshdown);
p = get(tsh,'Parent');
set(p,'ButtonDownFcn',@onvoiddown);


Here’s a little demo showing you what you can do with that:


function dragdemo(V,F)
figure
tsh = trisurf(F,V(:,1),V(:,2),V(:,3),'ButtonDownFcn',@onmeshdown);
p = get(tsh,'Parent');
set(p,'ButtonDownFcn',@onvoiddown);
view(2);
axis equal
axis manual
down_pos = [];
drag_pos = [];
down_V = [];

function onvoiddown(src,ev)
% do nothing
end

function onmeshdown(src,ev)
down_pos=get(gca,'currentpoint');
title(sprintf('(%1.0f,%1.0f)',down_pos(1,1,1),down_pos(1,2,1)));
set(gcf,'windowbuttonmotionfcn',@drag)
set(gcf,'windowbuttonupfcn',@up)
down_V = get(tsh,'Vertices');
end

function drag(src,ev)
drag_pos=get(gca,'currentpoint');
title( ...
sprintf( ...
'(%1.0f,%1.0f) -> (%1.0f,%1.0f)', ...
down_pos(1,1,1),down_pos(1,2,1), ...
drag_pos(1,1,1),drag_pos(1,2,1)));
set( ...
tsh, ...
'Vertices', ...
[ ...
(drag_pos(1,1,1) - down_pos(1,1,1)) + down_V(:,1), ...
(drag_pos(1,2,1) - down_pos(1,2,1)) + down_V(:,2), ...
down_V(:,3) ...
]);
end

function up(src,ev)
set(gcf,'windowbuttonmotionfcn','');
set(gcf,'windowbuttonupfcn','');
end

end


### Forced draw in Chess, or forced (infinite) move cycle

Friday, December 24th, 2010

I have been pondering for a while whether it is possible to fall into a position during a game of chess from both players are forced to make them same moves over and over again. In Article 5 section 2 of the World Chess Federation’s “Laws of Chess” a game is drawn under five circumstances:

a. The game is drawn when the player to move has no legal move and his king is not in check. The game is said to end in ‘stalemate’. This immediately ends the game, provided that the move producing the stalemate position was legal.
b. The game is drawn when a position has arisen in which neither player can checkmate the opponent’s king with any series of legal moves. The game is said to end in a ‘dead position’. This immediately ends the game, provided that the move producing the position was legal. (See Article 9.6)
c. The game is drawn upon agreement between the two players during the game. This immediately ends the game. (See Article 9.1)
d. The game may be drawn if any identical position is about to appear or has appeared on the chessboard at least three times. (See Article 9.2)
e. The game may be drawn if each player has made at least the last 50 consecutive moves without the movement of any pawn and without any capture. (See Article 9.3)

So my goal is to find positions in chess in which the players end up inevitably falling working toward circumstance d or e.

These are surprisingly hard to come up with: it’s very hard to create an arrangement of the pieces such that both opponents’ moves are forced.

Today I came up with a situation where both players can move only their kings and only to two different squares, thus creating a situation where draw circumstance d must occur.

Here’s the position that forces this draw:

And here’s an animation of one possible game that could lead up to this arrangement and the inevitable forced draw:

### Delaunay3 broken in matlab?

Sunday, December 19th, 2010

I’ve been messing around with 3D tetrahedral meshes in matlab. The function delaunay3 is supposed to take a list of 3D positions and return the delaunay tessellation of those points with tetrahedra. The output looks okay at first glance using the tetramesh function. But when I started processing these tet meshes I noticed significant problems. The tets don’t actually fill the volume and seem to wildly intersect eachother. You can see this in the tetramesh plot, too:

I am asking it to delaunay tessellate a regular grid where I break the rule of no more than 4 vertices on a sphere. Figures that matlab should handle this case gracefully, but okay. Good to know now that it doesn’t. It doesn’t even return tets flattened into triangles, its just contradictory, self-intersecting tets…

It seems matlab is already on this and I just need to upgrade. The newest version of matlab has replaced delaunay3 (QHull) with DelaunayTri (CGAL).

### Convert quad mesh OBJ to triangle mesh OBJ using regular expressions (with bash and sed)

Wednesday, December 15th, 2010

Here’s a simple bash one-liner that converts a quad mesh stored in an OBJ file into a triangle mesh stored in a OBJ file. It splits every quad (a,b,c,d) into two triangles (a,b,c) and (a,c,d).

I issue this in the bash shell of my mac (you may have to tweak the sed syntax):


cat quads.obj | sed -E -e "s/f ([0-9\/]+) ([0-9\/]+) ([0-9\/]+) ([0-9\/]+)/f \1 \2 \3echo -e "\r"f \1 \3 \4/g" > triangles.obj


Or I open the file with vim and use:


:%s/f $$[0-9\/]\+$$ $$[0-9\/]\+$$ $$[0-9\/]\+$$ $$[0-9\/]\+$$/f \1 \2 \3\rf \1 \
3 \4/g


### Hide bounding box and axes in MATLAB 3D plots

Monday, December 13th, 2010

Took me a while to track this down.

To get rid of the bounding box and axes that matlab puts under your 3D plot use:


set(gca, 'visible', 'off');


### Generate “regular” tetrahedral mesh in MATLAB

Monday, December 13th, 2010

Here is a simple function which generates “regular” tetrahedral meshes in MATLAB. I put regular in quotes because while the positions of the vertices are indeed those of a regular 3D grid, the connectivity is random. This is the result of not connecting the tets myself, instead I just feed the vertex positions to delaunayn, a MATLAB function which constructs the 3D delaunay mesh from a set of points. Perhaps later I will regularize the connectivity… especially because of how slow delaunayn is.

Save this in regular_tetrahedral_mesh.m:


function [T,V] = regular_tetrahedral_mesh(nx,ny,nz)
% REGULAR_TETRAHEDRAL_MESH
%
% [T,V] = regular_tetrahedral_mesh(nx,ny,nz)
%
% Generates a "regular" tetrahedral mesh with dimensions (nx,ny,nz)
%
% Input:
%   nx  number of points in x direction on grid
%   ny  number of points in y direction on grid
%   nz  number of points in z direction on grid
% Output:
%   T  tetrahedra list of indices into V
%   V  list of vertex coordinates in 3D
%
% Example
%    [T,V] = regular_tetrahedral_mesh(3,3,3);
%    tetramesh(T,V); %also try tetramesh(T,V,'FaceAlpha',0);
%
%

% Create a grid
[x,y,z] = meshgrid(linspace(0,1,nx),linspace(0,1,ny),linspace(0,1,nz));
V = [x(:) y(:) z(:)];
if(nx==2 && ny==2 && nz==2)
warning(['delaunayn vomits on 2x2x2...' ...
V = [V;0.5,0.5,0.5];
end
T = delaunayn(V);
end


Then you can call it like this:


[T,V] = regular_tetrahedral_mesh(3,3,3);
tetramesh(T,V);


which should produce this:

Update: Here’s the real regular tetrahedral mesh generator. Complete with regular connectivity. Next step is to keep track of surface faces…



function [T,V] = regular_tetrahedral_mesh(nx,ny,nz)
% REGULAR_TETRAHEDRAL_MESH
%
% [T,V] = regular_tetrahedral_mesh(nx,ny,nz)
%
% Generates a regular tetrahedral mesh with dimensions (nx,ny,nz)
%
% Input:
%   nx  number of points in x direction on grid
%   ny  number of points in y direction on grid
%   nz  number of points in z direction on grid
% Output:
%   T  tetrahedra list of indices into V
%   V  list of vertex coordinates in 3D
%
% Example
%    [T,V] = regular_tetrahedral_mesh(3,3,3);
%    tetramesh(T,V); %also try tetramesh(T,V,'FaceAlpha',0);
%
%

% Create a grid
[x,y,z] = meshgrid(linspace(0,1,nx),linspace(0,1,ny),linspace(0,1,nz));
V = [x(:) y(:) z(:)];
% meshgrid flips x and y ordering
idx = reshape(1:prod([ny,nx,nz]),[ny,nx,nz]);
v1 = idx(1:end-1,1:end-1,1:end-1);v1=v1(:);
v2 = idx(1:end-1,2:end,1:end-1);v2=v2(:);
v3 = idx(2:end,1:end-1,1:end-1);v3=v3(:);
v4 = idx(2:end,2:end,1:end-1);v4=v4(:);
v5 = idx(1:end-1,1:end-1,2:end);v5=v5(:);
v6 = idx(1:end-1,2:end,2:end);v6=v6(:);
v7 = idx(2:end,1:end-1,2:end);v7=v7(:);
v8 = idx(2:end,2:end,2:end);v8=v8(:);
T = [ ...
v1  v3  v8  v7; ...
v1  v8  v5  v7; ...
v1  v3  v4  v8; ...
v1  v4  v2  v8; ...
v1  v6  v5  v8; ...
v1  v2  v6  v8];
end


By the way, for a 20 by 20 by 20 grid the updated code is about 200 times faster than the code which calls delaunayn. For bigger grids it will only get better as connecting the tets procedurally is linear and vectorized and delaunayn uses QHull which is in general O(n*log(n)) and I’m not sure how parallelized it is.

I generate this 10 by 10 by 10 tet mesh in around a hundredth of a second (rendering it with tetramesh took about 4 seconds )

### cudaMemcpyToSymbol not copying data from host to __constant__ memory

Friday, December 10th, 2010

I had declared some constant memory variable at the top of my program like this:


__device__ __constant__ float constant_T[M];


Then I tried to fill it from my host array T with the following:


cudaMemcpyToSymbol("constant_T", T, sizeof(float)*M,cudaMemcpyHostToDevice);


Surprisingly I got no compilation errors and no cuda errors at runtime via


printf("Cuda errors: %s\n", cudaGetErrorString( cudaGetLastError() ) );


After some staring I finally figured it out, the copying doesn’t need the last parameter. It should be:


cudaMemcpyToSymbol("constant_T", T, sizeof(float)*M);


### Triangle mesh filled contour in MATLAB

Monday, December 6th, 2010

Here is matlab snippet that takes a 2D triangle mesh and a function and produces a pretty, filled contour plot in MATLAB. If your mesh vertices are stored in V, and your triangles in F, then use my previously posted code to get the boundary/outline edges of your mesh in O. Then you can create a contour plot of some function W defined over V using:


% resample W with a grid, THIS IS VERY SLOW if V is large
[Xr,Yr,Wr] = griddata(V(:,1),V(:,2),W,unique(V(:,1)),unique(V(:,2))');
% find all points inside the mesh, THIS IS VERY SLOW if V is large
IN = inpolygon(Xr,Yr,V(O,1),V(O,2));
% don't draw points outside the mesh
Wr(~IN) = NaN;
contourf(Xr,Yr,Wr)


And for a mesh like this one:

You get something like this:

You can also approach a continuous contour with this:


contourf(Xr,Yr,Wr,200)


Compare this to what you can get from the much faster:


trisurf(F,V(:,1),V(:,2),W,'FaceColor','interp','EdgeColor','interp')
view(2)


which produces:

The only problem with this method is that the boundary looks a little nasty because of the resampling. Nothing that can’t be fixed quickly in photoshop… Or even in MATLAB with something like:


line(V([O(:);O(1)],1),V([O(:);O(1)],2),'Color','k','LineWidth',6)


Which puts an ugly thick outline around the contour, but makes it look a little better…

Source

### Find outline of triangle mesh and plot in MATLAB

Monday, December 6th, 2010

Here’s a little snippet to determine the edges along the boundary of a triangle mesh in matlab. Given that your vertex values are in V and your triangle indices are in F, this will fill O with a list of edges make up the outline of your mesh:


% Find all edges in mesh, note internal edges are repeated
E = sort([F(:,1) F(:,2); F(:,2) F(:,3); F(:,3) F(:,1)]')';
% determine uniqueness of edges
[u,m,n] = unique(E,'rows');
% determine counts for each unique edge
counts = accumarray(n(:), 1);
% extract edges that only occurred once
O = u(counts==1,:);


Then you can view the outline with, in 2D:


plot([V(O(:,1),1) V(O(:,2),1)]',[V(O(:,1),2) V(O(:,2),2)]','-')


and in 3D:


plot3([V(O(:,1),1) V(O(:,2),1)]',[V(O(:,1),2) V(O(:,2),2)]',[V(O(:,1),3) V(O(:,2),3)]','-')