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 @ NCCA

n = 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 constraints
w=1;

% building the least-squares system matrix
A_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-primes
xyz_col = A_prime\rhs;
U = [xyz_col(1:n) xyz_col((n+1):(2*n)) xyz_col((2*n+1):(3*n))];