基于微分同胚的Demons配准算法

基于微分同胚的Demons配准算法，算法原型是基于Diffeomorphic demons efficient non-parametric image registration

% Demons Registration

function demons

figure（1）; clf; colormap gray; set（gcf‘renderer‘‘painter‘）

%% Parameters
niter           = 250;
sigma_fluid     = 2.0; % regularize update      field
sigma_diffusion = 2.0; % regularize deformation field
sigma_i         = 1.0; % weight on similarity term
sigma_x         = 1.0; % weight on spatial uncertainties （maximal step）
nlevel          = 3;   % multiresolution
do_display      = 1;   % display iterations

%% Load fixed image
F = double（imread（‘data/lenag2.png‘））;          % fixed   image
%F = double（imread（‘statue-rio.png‘））;           % fixed   image
F = imresize（F0.5）;
F = 256*（F-min（F（:）））/range（F（:））;              % normalize intensities

%% Load moving image
M = double（imread（‘data/lenag1.png‘））;          % moving  image
%M = double（imread（‘statue-rio-deformed.png‘））;  % moving  image
M = imresize（M0.5）;
M = 256*（M-min（M（:）））/range（M（:））;              % normalize intensities

% Translate
%shift = 3; tmp = zeros（size（M））; tmp（（1+shift）:end:） = M（1:（end-shift）:）; M = tmp;

%% Create random moving image
%[Fsx0sy0]  = randomdeform（F505）;
%figure（5）; showvector（sx0sy043lim）; drawnow;

if nlevel == 1
    
    %% Register
    disp（[‘Register...‘]）;
    opt = struct（‘niter‘niter ‘sigma_fluid‘sigma_fluid ‘sigma_diffusion‘sigma_diffusion ‘sigma_i‘sigma_i ‘sigma_x‘sigma_x ‘do_display‘do_display ‘do_plotenergy‘1）;
    [Mpsxsyvxvy] = register（FMopt）;

else
    
    %% Multiresolution
    vx = zeros（size（M））; % deformation field
    vy = zeros（size（M））;
    for k=nlevel:-1:1
        disp（[‘Register level: ‘ num2str（k） ‘...‘]）;

        % downsample
        scale = 2^-（k-1）;
        Fl = imresize（Fscale）;
        Ml = imresize（Mscale）;
        vxl = imresize（vx*scalescale）;
        vyl = imresize（vy*scalescale）;

        % register
        opt = struct（‘niter‘niter...
            ‘sigma_fluid‘sigma_fluid...
            ‘sigma_diffusion‘sigma_diffusion...
            ‘sigma_i‘sigma_i...
            ‘sigma_x‘sigma_x...
            ‘vx‘vxl ‘vy‘vyl...
            ‘do_display‘do_display ‘do_plotenergy‘1）;
        [Mpsxlsylvxlvyl] = register（FlMlopt）;

        % upsample
        vx = imresize（vxl/scalesize（M））;
        vy = imresize（vyl/scalesize（M））;
    end
    [sxsy] = expfield（vxvy）;
    
end

end

%% Register two images
function [Mpsxsyvxvy] = register（FMopt）

    if nargin<3;  opt = struct（）;  end;
    if ~isfield（opt‘sigma_fluid‘）;      opt.sigma_fluid     = 1.0;              end;
    if ~isfield（opt‘sigma_diffusion‘）;  opt.sigma_diffusion = 1.0;              end;
    if ~isfield（opt‘sigma_i‘）;          opt.sigma_i         = 1.0;              end;
    if ~isfield（opt‘sigma_x‘）;          opt.sigma_x         = 1.0;              end;
    if ~isfield（opt‘niter‘）;            opt.niter           = 250;              end;
    if ~isfield（opt‘vx‘）;              

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     目录           0  2013-01-28 22:49  demons\
     目录           0  2013-01-28 22:49  demons\demons2d\
     文件        6148  2013-01-28 22:49  demons\demons2d\.DS_Store
     目录           0  2013-01-28 22:49  demons\demons2d\data\
     文件         686  2013-01-28 22:49  demons\demons2d\data\heart-110.png
     文件         648  2013-01-28 22:49  demons\demons2d\data\heart-64.png
     文件        1253  2013-01-28 22:49  demons\demons2d\data\lace1.png
     文件        1908  2013-01-28 22:49  demons\demons2d\data\lace2.png
     文件       34678  2013-01-28 22:49  demons\demons2d\data\lenag1.png
     文件       36184  2013-01-28 22:49  demons\demons2d\data\lenag2.png
     文件       41466  2013-01-28 22:49  demons\demons2d\data\lenag3.png
     文件       10906  2013-01-28 22:49  demons\demons2d\data\statue-rio-deformed.png
     文件       10797  2013-01-28 22:49  demons\demons2d\data\statue-rio.png
     文件       14332  2013-01-28 22:49  demons\demons2d\demons.m
     目录           0  2013-01-28 22:49  demons\demons3d\
     文件        6148  2013-01-28 22:49  demons\demons3d\.DS_Store
     文件        1161  2013-01-28 22:49  demons\demons3d\compose.m
     目录           0  2013-01-28 22:49  demons\demons3d\data\
     文件       42198  2013-01-28 22:49  demons\demons3d\data\im1.mat
     文件       39311  2013-01-28 22:49  demons\demons3d\data\im2.mat
     文件        2030  2013-01-28 22:49  demons\demons3d\demons.m
     文件         523  2013-01-28 22:49  demons\demons3d\expfield.m
     文件        1327  2013-01-28 22:49  demons\demons3d\findupdate.m
     文件         808  2013-01-28 22:49  demons\demons3d\imgaussian.m
     文件         887  2013-01-28 22:49  demons\demons3d\iminterpolate.m
     文件         590  2013-01-28 22:49  demons\demons3d\jacobian.m
     文件         970  2013-01-28 22:49  demons\demons3d\piggyback.m
     文件        5504  2013-01-28 22:49  demons\demons3d\register.m
     文件         696  2013-01-28 22:49  demons\demons3d\resize.m
     文件          31  2013-01-28 22:49  demons\demons3d\runme.m
     文件         938  2013-01-28 22:49  demons\demons3d\showgrid.m
............此处省略3个文件信息