平面上三角网格每个结点(vertice)的位置可以用1-ring邻接结点的位置加权得到。 未知结点Vp可以用权重矩阵W和全部结点V的乘积来表达: Vp=W*V, 其中稀疏矩阵W表达相邻结点间的邻接关系的权重矩阵（参考图论的邻接矩阵Adjacency Matrix）, 右边的V既包含未知结点Vp，也包含已知结点Vh。迭代过程Vp收敛到固定点，这个过程也可以用Laplacian线性方程求解：L*V=0，只要把方程左边已知固定结点对应的列挪到右手边即可求解。

生成一个简单的三角剖分网格

close allclear all[x,y] = meshgrid(-1:0.25:1,-1:0.25:1);
[nr,nc] = size(x);
n = nr*nc;
x = x(:); y = y(:);
vp =[x,y];
% Add circle points to tail
r = 2;
%num of boundary points
nb = (nc+nc+1)*2;
th = (1:nb)'*(2*pi)/nb;
vh = r * [cos(th),sin(th)];
vertex = [vp;vh];

Delaunay方法生成三角网格 

faces = delaunay(vertex(:,1),vertex(:,2));
idx_inner = (1:n);
idx_cont = (n+1:n+nb);
figure,triplot(faces,vertex(:,1),vertex(:,2));
hold on,plot(vertex(idx_cont,1),vertex(idx_cont,2),'r.'),axis equal

邻接权重稀疏矩阵W

Collect list of half edges

edges = [faces(:,1),faces(:,2); faces(:,2),faces(:,3); faces(:,3),faces(:,1)];
% weight 
matrixw = zeros(n+nb,n+nb);
for i=1:size(edges,1)
    w(edges(i,1),edges(i,2)) = 1;    
    % turn to indirect edges
    w(edges(i,2),edges(i,1)) = 1;
end
d = sum(w,2) ;
id = diag(1./d);
w = id * w;

固定边界点进行加权平均

显示1,3,5,20步结果

eye_n = eye(n+nb,n+nb);
% fix some vertex
w(idx_cont,:) = eye_n(idx_cont,:);
V = vertex;
for iter = 1:20
   V = w * V;    
   if(iter==1 || iter==3 || iter==5 || iter==20)
       figure, triplot(faces,V(:,1),V(:,2));
       hold on,plot(V(idx_cont,1),V(idx_cont,2),'r.'),axis equal
   end
end

第1次迭代：

第3次迭代：

第5次迭代：

第20次迭代：

矩阵连乘

矩阵预先连乘也能达到同样效果 迭代过程中W每行的和是1，保持不变，但是越来越像右侧固定边界对应的部分集中. 如果只把矩阵右侧nb列矩阵按行求和，那么会越来越接近1

W = w*w;
figure(100),hold on,ylim([0,1.1])
for iter = 1:8
   W = W^2;
   s = sum(W(:,n+1:end),2);
   plot(s)
   title('Weight sum of fix point part ')
end

迭代过程“能量”向右侧的nb列集中：

权重矩阵的图像

从矩阵W可以看出，每行除了固定边界对应的部分其他都接近0。

figure,imshow(W,[]);

表示未知结点位置只和边界点有关,这样矩阵只保留右侧和边界结点对应的部分即可 

V = vertex;
W(idx_inner,idx_inner)=0;
V=W*V;
figure,triplot(faces,V(:,1),V(:,2));
hold on,plot(V(idx_cont,1),V(idx_cont,2),'r.'),axis equal

边界固定，内部结点一次迭代就能确定位置。

Laplacian 矩阵 L

矩阵每行都表示为当前结点和1-ring邻接结点的差的求和形式 for each row, sum(L(r,:))=0

L = eye(n+nb,n+nb) - w;
% remove known vertex lines from L
L(idx_cont,:)=[];
% move column of boundary vertex in LV=0 to right
Rhs = -L(:,idx_cont) * vh;
% Unknow vertex on left
Lhs = L(:,idx_inner);
Vp = Lhs\Rhs; %inv(Lhs)*Rhs
% copy back to vertex
vertex(1:n,:) = Vp;

用Laplacian矩阵直接求解的结果和上面W矩阵迭代收敛的结果一致

figure,triplot(faces,vertex(:,1),vertex(:,2));
hold on,plot(V(idx_cont,1),V(idx_cont,2),'r.'),axis equal

Matlab Code：

demo1.m

 Olga Sorkine 早年的的相关文章：

[1]. Least-squares Meshes, Olga Sorkine and Daniel Cohen-Or

LV=0方程组中，已知固定边界结点对应的行并不从方程中移走(硬约束hard constraints)，而是以w*I_n*V=w*Vh形式添加到方程组后面，就是所谓的Least-sqeare soft constraints。这样网格编辑的时候，已知结点并不总能精确地固定在给定的位置。

[2].Laplacian Mesh Processing, Olga Sorkine

介绍了W矩阵中wij的计算方法，不象上面是用wij=1，权重应该和邻域三角形的角度面积等相关。和Cotangent余切函数相关的权重在曲面保角变换(或称Conformal共形变换)变形中非常有用。

[3]. On Linear Variational Surface Deformation Methods, Mario Botsch and Olga Sorkine

很好的综述。Laplace-Beltrami算子是laplacian在离散化的三角曲面上的近似计算，尽管它本身不是旋转无关的(只是唯一无关)，但是配合局部的刚体变换就能比较完美地解决曲面变形问题。

[4].As-Rigid-As-Possible Surface Modeling, Olga Sorkine and Marc Alexa

根据已知对应点位置，根据LV=0，计算未知点的初始位置；再根据邻域（比如1-ring）计算局部刚体变换，相似变换或者仿射变换，调整局部frame的方向。这两个步骤的迭代是一种可行的解决surface deformation的思路。