机器视觉 增强现实分享 http://blog.sciencenet.cn/u/wanglin193 把算法当成业余爱好......

博文

三角剖分网格的邻接矩阵和Laplacian矩阵实验

已有 360 次阅读 2018-9-8 19:55 |个人分类:3D|系统分类:科研笔记| Laplacian矩阵, Tutte mapping

平面上三角网格每个结点(vertice)的位置可以用1-ring邻接结点的位置加权得到。 未知结点Vp可以用权重矩阵W和全部结点V的乘积来表达: Vp=W*V,$$ V_p = 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的思路。



http://blog.sciencenet.cn/blog-465130-1133653.html

上一篇:SLAM系统的滤波和优化方法笔记
下一篇:固定边界的三维人脸网格参数化

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2018-9-24 08:28

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部