Deformable 3D shape registration based on local similarity transforms

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

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

注意$t_k$公式在paper中是错误的
应该为 $t_k = c_k(t) - s_kR_kC^0_k$

因为$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);
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