function [tri,coord,nbr,normal]=MNIgetmesh(inputmesh)
% [tri,coord,nbr,normal]=MNIgetmesh(inputmesh)
% tri   = list of triangle element
% coord = coordinatates of vertices 
% nbr = neighboring vertex list
% normal = normal vectors
% inputmesh must be ASCII MNI format.
% there are total V + F rows.
% V = number of vertices
% F = number of faces
% it is assumed that the input mesh is topologically equvalent to a sphere, i.e.
% 2V - F = 4: Euler Characteristic equation.
% running time = 32.5770. see if you can speed up the algorithm.
%
% (C) 2004 Feb 1. Moo K. Chung, Sijie Tang
% Version 2.0.


fid=fopen(inputmesh);
frewind(fid);
fscanf(fid,'%c',1)
fscanf(fid,'%f',5)
n_points=fscanf(fid,'%i',1)
coord=fscanf(fid,'%f',[3,n_points]);
normal=fscanf(fid,'%f',[3,n_points]);
n_tri=fscanf(fid,'%i',1)
fscanf(fid,'%i',5+n_tri);
tri=fscanf(fid,'%i',[3,n_tri])'+1;
fclose(fid);


%---------------------------
% compute the maximum degree
degree=zeros(n_points,1);
for j=1:n_tri
degree(tri(j,:))=degree(tri(j,:))+1;
end
max_degree=max(degree);

%-------------------------
% compute the connectivity 
% can you come up with a simpler and faster algorithm.

nbr=zeros(n_points,max_degree);
for i_tri=1:n_tri
    for j=1:3
        cur_point = tri(i_tri,j);
        for k=1:3
            if (j ~= k)
                nbr_point= tri(i_tri,k);
                if find(nbr(cur_point,:)==nbr_point)
                    ;
                else
                    n_nbr = min(find(nbr(cur_point,:) == 0));
                    nbr(cur_point,n_nbr) = nbr_point;
                end;
            end;
        end;
    end;
end;

return;

