|||
平面上三角网格每个结点(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
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
边界固定,内部结点一次迭代就能确定位置。
矩阵每行都表示为当前结点和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:
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的思路。
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-8-17 06:20
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社