As-Conformal-As-Possible Surface Registration

106 篇文章 1 订阅

{3.2. Conformal stiffness energy

LConformal(T)=T11T222+T22T332+T33T112+T12+T212+T23+T322+T31+T132

where

T=T11T21T31T12T22T32T13T23T33

Minimizing this objective function is equivalent to the following equation set:
T11T22=0
T22T33=0
T33T11=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:

101000000100000001000100110000000010000001000010011000000000000000000000T11T21T31T12T22T32T13T23T33t1t2t3=000000

}

{consist energe
i,jN(i) , Ti(v0jv0i)+v0i+ti(v0j+tj)=0
writing in matrix form leads to:

Ti(11)Ti(21)Ti(31)Ti(12)Ti(22)Ti(32)Ti(13)Ti(23)Ti(33)v0j(1)v0i(1)v0j(2)v0i(2)v0j(3)v0i(3)+v0i(1)v0i(2)v0i(3)+ti(1)ti(2)ti(3)v0j(1)+tj(1)v0j(2)+tj(2)v0j(3)+tj(3)=Ti(11)(v0j(1)v0i(1))+Ti(12)(v0j(2)v0i(2))+Ti(13)(v0j(3)v0i(3))+v0i(1)+ti(1)v0j(1)tj(1)Ti(21)(v0j(1)v0i(1))+Ti(22)(v0j(2)v0i(2))+Ti(23)(v0j(3)v0i(3))+v0i(1)+ti(2)v0j(2)tj(2)Ti(31)(v0j(1)v0i(1))+Ti(32)(v0j(2)v0i(2))+Ti(33)(v0j(3)v0i(3))+v0i(1)+ti(3)v0j(3)tj(3)=000

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

v0j(1)v0i(1)000v0j(1)v0i(1)000v0j(1)v0i(1)v0j(2)v0i(2)000v0j(2)v0i(2)000v0j(2)v0i(2)v0j(3)v0i(3)000v0j(3)v0i(3)000v0j(3)v0i(3)100010001100010001Ti(11)Ti(21)Ti(31)Ti(12)Ti(22)Ti(32)Ti(13)Ti(23)Ti(33)ti(1)ti(2)ti(3)tj(1)tj(2)tj(3)=v0j(1)v0i(1)v0j(2)v0i(2)v0j(3)v0i(3)

这里写图片描述

}

{smooth energe
i,jN(i) , Ti(v0jv0i)+Tj(v0iv0j)=0
writing in matrix form leads to:

Ti(11)Ti(21)Ti(31)Ti(12)Ti(22)Ti(32)Ti(13)Ti(23)Ti(33)v0j(1)v0i(1)v0j(2)v0i(2)v0j(3)v0i(3)+Tj(11)Tj(21)Tj(31)Tj(12)Tj(22)Tj(32)Tj(13)Tj(23)Tj(33)v0i(1)v0j(1)v0i(2)v0j(2)v0i(3)v0j(3)=Ti(11)(v0j(1)v0i(1))+Ti(12)(v0j(2)v0i(2))+Ti(13)(v0j(3)v0i(3))+Tj(11)(v0i(1)v0j(1))+Tj(12)(v0i(2)v0j(2))+Tj(13)(v0i(3)v0j(3))Ti(21)(v0j(1)v0i(1))+Ti(22)(v0j(2)v0i(2))+Ti(23)(v0j(3)v0i(3))+Tj(21)(v0i(1)v0j(1))+Tj(22)(v0i(2)v0j(2))+Tj(23)(v0i(3)v0j(3))Ti(31)(v0j(1)v0i(1))+Ti(32)(v0j(2)v0i(2))+Ti(33)(v0j(3)v0i(3))+Tj(31)(v0i(1)v0j(1))+Tj(32)(v0i(2)v0j(2))+Tj(33)(v0i(3)v0j(3))=000

v0j(1)v0i(1)000v0j(1)v0i(1)000v0j(1)v0i(1)v0j(2)v0i(2)000v0j(2)v0i(2)000v0j(2)v0i(2)v0j(3)v0i(3)000v0j(3)v0i(3)000v0j(3)v0i(3)v0i(1)v0j(1)000v0i(1)v0j(1)000v0i(1)v0j(1)v0i(2)v0j(2)000v0i(2)v0j(2)000v0i(2)v0j(2)v0i(3)v0j(3)000v0i(3)v0j(3)000v0i(3)v0j(3)Ti(11)Ti(21)Ti(31)Ti(12)Ti(22)Ti(32)Ti(13)Ti(23)Ti(33)Tj(11)Tj(21)Tj(31)Tj(12)Tj(22)Tj(32)Tj(13)Tj(23)Tj(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 positions
F=X0.tri'; % list of triangles indices
N=X0.vertex_mean_normal';
n = size(V,1); %number of vertices


if( 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 energe
count = 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
    end
end
t = 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 energe
count = 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
    end
end
t = toc;
fprintf('construction for smooth energe completes: %gs\n',t);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%conformal%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(conformal)
tic
%construct matrix for linear conformal energe
E_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];    

end
t = toc;
fprintf('construction for linear conformal energe completes: %gs\n',t);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%feature%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(feature)
tic
%construct matrix for feature point constraints


feature_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,:);
end
t = toc;
fprintf('construction for feature point constraints completes: %gs\n',t);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%closest%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(closest)
tic
%construct matrix for closest point constraints
closest_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;
end
t = toc;
fprintf('construction for cloest point constraints completes: %gs\n',t);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%total%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
%construct the total system
A=[];
b=[];
if(consist)
    A = [A; E_consist_A];
    b = [b; E_consist_b];
end
if(smooth)
    A = [A; E_smooth_A];
    b = [b; E_smooth_b];
end
if(conformal)
    A = [A; E_lconformal_A];
    b = [b; E_lconformal_b];    
end
if(feature)
    A = [A; E_feature_A];
    b = [b; E_feature_b];    
end
if(closest)
    A = [A; E_closest_A];
    b = [b; E_closest_b];    
end
t = toc;
fprintf('construction for total system completes: %gs\n',t);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%solving%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
%solve the whole system
%T=linsolve(A'*A,A'*b);%7.13842s

T = (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.724697s


t = 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 on
trimesh(F,V(:,1),V(:,2),V(:,3),'EdgeColor','b'); 
%trimesh(F_target,V_target(:,1),V_target(:,2),V_target(:,3),'EdgeColor','r'); 
axis equal
view(3)
disp('done!')
>> acap
construction for consist energe completes: 0.0397842s
construction for smooth energe completes: 0.0420681s
construction for linear conformal energe completes: 0.0217301s
construction for feature point constraints completes: 0.0140695s
construction for cloest point constraints completes: 0.0417374s
construction for total system completes: 1.08451s
solving completes: 4.71639s
done!
>> 

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

matlab 源代码

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值