本文共 4248 字，大约阅读时间需要 14 分钟。

本质上是shape matching的as similar as possible版本

它的主要贡献在于, 找到一种方法使displacement field尽量平滑, 具体的做法就是公式7, 使每个点的偏移量, 尽量接近它周围邻居偏移量的平均值.

注意$t_k$公式在paper中是错误的

因为$r_k = s_kR_kx^0_k + t_k = s_kR_kx^0_k + c_k(t) - s_kR_kC^0_k = s_kR_k(x^0_k-C^0_k) + c_k(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);endfunction 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));endfunction 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));endfunction 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});endfunction r = compute_r(x, v0, s, R, t)    r = s(x,:) * R(:,:,x) * v0(x,:)' + t(x,:)';end

转载地址：http://bvxdi.baihongyu.com/