As-Conformal-As-Possible Surface Registration

来源：互联网发布：python获取字符串长度编辑：程序博客网时间：2024/05/16 04:46

{3.2. Conformal stiffness energy

LConformal(T)=∥T11−T22∥2+∥T22−T33∥2+∥T33−T11∥2+∥T12+T21∥2+∥T23+T32∥2+∥T31+T13∥2

where

T = ⎡ ⎣ ⎢ T 11 T 21 T 31 T 12 T 22 T 32 T 13 T 23 T 33 ⎤ ⎦ ⎥

Minimizing this objective function is equivalent to the following equation set:
∥T11−T22∥=0
∥T22−T33∥=0
∥T33−T11∥=0
∥T12+T21∥=0
∥T23+T32∥=0
∥T31+T13∥=0
[see the fourth one in http://blog.csdn.net/seamanj/article/details/51803822 , which can be applied in following derivations]

which can be written in matrix form:

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 10 - 1 000 000100000001000100 - 1 10000 000010000001000010 0 - 1 1000 000000000000000000 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ T 11 T 21 T 31 T 12 T 22 T 32 T 13 T 23 T 33 t 1 t 2 t 3 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 000000 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

}

{consist energe
∀i,∀j∈N(i), ∥Ti(v0j−v0i)+v0i+ti−(v0j+tj)∥=0
writing in matrix form leads to:

⎡ ⎣ ⎢ T i (11) T i (21) T i (31) T i (12) T i (22) T i (32) T i (13) T i (23) T i (33) ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎢ v 0 j (1) - v 0 i (1) v 0 j (2) - v 0 i (2) v 0 j (3) - v 0 i (3) ⎤ ⎦ ⎥ ⎥ + ⎡ ⎣ ⎢ ⎢ v 0 i (1) v 0 i (2) v 0 i (3) ⎤ ⎦ ⎥ ⎥ + ⎡ ⎣ ⎢ t i (1) t i (2) t i (3) ⎤ ⎦ ⎥ - ⎡ ⎣ ⎢ ⎢ v 0 j (1) + t j (1) v 0 j (2) + t j (2) v 0 j (3) + t j (3) ⎤ ⎦ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ T i (11) (v 0 j (1) - v 0 i (1)) + T i (12) (v 0 j (2) - v 0 i (2)) + T i (13) (v 0 j (3) - v 0 i (3)) + v 0 i (1) + t i (1) - v 0 j (1) - t j (1) T i (21) (v 0 j (1) - v 0 i (1)) + T i (22) (v 0 j (2) - v 0 i (2)) + T i (23) (v 0 j (3) - v 0 i (3)) + v 0 i (1) + t i (2) - v 0 j (2) - t j (2) T i (31) (v 0 j (1) - v 0 i (1)) + T i (32) (v 0 j (2) - v 0 i (2)) + T i (33) (v 0 j (3) - v 0 i (3)) + v 0 i (1) + t i (3) - v 0 j (3) - t j (3) ⎤ ⎦ ⎥ ⎥ = ⎡ ⎣ ⎢ 000 ⎤ ⎦ ⎥

Collecting all the unknown variables into a column vector yields us:

⎡ ⎣ ⎢ ⎢ v 0 j (1) - v 0 i (1) 00 0 v 0 j (1) - v 0 i (1) 0 00 v 0 j (1) - v 0 i (1) v 0 j (2) - v 0 i (2) 00 0 v 0 j (2) - v 0 i (2) 0 00 v 0 j (2) - v 0 i (2) v 0 j (3) - v 0 i (3) 00 0 v 0 j (3) - v 0 i (3) 0 00 v 0 j (3) - v 0 i (3) 100010001 ⎤ ⎦ ⎥ ⎥ \dots ⎡ ⎣ ⎢ - 1 00 0 - 1 0 00 - 1 ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ T i (11) T i (21) T i (31) T i (12) T i (22) T i (32) T i (13) T i (23) T i (33) t i (1) t i (2) t i (3) ⋮ t j (1) t j (2) t j (3) ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ v 0 j (1) - v 0 i (1) v 0 j (2) - v 0 i (2) v 0 j (3) - v 0 i (3) ⎤ ⎦ ⎥ ⎥

}

{smooth energe
∀i,∀j∈N(i), ∥Ti(v0j−v0i)+Tj(v0i−v0j)∥=0
writing in matrix form leads to:

⎡ ⎣ ⎢ T i (11) T i (21) T i (31) T i (12) T i (22) T i (32) T i (13) T i (23) T i (33) ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎢ v 0 j (1) - v 0 i (1) v 0 j (2) - v 0 i (2) v 0 j (3) - v 0 i (3) ⎤ ⎦ ⎥ ⎥ + ⎡ ⎣ ⎢ T j (11) T j (21) T j (31) T j (12) T j (22) T j (32) T j (13) T j (23) T j (33) ⎤ ⎦ ⎥ ⎡ ⎣ ⎢ ⎢ v 0 i (1) - v 0 j (1) v 0 i (2) - v 0 j (2) v 0 i (3) - v 0 j (3) ⎤ ⎦ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ T i (11) (v 0 j (1) - v 0 i (1)) + T i (12) (v 0 j (2) - v 0 i (2)) + T i (13) (v 0 j (3) - v 0 i (3)) + T j (11) (v 0 i (1) - v 0 j (1)) + T j (12) (v 0 i (2) - v 0 j (2)) + T j (13) (v 0 i (3) - v 0 j (3)) T i (21) (v 0 j (1) - v 0 i (1)) + T i (22) (v 0 j (2) - v 0 i (2)) + T i (23) (v 0 j (3) - v 0 i (3)) + T j (21) (v 0 i (1) - v 0 j (1)) + T j (22) (v 0 i (2) - v 0 j (2)) + T j (23) (v 0 i (3) - v 0 j (3)) T i (31) (v 0 j (1) - v 0 i (1)) + T i (32) (v 0 j (2) - v 0 i (2)) + T i (33) (v 0 j (3) - v 0 i (3)) + T j (31) (v 0 i (1) - v 0 j (1)) + T j (32) (v 0 i (2) - v 0 j (2)) + T j (33) (v 0 i (3) - v 0 j (3)) ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ 000 ⎤ ⎦ ⎥

⎡ ⎣ ⎢ ⎢ v 0 j (1) - v 0 i (1) 00 0 v 0 j (1) - v 0 i (1) 0 00 v 0 j (1) - v 0 i (1) v 0 j (2) - v 0 i (2) 00 0 v 0 j (2) - v 0 i (2) 0 00 v 0 j (2) - v 0 i (2) v 0 j (3) - v 0 i (3) 00 0 v 0 j (3) - v 0 i (3) 0 00 v 0 j (3) - v 0 i (3) ⎤ ⎦ ⎥ ⎥ \dots ⎡ ⎣ ⎢ ⎢ v 0 i (1) - v 0 j (1) 00 0 v 0 i (1) - v 0 j (1) 0 00 v 0 i (1) - v 0 j (1) v 0 i (2) - v 0 j (2) 00 0 v 0 i (2) - v 0 j (2) 0 00 v 0 i (2) - v 0 j (2) v 0 i (3) - v 0 j (3) 00 0 v 0 i (3) - v 0 j (3) 0 00 v 0 i (3) - v 0 j (3) ⎤ ⎦ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ T i (11) T i (21) T i (31) T i (12) T i (22) T i (32) T i (13) T i (23) T i (33) ⋮ T j (11) T j (21) T j (31) T j (12) T j (22) T j (32) T j (13) T j (23) T j (33) ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ 000 ⎤ ⎦ ⎥

这里写图片描述

}

At last, I would like release the matlab source of this paper:

clear;consist = 1;smooth = 1;conformal = 1;feature = 1;closest = 1;X0=read_obj('sphere_before.obj');V=X0.xyz';  % list of vertex positionsF=X0.tri'; % list of triangles indicesN=X0.vertex_mean_normal';n = size(V,1); %number of verticesif( feature || closest )     X1=read_obj('sphere_ffd.obj');      V_target=X1.xyz';  % list of vertex positions     F_target=X1.tri';end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%finding neighbour%%%%%%%%%%%%%%%%%%%%%%%%%%%%neighbour_matrix = zeros(n,n);for i=1:n    neighbour = find(F(:,1)==i | F(:,2)==i | F(:,3)==i);    neighbour = F(neighbour,:);    neighbour = setdiff(neighbour(:),i);    neighbour_matrix(i,neighbour(:)) = ones(1,size(neighbour,2));end neighbour_count = sum(sum(neighbour_matrix,2),1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%consist%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(consist)E_consist_A=zeros(3*neighbour_count, 12*n);E_consist_b=zeros(3*neighbour_count, 1);tic%%construct matrix for consist energecount = 1;for i = 1:n    for j = 1:n        if( neighbour_matrix(i,j) && i~=j)             E_consist_A((count-1)*3+1 : count*3, (i-1)*12 + 1 : i*12) = ...                 [eye(3)*(V(j,1)-V(i,1)) eye(3)*(V(j,2)-V(i,2)) eye(3)*(V(j,3)-V(i,3)) eye(3)];             E_consist_A((count-1)*3+1 : count*3, (j-1)*12 + 10 : j*12) = ...                  -1 * eye(3);             E_consist_b((count-1)*3+1 : count*3,1) = V(j,:)' - V(i,:)';             count = count + 1;        end    endendt = toc;fprintf('construction for consist energe completes: %gs\n',t);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%smooth%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(smooth)E_smooth_A=zeros(3*neighbour_count, 12*n);E_smooth_b=zeros(3*neighbour_count, 1);    tic%construct matrix for smooth energecount = 1;for i = 1:n    for j = 1:n        if( neighbour_matrix(i,j) && i~=j)            % E_smooth_A_temp = sparse(3, 12*n);            E_smooth_A( (count-1)*3+1 : count*3, (i-1)*12 + 1 : (i-1)*12 + 9) = ...            [eye(3)*(V(j,1)-V(i,1)) eye(3)*(V(j,2)-V(i,2)) eye(3)*(V(j,3)-V(i,3))];            E_smooth_A( (count-1)*3+1 : count*3, (j-1)*12 + 1 : (j-1)*12 + 9) = ...            [eye(3)*(V(i,1)-V(j,1)) eye(3)*(V(i,2)-V(j,2)) eye(3)*(V(i,3)-V(j,3))];            count = count + 1;        end    endendt = toc;fprintf('construction for smooth energe completes: %gs\n',t);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%conformal%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(conformal)tic%construct matrix for linear conformal energeE_lconformal_A = zeros(6*n, 12*n);E_lconformal_b = zeros(6*n, 1);for i = 1:n             E_lconformal_A( (i-1)*6 + 1 : i*6, (i-1)*12 + 1 : (i-1)*12 + 9) = ...             [ 1 0 0   0 -1 0   0 0  0;               0 0 0   0  1 0   0 0 -1;              -1 0 0   0  0 0   0 0  1;               0 1 0   1  0 0   0 0  0;               0 0 0   0  0 1   0 1  0;               0 0 1   0  0 0   1 0  0];    endt = toc;fprintf('construction for linear conformal energe completes: %gs\n',t);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%feature%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(feature)tic%construct matrix for feature point constraintsfeature_matching_pairs = [1:200; 1:200]';E_feature_A = zeros(3*n, 12*n);E_feature_b = zeros(3*n, 1);for i = 1:size(feature_matching_pairs,1)    template_index = feature_matching_pairs(i,1);    target_index = feature_matching_pairs(i,2);    E_feature_A( (template_index-1)*3+1 : template_index*3, (template_index-1)*12+10 : template_index*12) = ...    eye(3);    E_feature_b((template_index-1)*3+1 : template_index*3) = V_target(target_index,:) - V(template_index,:);endt = toc;fprintf('construction for feature point constraints completes: %gs\n',t);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%closest%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(closest)tic%construct matrix for closest point constraintsclosest_matching_pairs = [1:200; 1:200]';E_closest_A = zeros(3*n, 12*n);E_closest_b = zeros(3*n, 1);for i = 1:size(closest_matching_pairs,1)    template_index = closest_matching_pairs(i,1);    target_index = closest_matching_pairs(i,2);    E_closest_A( (template_index-1)*3+1 : template_index*3, (template_index-1)*12+10 : template_index*12) = ...    eye(3);    VP = V_target(target_index,:) - V(template_index,:);    Normal = N(template_index,:);    E_closest_b((template_index-1)*3+1 : template_index*3) = dot(VP,Normal).*Normal;endt = toc;fprintf('construction for cloest point constraints completes: %gs\n',t);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%total%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%tic%construct the total systemA=[];b=[];if(consist)    A = [A; E_consist_A];    b = [b; E_consist_b];endif(smooth)    A = [A; E_smooth_A];    b = [b; E_smooth_b];endif(conformal)    A = [A; E_lconformal_A];    b = [b; E_lconformal_b];    endif(feature)    A = [A; E_feature_A];    b = [b; E_feature_b];    endif(closest)    A = [A; E_closest_A];    b = [b; E_closest_b];    endt = toc;fprintf('construction for total system completes: %gs\n',t);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%solving%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%tic%solve the whole system%T=linsolve(A'*A,A'*b);%7.13842sT = (A'*A)\(A'*b);%6.63639s%T=pcg(A'*A,A'*b,1e-12,2000);%9.59388s %pcg converged at iteration 649 to a solution with relative residual 8.8e-13.%  A = sparse(A);%  b = sparse(b);%  [L,U,P,Q,R] = lu(A);%  T = Q*(U\(L\(P*(R\b))));%0.724697st = toc;fprintf('solving completes: %gs\n',t);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%update%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%for i = 1:n    V(i,:) = V(i,:) + T((i-1)*12+10 : i*12)';end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%hold ontrimesh(F,V(:,1),V(:,2),V(:,3),'EdgeColor','b'); %trimesh(F_target,V_target(:,1),V_target(:,2),V_target(:,3),'EdgeColor','r'); axis equalview(3)disp('done!')

>> acapconstruction for consist energe completes: 0.0397842sconstruction for smooth energe completes: 0.0420681sconstruction for linear conformal energe completes: 0.0217301sconstruction for feature point constraints completes: 0.0140695sconstruction for cloest point constraints completes: 0.0417374sconstruction for total system completes: 1.08451ssolving completes: 4.71639sdone!>>

the figure after running seems like:
这里写图片描述

matlab 源代码

0 0