本文共 9040 字,大约阅读时间需要 30 分钟。
{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
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 , which can be applied in following derivations]which can be written in matrix form:
}
{consist energe
∀i,∀j∈N(i) , ∥Ti(v0j−v0i)+v0i+ti−(v0j+tj)∥=0 writing in matrix form leads to:Collecting all the unknown variables into a column vector yields us:
}
{smooth energe
∀i,∀j∈N(i) , ∥Ti(v0j−v0i)+Tj(v0i−v0j)∥=0 writing in matrix form leads to:
}
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: