mesh_laplacian - Laplacian of irregular triangular mesh Usage: [lap,edge] = mesh_laplacian(vertex,face) Returns 'lap', the Laplacian (2nd spatial derivative) of an irregular triangular mesh, and 'edge', the linear distances between vertices of 'face'. 'lap' and 'edge' are square, [Nvertices,Nvertices] in size, sparse in nature. It is assumed that 'vertex' contains the (x,y,z) Cartesian coordinates of each vertex and that 'face' contains the triangulation of vertex with indices into 'vertex' that are numbered from 1:Nvertices. For information about triangulation, see 'help convhull' or 'help convhulln'. The neighbouring vertices of vertex 'i' is given by: k = find(edge(i,:)); The math of this routine is given by: Oostendorp, Oosterom & Huiskamp (1989), Interpolation on a triangulated 3D surface. Journal of Computational Physics, 80: 331-343. See also, eeg_interp_scalp_mesh
0001 function [lap,edge] = mesh_laplacian(vertex,face) 0002 0003 % mesh_laplacian - Laplacian of irregular triangular mesh 0004 % 0005 % Usage: [lap,edge] = mesh_laplacian(vertex,face) 0006 % 0007 % Returns 'lap', the Laplacian (2nd spatial derivative) of an 0008 % irregular triangular mesh, and 'edge', the linear distances 0009 % between vertices of 'face'. 'lap' and 'edge' are square, 0010 % [Nvertices,Nvertices] in size, sparse in nature. 0011 % 0012 % It is assumed that 'vertex' contains the (x,y,z) Cartesian 0013 % coordinates of each vertex and that 'face' contains the 0014 % triangulation of vertex with indices into 'vertex' that 0015 % are numbered from 1:Nvertices. For information about 0016 % triangulation, see 'help convhull' or 'help convhulln'. 0017 % 0018 % The neighbouring vertices of vertex 'i' is given by: 0019 % 0020 % k = find(edge(i,:)); 0021 % 0022 % The math of this routine is given by: 0023 % 0024 % Oostendorp, Oosterom & Huiskamp (1989), 0025 % Interpolation on a triangulated 3D surface. 0026 % Journal of Computational Physics, 80: 331-343. 0027 % 0028 % See also, eeg_interp_scalp_mesh 0029 % 0030 0031 % $Revision: 1.2 $ $Date: 2005/08/14 21:25:08 $ 0032 0033 % Licence: GNU GPL, no implied or express warranties 0034 % History: 04/2002, Darren.Weber_at_radiology.ucsf.edu 0035 % - initial version was inefficient and incorrect 0036 % at one point. 0037 % (c) 04/2002 Robert Oostenveld 0038 % - completely revised/reconstructed code (lapcal.m) 0039 % - agreed to release into eeg_toolbox under GNU GPL 0040 % 04/2002, Darren.Weber_at_radiology.ucsf.edu 0041 % - modified edge initialization from sparse to 0042 % full matrix and slightly improved speed of 0043 % calculation for edge norms 0044 % - added tic/toc timing report 0045 % 07/2002, Darren.Weber_at_radiology.ucsf.edu 0046 % - added check for duplicate vertices (lines 78-79) 0047 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0048 0049 nvertex = size(vertex,1); 0050 nface = size(face,1); 0051 0052 fprintf('MESH_LAPLACIAN: Calc Laplacian matrix for %5d vertices...',nvertex); 0053 tic 0054 0055 % the matrix 'edge' is the connectivity of all vertices 0056 edge = sparse(nvertex,nvertex); 0057 for i=1:nface, 0058 0059 % compute the length of all triangle edges (Diff is [3x3]) 0060 Diff = [vertex(face(i,[1 2 3]),:) - vertex(face(i,[2 3 1]),:)]; 0061 Norm = sqrt( sum(Diff.^2, 2) ); 0062 0063 edge(face(i,1),face(i,2)) = Norm(1); 0064 edge(face(i,2),face(i,3)) = Norm(2); 0065 edge(face(i,3),face(i,1)) = Norm(3); 0066 0067 % make sure that all edges are symmetric 0068 edge(face(i,2),face(i,1)) = Norm(1); 0069 edge(face(i,3),face(i,2)) = Norm(2); 0070 edge(face(i,1),face(i,3)) = Norm(3); 0071 end 0072 0073 % Using edge to identify nearest vertices, calculate 0074 % the Laplacian for an irregular mesh 0075 lap = sparse(nvertex,nvertex); 0076 for i=1:nvertex, 0077 0078 k = find(edge(i,:)); % the indices of the neighbours 0079 0080 % remove any duplicate neighbour vertices 0081 [vert, m, n] = unique(vertex(k,:),'rows'); 0082 k = k(m); 0083 0084 ni = length(k); % the number of neighbours 0085 0086 hi = mean(edge(i,k)); % the average distance to the neighbours 0087 invhi = mean(1./edge(i,k)); % the average inverse distance to the neighbours 0088 0089 lap(i,i) = -(4/hi) * invhi; % Laplacian of vertex itself 0090 0091 lap(i,k) = (4/(hi*ni)) * 1./edge(i,k); % Laplacian of direct neighbours 0092 0093 % Laplacian is zero for all indirect neighbours 0094 % See Oostendorp, Oosterom & Huiskamp (1989, pp. 334-335) 0095 end 0096 0097 t = toc; 0098 fprintf('done (%6.2f sec).\n',t); 0099 0100 return