高效的计算Laplacian 矩阵

function L1 = CP_Laplacian( M )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
% number of vertices (size of matrix)
nverts = size(M.vertices,1);

% assume manifold and faces oriented uniformly
edges = [M.faces(:,[1,2]); M.faces(:,[2,3]); M.faces(:,[3,1])];

% extract vertices from matrix
v1Idxs = edges(:,1);
v2Idxs = edges(:,2);
v1s = M.vertices(v1Idxs,:);
v2s = M.vertices(v2Idxs,:);

% compute edge lengths and weights
eps = 1e-8;
edgel = sqrt( sum( (v1s-v2s).^2, 2 ) );
edgew = 1 ./ (eps + edgel );

% construct sparse matrix
nnzeros = size(edges,1) / 2;
L1 = sparse( v1Idxs, v2Idxs, edgew, nverts, nverts );

% compute diagonal by summing column-major (faster)
D = -full( sum( L1,1 ) )‘;

% insert D in 0-diagonal of L1
L1 = spdiags(D,0,L1);

end