本文共 2914 字,大约阅读时间需要 9 分钟。
先讲下各个公式
公式(2)是laplacian coordinates的定义
公式(5)第一部分是保证laplacian coordinate坐标的一致性, 后面是handle点的约束
公式(8)是T的构造
公式(12)是求T中的系数
2D MATLAB 版在它的project里面自带, 这里将其拓展成3D版本, 源代码如下:
function U = laplacian_surface_editing_3D(vertex,faces,BI,BC)%The file is an implementation of 'Laplacian surface editing' in 3D. The%original accompany code with this paper is in 2D. Based on this code, I%build the 3D version.%Inputs: vertex % vertex #vertex by 3 list of rest domain positions % faces #faces by 3 list of triangle indices into vertex % b #b list of indices of constraint (boundary) vertices % bc #b by 3 list of constraint positions for b %Output: % U #V by dim list of new positions% By seamanj @ NCCAn = length(vertex);options.symmetrize=0;options.normalize=1;L = compute_mesh_laplacian(vertex, faces, 'combinatorial', options );delta = L * vertex; %delta为拉普拉斯坐标L_prime = [ L zeros(n) zeros(n) % the x-part zeros(n) L zeros(n) zeros(n) zeros(n) L ]; neighbors = compute_vertex_ring(faces); for i = 1:n ring = [i neighbors{i}]; V = vertex(ring,:)'; V = [V ones(1,length(ring))];%这里为1, 乘上v'就变成公式(10)中的A了,这里这么写就是把v'合并了 %V第一行为x,第二行为y,第三行为z,第四行为1 C = zeros(length(ring) * 3, 7); % ... Fill C in for r=1:length(ring) C(r,:) = [V(1,r) 0 V(3,r) (-1)*V(2,r) V(4,r) 0 0]; %V第一行为x,第二行为y,第三行为z,第四行为1 C(length(ring)+r,:) = [V(2,r) (-1)*V(3,r) 0 V(1,r) 0 V(4,r) 0]; C(2*length(ring)+r,:) = [V(3,r) V(2,r) (-1)*V(1,r) 0 0 0 V(4,r)]; end; Cinv = pinv(C); s = Cinv(1,:); h1 = Cinv(2,:); h2 = Cinv(3,:); h3 = Cinv(4,:); delta_i = delta(i,:)'; delta_ix = delta_i(1); delta_iy = delta_i(2); delta_iz = delta_i(3); % T*delta gives us an array of coefficients % T*delta*V'等于公式(5)的T(V')*delta, 注意这里我们要求的是V' Tdelta = [delta_ix*s + delta_iy*(-1)*h3 + delta_iz*h2 delta_ix*h3 + delta_iy*s + delta_iz*(-1)*h1 delta_ix*(-1)*h2 + delta_iy*h1 + delta_iz*s]; % updating the weights in Lx_prime, Ly_prime, Lw_prime %L_prime里面原本有L,这里就是公式(5)中的T(V')*delta - L(V') L_prime(i,[ring (ring + n) (ring + 2*n)]) = ... L_prime(i,[ring (ring + n) (ring + 2*n)]) + (-1)*Tdelta(1,:); L_prime(i+n,[ring (ring + n) (ring + 2*n)]) = ... L_prime(i+n,[ring (ring + n) (ring + 2*n)]) + (-1)*Tdelta(2,:); L_prime(i+n*2,[ring (ring + n) (ring + 2*n)]) = ... L_prime(i+n*2,[ring (ring + n) (ring + 2*n)]) + (-1)*Tdelta(3,:); end % weight for the constraintsw=1;% building the least-squares system matrixA_prime = L_prime;rhs = zeros(3*n,1);for j=1:length(BI) A_prime = [A_prime w*((1:(3*n))==BI(j)) w*((1:(3*n))==(BI(j)+n)) w*((1:(3*n))==(BI(j)+2*n))]; rhs = [rhs w*BC(j,1) w*BC(j,2) w*BC(j,3)];end;% solving for v-primesxyz_col = A_prime\rhs;U = [xyz_col(1:n) xyz_col((n+1):(2*n)) xyz_col((2*n+1):(3*n))];
转载地址:http://yvxdi.baihongyu.com/