Deformable 3D shape registration based on local similarity transforms

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

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

注意 tk 公式在paper中是错误的
应该为 tk=ck(t)skRkC0k

因为 rk=skRkx0k+tk=skRkx0k+ck(t)skRkC0k=skRk(x0kC0k)+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
形状匹配 医学图像 Skuller A Volumetric Shape Registration Algorithm for Modeling Skull Deformities Medical Image Analysis Yusuf Sahilliog˘lu Ladislav Kavan We present an algorithm for volumetric registration of 3D solid shapes. In comparison to previous work on image based registration, our technique achieves higher efficiency by leveraging a template tetrahedral mesh. In contrast to point- and surface-based registration techniques, our method better captures volumetric nature of the data, such as bone thickness. We apply our algorithm to study pathological skull deformities caused by a particular condition, i.e., craniosynostosis. The input to our system is a pair of volumetric 3D shapes: a tetrahedral mesh and a voxelized object represented by a set of voxel cells segmented from computed tomography (CT) scans. Our general framework first performs a global registration and then launches a novel elastic registration process that uses as much volumetric information as possible while deforming the generic template tetrahedral mesh of a healthy human skull towards the underlying geometry of the voxel cells. Both data are high-resolution and differ by large non-rigid deformations. Our fully-automatic solution is fast and accurate, as compared with the state of the arts from the reconstruction and medical image registration fields. We use the resulting registration to match the ground-truth surfaces extracted from the medical data as well as to quantify the severity of the anatomical deformity
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值