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))];