Graph theory toolbox
Copyright (c) 2007 Gabriel Peyre
This toolbox contains useful functions to deal with graph and triangulation.
The basic representation of a graph of n vertices is the adjacency matrix A where A(i,j)=1 if vertex i is linked to vertex j. A graph often comes with a geometric realization in R^d which an (d,n) matrix where vertex(:,i) is the position of the ith vertex.
A triangulation of m faces and n vertex is represented through:
* a set of faces which is a (3,m) matrix where face(:,i) are the vertices indexes of the ith face.
* a set of vertex which is a (d,n) matrix.
The toolbox contains function to deal more easily with a triangulation data structure, and allows to retrieve vertex and face 1ring and switch from adjacency to faces.
The graph part of the toolbox contains function to creates synthetic graph and compute shortest path (dijkstra and isomap algorithm).
This toolbox contains a lot of function to deal with spectral theory of triangulation. You can load triangulations from files and then display the resulting mesh. It allows to compute various laplacian operator, and the to compute parameterization using spectral decomposition, harmonic mapping, free boundary harmonic mapping, and isomap.
1.4  Modified license. 

1.2  Update of Licence 

1.1  n/a 

Added html help files. 

Added image approximation using refinement. 

Added support for tetrahedral meshes. 

Fixed a few bugs. 

Added curvature computation on meshes. 

Added heat diffusion on meshes, clean definition of laplacians. 

GPL license. 

Added free boundary parameterization. 

Added support for various file formats (VRML,PLY,SMF).


Now in zip format. 

Added dual graph support. 
Inspired: gptoolbox, simulation of breast cancer, Code to realize Olga Sorkine paper, ICP Registration using Efficient Variants and MultiResolution Scheme, Optimal Step Nonrigid ICP, laplacian_surface_editing_3D, Local Depth SIFT and Scale Invariant Spin Image local features for 3D meshes
filipp (view profile)
Hi there, even though compute_curvature is very easy to implement it seems that its output is not completely reliable. I computed curvature for a cube and get different curvature values for parallel edges of the cube  even though they are exactly the same (just translated in space). Other curvature algorithms (that are however showing other problems) produce the same values for these parallel edges. Any idea why this happens?
Thanks!
Berthe Emmanuelle Kamga Monayoum (view profile)
i get this error by vtk function ; 'Warning: Problem in reading vertices.'
Suwoong Heo (view profile)
Cool! Can someone recommend me a textbook or paper related to the theory of this toolbox?
Thomas AttaFosu (view profile)
Here's a workaround to the error in 'compute_curvature' at line 75: "dp = sum( normal(:,E(:,1)) .* normal(:,E(:,2)), 1 );
Copy and paste the following code into 'compute_curvature' right after line 47, i.e after the line "s = [1:m 1:m 1:m];" and right before "A = sparse(i,j,s,n,n)"
[u,ui,~] = unique([i' j'],'rows','stable');
s = s(:,ui);
i = u(:,1);
j = u(:,2);
isfzade (view profile)
Sébastien Herbert (view profile)
Hi, thanks for your submission.
I think there is a small bug in the read_ply.m function. When I run it, I get this error
Reference to nonexistent field 'vertex_indices'.
Error in read_ply (line 16)
vi = d.face.vertex_indices;
It appears several field names are possible line 460 but not tested afterward
PropertyName = {'vertex_indices','vertex_indexes','vertex_index','indices','indexes'};
In my case, the error could easily be solved by changing the line 26 to
if isfield(d.face,'vertex_indices')
vi = d.face.vertex_indices;
elseif isfield(d.face,'vertex_index')
vi = d.face.vertex_index;
end
Sébastien Herbert (view profile)
Jinny (view profile)
For read_vtk function, what should be for 'verbose' input? If I don't put any value for I gets 'Warning: Problem in reading vertices.'
chao lei (view profile)
@az why i modify the code using your method,but the compute_normal made some wrong?
Maria Charalampidou (view profile)
marc (view profile)
Seems good !
Jingwei Zhuang (view profile)
lilian yang (view profile)
Thanks for this library.
But I met question when running the following code:
options.method = 'flattening';
xy2 = compute_parameterization(vertex,faces,options);
options.method = 'isomap';
xy3 = compute_parameterization(vertex,faces,options);
% display
clf;
subplot(1,2,1);
plot_graph(A,xy2,'k.'); axis tight;
title('Laplacian eigenmaps');
subplot(1,2,2);
plot_graph(A,xy3,'k.'); axis tight;
title('Isomap');
The error is:
The subscript index must be a positive integer type or a logical type.
Error perform_dijkstra_propagation_slow>dijkstra_init (line 110)
data.A(start_verts) = 0;
Error perform_dijkstra_propagation_slow (line 42)
data = dijkstra_init(W, start_verts);
Error perform_dijkstra (line 66)
[D,S] = perform_dijkstra_propagation_slow(W',start_points1,end_points1,nb_iter_max, H); % use transposate
Error compute_distance_graph (line 34)
[d,S] = perform_dijkstra(W, point_list(i));
Error isomap (line 101)
D = compute_distance_graph(Ws, points_list);
Error compute_parameterization (line 78)
vertex1 = isomap(D,ndim, options);
lilian yang (view profile)
Zhenyi Wang (view profile)
Alexander Naitzat (view profile)
fxy (view profile)
az (view profile)
I solved few bugs (out of bounds and divide by zero errors) that were in "compute_curvature".
Copy the following lines instead of lines 4271:
% associate each edge to a pair of faces
i = [face(1,:) face(2,:) face(3,:)];
j = [face(2,:) face(3,:) face(1,:)];
L = sub2ind([n, n], i, j);
[a,b]=hist(L,unique(L));
Q = ismember(L,b(a>1));
[i1, j1] = ind2sub([n, n], L(Q));
W = [i' j'];
R = ismember(W, [i1' j1'], 'rows');
J = find(R);
[Q, K] = sortrows(W(R,:));
J = J(K);
[~, ia, ~] = unique(Q, 'rows');
R = true(size(J));
R(ia) = false;
J = J(R);
s = [1:m 1:m 1:m];
s(J) = 0;
A = sparse(i,j,s,n,n); % if (i,j) are an edge at some face then A(i,j) is
% this face index.
[i,j,~] = find(A); % direct link
X = [i j];
Y = [j i];
[~, I1, ~]=intersect(X, Y, 'rows');
I = X(I1, :);
K1 = I(I(:,1) < I(:,2),:); % only directed edges
K2 = [K1(:,2) K1(:,1)];
[~, ~, s1]= find(A(sub2ind(size(A),K1(:,1), K1(:,2))));
[~, ~, s2]= find(A(sub2ind(size(A),K2(:,1), K2(:,2))));
% links edge>faces
E = [s1 s2];
i=K2(:,1);
j=K2(:,2);
ne = length(i); % number of directed edges
% normalized edge
e = vertex(:,j)  vertex(:,i);
d = sqrt(sum(e.^2,1));
d(d<eps) = 1;
e = e ./ repmat(d,3,1);
% avoid too large numerics
d = d./mean(d);
hariharan mahalingam (view profile)
Daniel Frisch (view profile)
The test files don't run due to changed function names or syntax, as well as missing functions or input files. This is sad because these test files could give the users a quick overview about the capability of the toolbox.
Teng Long (view profile)
Hi I tried to plot a mesh with its normals by 'plot_mesh(vertex,faces,options);' It returns an error saying"plot_mesh (line 145) Z and U must be the same size". I realized that the normals vector returned from compute_normal functions is of size 3*n. Currently I modified line 145 in plot_mesh function to "quiver3(vertex(sel,1),vertex(sel,2),vertex(sel,3),normal(1,sel)',normal(2,sel)',normal(3,sel)',normal_scaling);" and it worked. I'm just wondering why you have to return a list of normals from compute_normal function instead of a vector of normal?
Teng Long (view profile)
Duc Fehr (view profile)
Ghazal (view profile)
Super great Toolbox!
I need to have colored meshes on Matlab. I tried to read an object file using read_obj.m but don't succeed to plot it.
Do you have any idea how I can have a 3D mesh with a color value assigned to each face on Matlab?
iMorrun (view profile)
Pretty nice code.but I still met question when running the code.when I wanted to read a .ply file by using read_ply function,the system always gave error like this:error in read_ply at 16 vi = d.face.vertex.indices;...I am in pretty hurry,can any one help me solve this problem?By the way,does anyone know the website which can download .off file?
Ego (view profile)
Fritz (view profile)
@az: test_curvature.m works fine for me with the 'mushroom', see results:
http://imgur.com/a/OA4xY
Probably it's related to your mesh? Can you upload an example?
az (view profile)
Has anyone solved the bug in "compute_curvature"  at line 75?
Shahnawaz (view profile)
I solved this problem :) One more issue. Can I read .ply file as mesh instead of .off file?
Broni (view profile)
Hi! thanks for the nice toolbox.
I want to create a graph for HD image (1080x1920), but the build_graph... returns error because it goes out of memory when initializing the very large graph matrix.
Do you maybe have any suggestions how to tackle this issue. I want to implement a gaphcuts algorithm for HD images or larger resolutions.
EL OIRRAK (view profile)
How can iniatilize OPTIONS ???
kijo maken (view profile)
to fix the problem with quiver at plot_mesh, you just need to scroll to line 139 and change the code to quiver3(vertex(sel,1),vertex(sel,2),vertex(sel,3),normal(1,sel)',normal(2,sel)',normal(3,sel)',normal_scaling);
Richard Kennaway (view profile)
The function select3d contained in this toolbox no longer works in 2014b and 2015a (prerelease). It contains a function local_Data2PixelTransform, which tries to obtain the properties 'x_RenderTransform', 'x_RenderOffset', and 'x_RenderScale' of an axes object. These properties no longer exist.
EL OIRRAK (view profile)
how to use compute_curvature
Helen (view profile)
As others have pointed out, the compute_curvature function results in an index out of bounds error in some cases. So far, my untested and kind of random 'fix' was to comment out the code
%normal = normal'; (line 72) and
insert it again as
normal = normal'; (at line 146, after compute_normal()).
This seems to avoid the out of index error BUT I have no idea if it a right thing to do otherwise.
Some support from Gabriel, would be highly appreciated.
Thanks a lot for your efforts
(on line 75 "dp = sum( normal(:,E(:,1)) .* normal(:,E(:,2)), 1 );")
Alec Jacobson (view profile)
Thanks for this library.
Just spent quite some time tracking down the same bug in `write_ply.m` as Chris Rorden found a year ago. Please fix it.
dhara (view profile)
great tool.. thnx...
but i m able to read this .wrl file cara1_abajo.wrl.. plz anybody help me.. it's urgent... i tried a lot but getting error in reshape command
David (view profile)
It seems to be a useful tool.
There might be a bug to generate ply files. I there is always a bug as follow.
Error using write_ply (line 22)
vertex does not have correct format.
However, I have strictly follow the input format for using this function. The vertex is 'nb.vert x 3' array, can you help me about that?
Xu (view profile)
Good tools.
HAYOUNG (view profile)
Many thanks to this great toolbox.
FYI, he also made comprehensive manual in his website (https://www.ceremade.dauphine.fr/~pey re/matlab/graph/content.html), which is also nice.
HAYOUNG (view profile)
Many thanks to great toolbox.
ness (view profile)
Hi
Thanks for your code.
I have a question when i use plot_mesh i can visualize the 3D mesh, can i hold on the same figure a 2d image. when i use 'hold on' then 'imshow' it doesn't work.
Ben (view profile)
anyone solved the bug in the compute_curvature function (on line 75 "dp = sum( normal(:,E(:,1)) .* normal(:,E(:,2)), 1 );")? Thanks!
Anass (view profile)
Hi,
I'm asking if the function " red_obj()" can read a model with texture?
I tried but it did not be able to read it.
Can you help please?
David (view profile)
There is some code to display geometry images. Does this toolbox provide code for producing a geometry image from a mesh?
Chris Rorden (view profile)
Lots of nice code. I did find a small bug with write_ply when creating files with more than 32767 faces. To fix this you need to change the lines
IntegerDataMin = [128,0,2^15,2^31,0];
IntegerDataMax = [127,255,2^161,2^311,2^321];
to read
IntegerDataMin = [128,0,2^15,0,2^31,0];
IntegerDataMax = [127,255,2^151,2^161, 2^311,2^321];
Eitan (view profile)
Very useful tool.
I noticed that compute_curvature tries to enforce Cmin<=Cmax. In case of negative Cmax, it can cause abs(Cmin)>abs(Cmax), which I think is not right. IMHO these lines can be removed from the function.
wikimenWoko wikk (view profile)
Julian Xue (view profile)
is there a method here to just plot an adjacency matrix?
Ged Ridgway (view profile)
Looks really useful. However, isomap seems to fail:
>> xy = isomap(dat);
??? Undefined command/function 'compute_nn_graph'.
Error in ==> isomap at 85
D = compute_nn_graph(X,options);
Peter (view profile)
kamel (view profile)
Very useful
compute_curvature will generate an error on line 75 ("dp = sum( normal(:,E(:,1)) .* normal(:,E(:,2)), 1 );") for SOME surfaces.
what's the solution of this bug
thanks
Raymond Cheng (view profile)
Thanks for your sharing.
YUANCHUN WANG (view profile)
Very useful.
But I have a question: when I tried to use the function "perform_mesh_simplification" to make the simplificatino for mesh. Matlab said that QSlim is not a command. And here is the command where situates the problem: "system(['QSlim tmp.smf o tmp1.smf s ' num2str(nface)]);"
I'm afraid this is because the qslim.exe was missed?
thanks for your explaination!
europe lee (view profile)
test_diffusion_wavelet.m file did not work!
Because it miss some functions, such as MakeDiffusion function.
Although I download and install diffusion wavelet code from MAURO MAGGIONI, I find that MakeDiffusion in his code is different from yours,like parameter called 'gauss' in your MakeDiffusion was not included in Mauro's MakeDiffusion.
Ofra (view profile)
Very useful.
Is there a Matlab 2008b version ?
Axel zheng (view profile)
Cool!!!
David S (view profile)
Just wanted to add that the error I mentioned is caused by the flipping of the sign of the surface normal of just one face by rearranging the order of the vertices in the face matrix
David S (view profile)
Very useful, except for a bug in the compute_curvature function:
compute_curvature will generate an error on line 75 ("dp = sum( normal(:,E(:,1)) .* normal(:,E(:,2)), 1 );") for SOME surfaces. The error stems from E containing indices that are out of range which is caused by line 48 ("A = sparse(double(i),double(j),s,n,n);") where A's values eventually entirely make up the E matrix. The problem occurs when the i and j vectors create the same ordered pair twice in which case the sparse function adds the two s vector elements together for that matrix location resulting in a value that is too large to be used as an index on line 75. For example, if i = [1 1] and j = [2 2] and s = [3 4] then A(1,2) will equal 3 + 4 = 7.
The i and j vectors are created here:
i = [face(1,:) face(2,:) face(3,:)];
j = [face(2,:) face(3,:) face(1,:)];
The fact that your code seems to depend on the order of the vertices in the faces matrix worries me because the curvature should be the same regardless of the order, obviously. To be fair, I don't completely understand how your code works so perhaps the way it is written it works out to not matter except that it does certainly matter when it results in an index out of bounds error as previously described.
very useful
very useful
Very helpful toolbox!!
good
Very useful
Perfect