本文共 4731 字，大约阅读时间需要 15 分钟。
As to the introduction of laplacian matrix, please refer to http://blog.csdn.net/seamanj/article/details/53026157
http://www.numerical-tours.com/matlab/
C:\TestMatlab\numerical-tours-master\matlab\toolbox_graph
compute_mesh_laplacian.m
function L = compute_mesh_laplacian(vertex,face,type,options)% compute_mesh_laplacian - compute a laplacian matrix%% L = compute_mesh_laplacian(vertex,face,type,options);%% If options.symmetrize=1 and options.normalize=0 then % L = D-W% If options.symmetrize=1 and options.normalize=1 then % L = eye(n)-D^{-1/2}*W*D^{-1/2}% If options.symmetrize=0 and options.normalize=1 then % L = eye(n)-D^{-1}*W.% where D=diag(sum(W,2)) and W is the unormalized weight matrix % (see compute_mesh_weight).%% type can 'combinatorial', 'distance', 'conformal'.%% See also compute_mesh_weight.%% Copyright (c) 2007 Gabriel Peyreoptions.null = 0;if isfield(options, 'normalize') normalize = options.normalize;else normalize = 1;endif isfield(options, 'symmetrize') symmetrize = options.symmetrize;else symmetrize = 1;endoptions.normalize = 0;W = compute_mesh_weight(vertex,face,type,options);n = size(W,1);if symmetrize==1 && normalize==0 L = diag(sum(W,2)) - W;elseif symmetrize==1 && normalize==1 L = speye(n) - diag(sum(W,2).^(-1/2)) * W * diag(sum(W,2).^(-1/2));elseif symmetrize==0 && normalize==1 L = speye(n) - diag(sum(W,2).^(-1)) * W;else error('Does not work with symmetrize=0 and normalize=0'); end compute_mesh_weight.m function W = compute_mesh_weight(vertex,face,type,options)% compute_mesh_weight - compute a weight matrix%% W = compute_mesh_weight(vertex,face,type,options);%% W is sparse weight matrix and W(i,j)=0 is vertex i and vertex j are not% connected in the mesh.%% type is either % 'combinatorial': W(i,j)=1 is vertex i is conntected to vertex j.% 'distance': W(i,j) = 1/d_ij^2 where d_ij is distance between vertex% i and j.% 'conformal': W(i,j) = cot(alpha_ij)+cot(beta_ij) where alpha_ij and% beta_ij are the adjacent angle to edge (i,j)%% If options.normalize=1, the the rows of W are normalize to sum to 1.%% Copyright (c) 2007 Gabriel Peyreoptions.null = 0;[vertex,face] = check_face_vertex(vertex,face);nface = size(face,1);n = max(max(face));verb = getoptions(options, 'verb', n>5000);if nargin<3 type = 'conformal';endswitch lower(type) case 'combinatorial' W = triangulation2adjacency(face); case 'distance' W = my_euclidean_distance(triangulation2adjacency(face),vertex); W(W>0) = 1./W(W>0); W = (W+W')/2; case 'conformal' % conformal laplacian W = sparse(n,n); for i=1:3 i1 = mod(i-1,3)+1; i2 = mod(i ,3)+1; i3 = mod(i+1,3)+1; pp = vertex(:,face(i2,:)) - vertex(:,face(i1,:)); qq = vertex(:,face(i3,:)) - vertex(:,face(i1,:)); % normalize the vectors pp = pp ./ repmat( sqrt(sum(pp.^2,1)), [3 1] ); qq = qq ./ repmat( sqrt(sum(qq.^2,1)), [3 1] ); % compute angles ang = acos(sum(pp.*qq,1)); W = W + sparse(face(i2,:),face(i3,:),cot(ang),n,n); W = W + sparse(face(i3,:),face(i2,:),cot(ang),n,n); end if 0 %% OLD CODE W = sparse(n,n); ring = compute_vertex_face_ring(face); for i = 1:n if verb progressbar(i,n); end for b = ring{i} % b is a face adjacent to a bf = face(:,b); % compute complementary vertices if bf(1)==i v = bf(2:3); elseif bf(2)==i v = bf([1 3]); elseif bf(3)==i v = bf(1:2); else error('Problem in face ring.'); end j = v(1); k = v(2); vi = vertex(:,i); vj = vertex(:,j); vk = vertex(:,k); % angles alpha = myangle(vk-vi,vk-vj); beta = myangle(vj-vi,vj-vk); % add weight W(i,j) = W(i,j) + cot( alpha ); W(i,k) = W(i,k) + cot( beta ); end end end otherwise error('Unknown type.')endif isfield(options, 'normalize') && options.normalize==1 W = diag(sum(W,2).^(-1)) * W;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function beta = myangle(u,v);du = sqrt( sum(u.^2) );dv = sqrt( sum(v.^2) );du = max(du,eps); dv = max(dv,eps);beta = acos( sum(u.*v) / (du*dv) );%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%function W = my_euclidean_distance(A,vertex)if size(vertex,1)