Laplacian surface editing

先讲下各个公式

公式(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  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))];