本质上是shape matching的as similar as possible版本
它的主要贡献在于, 找到一种方法使displacement field尽量平滑, 具体的做法就是公式7, 使每个点的偏移量, 尽量接近它周围邻居偏移量的平均值.
注意
tk
公式在paper中是错误的
应该为
tk=ck(t)−skRkC0k
因为 rk=skRkx0k+tk=skRkx0k+ck(t)−skRkC0k=skRk(x0k−C0k)+ck(t) 这样才符合shape matching中的逻辑
function step3_fitting
debug = 1;
step2_find_features(1,0,0.1,'interp');
load('Target_vertex','Target_vertex');
load('Target_faces','Target_faces');
vertex = load('coarse_vertex','coarse_vertex');
vertex = vertex.coarse_vertex;
v0 = vertex;
v3 = vertex;
n = length(v0);
faces = load('coarse_faces','coarse_faces');
faces = faces.coarse_faces;
neighbors = compute_vertex_ring(faces);
% cal c0
c0 = arrayfun(@(x) mean(v0([x neighbors{x}],:)), 1:n, 'UniformOutput', 0);
c0 = reshape(cell2mat(c0),3,[]);
c0 = c0';
A = zeros(3,3,n);
R = zeros(3,3,n);
s = zeros(n,1);
t = zeros(n,1);
D = zeros(n,1);
alpha = 0.95;
threshold = 0.2;
while alpha >= 0.5
[IDX,~] = knnsearch(Target_vertex, v3, 'k',1, 'NSMethod', 'kdtree');
v3_old = zeros(size(v3));
while sum( normrow(v3 - v3_old), 1) > 1
c3 = arrayfun(@(x) mean(v3([x neighbors{x}],:)), 1:n, 'UniformOutput', 0);
c3 = reshape(cell2mat(c3),3,[]);
c3 = c3';
f = @(x)compute_A(x, v3, c3, v0, c0, neighbors);
A = reshape(cell2mat(arrayfun(f,1:n, 'UniformOutput', 0)), [3 3 n]);
for ii = 1:n
% svd
[su,ss,sv]=svd(A(:,:,ii));
Ri = su*sv';
% if reflection then flip last column
if( det(Ri) < 0 )
su(:,end) = -su(:,end);
Ri = su*sv';
end
R(:,:,ii) = Ri;
end
f2 = @(x)compute_s(x, v3, c3, v0, c0, neighbors);
s = cell2mat(arrayfun(f2,1:n, 'UniformOutput', 0))';
t = arrayfun(@(x) c3(x,:)' - s(x,:) * R(:,:,x) * c0(x,:)' , 1:n, 'UniformOutput', 0);
t = cell2mat(t)';
f3 = @(x)compute_D(x, v3, v0, neighbors);
D = cell2mat(arrayfun(f3,1:n, 'UniformOutput', 0))';
v3_old = v3;
options.delete_patch = 'true';
[options] = my_plot_mesh(v3, faces, options);
IDX2 = find((D <= threshold));
if debug
num = size(IDX2, 1);
colorvector = (1:num)'/ num;
hold on
if exist('h11')~=0
delete(h11);
end
h11 = scatter3(v3(IDX2,1),v3(IDX2,2),v3(IDX2,3),20,colorvector);
if exist('h22')~=0
delete(h22);
end
h22 = scatter3(Target_vertex(IDX(IDX2),1),Target_vertex(IDX(IDX2),2),Target_vertex(IDX(IDX2),3),20,colorvector,'filled');
colormap jet(256);
%delete(h11);
%delete(h22);
end
f4 = @(x)compute_r(x, v0, s, R, t);
v3 = alpha * cell2mat(arrayfun(f4,1:n, 'UniformOutput', 0))' + bsxfun(@times, (1 - alpha) .* (D <= threshold), Target_vertex(IDX,:)) + bsxfun(@times, (1 - alpha) .* (D > threshold), cell2mat(arrayfun(f4,1:n, 'UniformOutput', 0))');
%v3 = alpha * cell2mat(arrayfun(f4,1:n, 'UniformOutput', 0))' + bsxfun(@times, (1 - alpha) , Target_vertex(IDX,:)) ;
sum( normrow(v3 - v3_old), 1)
%v3 = cell2mat(arrayfun(f4,1:n, 'UniformOutput', 0))';
getframe;
end
alpha = alpha - 0.05
end
write_obj('v3_mesh.obj', v3, faces);
end
function A = compute_A(x, v3, c3, v0, c0, neighbors)
A = bsxfun(@minus, v3([x neighbors{x}],:), c3(x))' * bsxfun(@minus, v0([x neighbors{x}],:), c0(x));
end
function s = compute_s(x, v3, c3, v0, c0, neighbors)
s = realsqrt(sum(normrow(bsxfun(@minus, v3([x neighbors{x}],:), c3(x))).^2,1) / sum(normrow(bsxfun(@minus, v0([x neighbors{x}],:), c0(x))).^2,1));
end
function D = compute_D(x, v3, v0, neighbors)
D = sum(abs(normrow(bsxfun(@minus, v3(neighbors{x},:), v3(x)))-normrow(bsxfun(@minus, v0(neighbors{x},:), v0(x)))) ./ normrow(bsxfun(@minus, v0(neighbors{x},:), v0(x))),1) / length(neighbors{x});
end
function r = compute_r(x, v0, s, R, t)
r = s(x,:) * R(:,:,x) * v0(x,:)' + t(x,:)';
end