Code covered by the BSD License  

Highlights from
B-spline Grid, Image and Point based Registration

image thumbnail

B-spline Grid, Image and Point based Registration

by

 

26 May 2008 (Updated )

B-spline registration of two 2D / 3D images or corrsp. points, affine and with smooth b-spline grid.

[O_error, O_grad]=bspline_registration_gradient(O_grid,sizes,Spacing,I1,I2,options,MaskI1,MaskI2,Points1,Points2,PStrength)
function [O_error, O_grad]=bspline_registration_gradient(O_grid,sizes,Spacing,I1,I2,options,MaskI1,MaskI2,Points1,Points2,PStrength)
% Function registration_gradient. This function will calculate a registration
% error value and gradient after b-spline non-rigid registration
% of two images / volumes.
%
% [error, errorgrad]=bspline_registration_gradient(grid,sizes,Spacing,I1,I2,options,MaskI1,MaskI2)
%
% inputs,
%   grid: b-spline control grid created by make_init_grid.m, (and reshaped
%         to one long vector.)
%   sizes: sizes need tot reshape the grid to matrix format.
%   Spacing: The spacing in x,y (and z) direction of the b-spline grid
%           knots
%   I1 and I2: The input images, of which I1 is transformed.
% (optional)
%   options: Struct with options
%       ->type: Type of image similarity(error) measure used
%               (see image_difference.m) (default 'sd')
%       ->penaltypercentage Percentage of penalty smoothness (bending energy
%               of thin sheet of metal) (default 0.01)
%       ->step: Delta step used for error gradient (default 0.01)
%       ->centralgrad: Calculate derivatives with central instead of forward
%               gradient (default true)
%       ->interpolation Type of interpolation used, linear (default) or
%               cubic.
%       ->scaling: Scaling of dimensions 1x2 or 1x3 (like mm/px), to
%               be used in smoothness penalty
%       ->verbose : Display information (default false)
%   MaskI1: Image/volume which transformed in the same way as I1 and
%           is multiplied with the individual pixel errors
%         before calculation of the te total (mean) similarity error.
%   MaskI2: Also a Mask but is used  for I2
%
% outputs,
%   error: The registration error value
%   errorgrad: The registration error gradient
%
% Note :
%   Control points only influence their neighbor region. Thus when calculating
%   the gradient with finited difference, multiple control points (with a spacing
%   between them of 4) can be moved in one finite difference deformation step, and
%   only the pixels in the neigbourhood of a control point  are used to calculate
%   the finite difference.
%
% example,
%   I1=im2double(imread('lenag1.png'));
%   I2=im2double(imread('lenag2.png'));
%   O_trans=make_init_grid([32 32],size(I1)); sizes=size(O_trans);
%   O_trans = fminsd(@(x)bspline_registration_gradient(x,sizes,[32 32],I1,I2,'sd'), O_trans(:), ...
%               optimset('Display','iter','GradObj','on','MaxIter',20,'DiffMinChange',0.1));
%   Icor=bspline_transform(reshape(O_trans,sizes),I1,Spacing);
%   figure, imshow(I1), figure, imshow(I2), figure, imshow(Icor)
%
%
% This function is written by D.Kroon University of Twente (April 2009)

% Check/set input options
defaultoptions=struct('type','sd', 'penaltypercentage',0.01,'step',0.01,'centralgrad',true,'interpolation','linear','scaling',[1 1 1],'verbose',false);
if(~exist('options','var')),
    options=defaultoptions;
else
    tags = fieldnames(defaultoptions);
    for i=1:length(tags)
        if(~isfield(options,tags{i})),  options.(tags{i})=defaultoptions.(tags{i}); end
    end
    if(length(tags)~=length(fieldnames(options))),
        warning('BsplineRegistrationGradient:unknownoption','unknown options found');
    end
end
% Check if there is any Mask input
if(~exist('MaskI1','var')), MaskI1=[]; end
if(~exist('MaskI2','var')), MaskI2=[]; end
% Check on points input
if(~exist('Points1','var')), Points1=[]; end
if(~exist('Points2','var')), Points2=[]; end
if(~exist('PStrength','var')), PStrength=[]; end

% Type of image similarity(error) measure used
type=options.type;
% Percentage of penalty smoothness (bending energy of thin sheet of metal)
penaltypercentage=options.penaltypercentage;
% Delta step used for error gradient
step=options.step;
% Central gradient or (faster) forward gradient
centralgrad=options.centralgrad;
% Interpolation used linear or cubic
if(strcmpi(options.interpolation(1),'l')), interpolation_mode=0; else interpolation_mode=2; end
% Convert Grid vector to grid matrix
O_grid=reshape(O_grid,sizes);

if(penaltypercentage>0)
    % Calculate penalty smoothness (bending energy of thin sheet of metal)
    if ( nargout > 1 )
        % If penalty gradient needed, also determine the gradient.
        [SO_error, SO_grad]=penalties_smoothness(O_grid,size(I1),options.scaling);
    else
        SO_error=penalties_smoothness(O_grid,size(I1),options.scaling);
    end
end

% Transform MaskI1 with the b-spline transformation.
if(~isempty(MaskI1)), MaskI1=bspline_transform(O_grid,MaskI1,Spacing); end
% Calculate a combine masked for excluding some regions from image
if(~isempty(MaskI1)), if(~isempty(MaskI2)), Mask=MaskI1.*MaskI2; else Mask=MaskI1; end
elseif(~isempty(MaskI2)), Mask=MaskI2;
else Mask=[];
end

% Quick SSD error and error gradient calculation
if(strcmpi(type,'sd')&&isempty(Points1)&&isempty(Mask));
    try
        if ( nargout == 1 )
            if(size(I1,3)<4)
                O_error=bspline_error_2d_double(double(O_grid(:,:,1)),double(O_grid(:,:,2)),double(I1),double(I2),double(Spacing(1)),double(Spacing(2)));
            else
                if(isa(I1,'double'))
                    O_error=bspline_error_3d_double(double(O_grid(:,:,:,1)),double(O_grid(:,:,:,2)),double(O_grid(:,:,:,3)),double(I1),double(I2),double(Spacing(1)),double(Spacing(2)),double(Spacing(3)));
                else
                    O_error=bspline_error_3d_single(single(O_grid(:,:,:,1)),single(O_grid(:,:,:,2)),single(O_grid(:,:,:,3)),single(I1),single(I2),single(Spacing(1)),single(Spacing(2)),single(Spacing(3)));
                end
            end
            % Add penalty to total error
            if(penaltypercentage>0), O_error=O_error+SO_error*penaltypercentage; end
        else
            if(size(I1,3)<4)
                [O_error,O_grad]=bspline_error_2d_double(double(O_grid(:,:,1)),double(O_grid(:,:,2)),double(I1),double(I2),double(Spacing(1)),double(Spacing(2)));
            else
                if(isa(I1,'double'))
                    [O_error,O_grad]=bspline_error_3d_double(double(O_grid(:,:,:,1)),double(O_grid(:,:,:,2)),double(O_grid(:,:,:,3)),double(I1),double(I2),double(Spacing(1)),double(Spacing(2)),double(Spacing(3)));
                else
                    [O_error,O_grad]=bspline_error_3d_single(single(O_grid(:,:,:,1)),single(O_grid(:,:,:,2)),single(O_grid(:,:,:,3)),single(I1),single(I2),single(Spacing(1)),single(Spacing(2)),single(Spacing(3)));
                end
            end
            % Add penalty to total error
            if(options.verbose)
                disp(['Error' num2str(O_error)]);
                disp(['Smoothness Error' num2str(SO_error)]);
                disp(['Abs Mean Error Gradient' num2str(mean(abs(O_grad(:))))]);
                disp(['Abs Max  Error Gradient' num2str(max(abs(O_grad(:))))]);
                disp(['Abs Max  Smoothness Gradient' num2str(max(abs(SO_grad(:))))]);
                disp(['Abs Mean Smoothness Gradient' num2str(mean(abs(SO_grad(:))))]);
            end


            if(penaltypercentage>0), O_error=O_error+SO_error*penaltypercentage;
                O_grad=O_grad+SO_grad*penaltypercentage;
            end
        end
        O_grad(isnan(O_grad))=0;
        return
    catch
    end
end

% Transform the image with the B-spline grid
if(~isempty(Points1))
    % In this case we need the transformation fields to get the movements
    % of the landmarks
    if(size(I1,3)<4)
        [I_init,B]=bspline_transform(O_grid,I1,Spacing,interpolation_mode);
        Bx=B(:,:,1);
        By=B(:,:,2);
    else
        [I_init,B]=bspline_transform(O_grid,I1,Spacing,interpolation_mode);
        Bx=B(:,:,:,1);
        By=B(:,:,:,2);
        Bz=B(:,:,:,3);
    end
else
    Bx=[]; By=[]; Bz=[];
    I_init=bspline_transform(O_grid,I1,Spacing,interpolation_mode);
end

% Calculate the current registration error
if(strcmpi(type,'mi'))
    O_uniform=make_init_grid(Spacing,size(I1));
    O_error=0;
    MaskNum=numel(I1);
    if(size(I1,3)<4)
        for i=1:size(O_uniform,1)
            for j=1:size(O_uniform,2)
                [regAx,regAy,regBx,regBy]=regioninfluenced2D(i,j,O_uniform,size(I1));
                % Mask some regions from the registration error measure
                if(~isempty(Mask)), MaskRegion=Mask(regAx:regBx,regAy:regBy); else MaskRegion=[]; end
                % Determine the registration error in the region
                O_error=O_error+image_difference(I_init(regAx:regBx,regAy:regBy),I2(regAx:regBx,regAy:regBy),type,MaskRegion,MaskNum);
            end
        end
    else
        for i=1:size(O_uniform,1)
            for j=1:size(O_uniform,2)
                for k=1:size(O_uniform,3)
                    [regAx,regAy,regAz,regBx,regBy,regBz]=regioninfluenced3D(i,j,k,O_uniform,size(I1));
                    % Mask some regions from the registration error measure
                    if(~isempty(Mask)), MaskRegion=Mask(regAx:regBx,regAy:regBy,regAz:regBz); else MaskRegion=[]; end
                    % Determine the registration error in the region
                    O_error=O_error+image_difference(I_init(regAx:regBx,regAy:regBy,regAz:regBz),I2(regAx:regBx,regAy:regBy,regAz:regBz),type,MaskRegion,MaskNum);
                end
            end
        end
    end
    O_error=O_error/numel(O_grid);
else
    O_error=image_difference(I_init,I2,type,Mask);
end

if(options.verbose&&(nargout > 1 ))
    disp(['Error' num2str(O_error)]);
    disp(['Smoothness Error' num2str(SO_error)]);
end

% Add penalty to total error
if(penaltypercentage>0), O_error=O_error+SO_error*penaltypercentage; end

% Add the Point error
if(~isempty(Points1))
    reg_start=[1 1 1]; reg_end=size(I1);
    if(size(I1,3)<4)
        O_error=O_error+point_error(I1,reg_start,reg_end,Points1,Points2,PStrength,Bx,By);
    else
        O_error=O_error+point_error(I1,reg_start,reg_end,Points1,Points2,PStrength,Bx,By,Bz);
    end
end

%
% Code below is only needed, when the error gradient is asked by the optimizer
%

% If gradient needed, also determine the gradient.
if ( nargout > 1 )
    if(size(I1,3)<4)
        O_grad=error_gradient2d(Spacing,I1,I2,I_init,O_grid,Mask,Points1,Points2,PStrength,Bx,By,interpolation_mode,step,centralgrad,type);
    else
        O_grad=error_gradient3d(Spacing,I1,I2,I_init,O_grid,Mask,Points1,Points2,PStrength,Bx,By,Bz,interpolation_mode,step,centralgrad,type);
    end
    
    if(strcmpi(type,'mi'))
        O_grad=O_grad/numel(O_grad);
    end
    
    if(options.verbose)
        disp(['Abs Mean Error Gradient' num2str(mean(abs(O_grad(:))))]);
        disp(['Abs Max  Error Gradient' num2str(max(abs(O_grad(:))))]);
        disp(['Abs Max  Smoothness Gradient' num2str(max(abs(SO_grad(:))))]);
        disp(['Abs Mean Smoothness Gradient' num2str(mean(abs(SO_grad(:))))]);
    end

    % Add smoothness penalty gradient
    if(penaltypercentage>0), O_grad=O_grad+SO_grad*penaltypercentage; end
    O_grad=O_grad(:);
end

function [regAx,regAy,regBx,regBy]=regioninfluenced2D(i,j,O_uniform,sizeI)
% Calculate pixel region influenced by a grid node
irm=i-2; irp=i+2;
jrm=j-2; jrp=j+2;
irm=max(irm,1); jrm=max(jrm,1);
irp=min(irp,size(O_uniform,1)); jrp=min(jrp,size(O_uniform,2));

regAx=O_uniform(irm,jrm,1); regAy=O_uniform(irm,jrm,2);
regBx=O_uniform(irp,jrp,1); regBy=O_uniform(irp,jrp,2);

if(regAx>regBx), regAxt=regAx; regAx=regBx; regBx=regAxt; end
if(regAy>regBy), regAyt=regAy; regAy=regBy; regBy=regAyt; end

regAx=max(regAx,1); regAy=max(regAy,1);
regBx=max(regBx,1); regBy=max(regBy,1);
regAx=min(regAx,sizeI(1)); regAy=min(regAy,sizeI(2));
regBx=min(regBx,sizeI(1)); regBy=min(regBy,sizeI(2));

function [regAx,regAy,regAz,regBx,regBy,regBz]=regioninfluenced3D(i,j,k,O_uniform,sizeI)
% Calculate pixel region influenced by a grid node
irm=i-2; irp=i+2;
jrm=j-2; jrp=j+2;
krm=k-2; krp=k+2;
irm=max(irm,1); jrm=max(jrm,1); krm=max(krm,1);
irp=min(irp,size(O_uniform,1)); jrp=min(jrp,size(O_uniform,2)); krp=min(krp,size(O_uniform,3));

regAx=O_uniform(irm,jrm,krm,1);
regAy=O_uniform(irm,jrm,krm,2);
regAz=O_uniform(irm,jrm,krm,3);
regBx=O_uniform(irp,jrp,krp,1);
regBy=O_uniform(irp,jrp,krp,2);
regBz=O_uniform(irp,jrp,krp,3);

if(regAx>regBx), regAxt=regAx; regAx=regBx; regBx=regAxt; end
if(regAy>regBy), regAyt=regAy; regAy=regBy; regBy=regAyt; end
if(regAz>regBz), regAzt=regAz; regAz=regBz; regBz=regAzt; end

regAx=max(regAx,1); regAy=max(regAy,1); regAz=max(regAz,1);
regBx=max(regBx,1); regBy=max(regBy,1); regBz=max(regBz,1);
regAx=min(regAx,sizeI(1)); regAy=min(regAy,sizeI(2)); regAz=min(regAz,sizeI(3));
regBx=min(regBx,sizeI(1)); regBy=min(regBy,sizeI(2)); regBz=min(regBz,sizeI(3));



function O_grad=error_gradient2d(Spacing,I1,I2,I_init,O_grid,Mask,Points1,Points2,PStrength,Bx,By,interpolation_mode,step,centralgrad,type)
O_uniform=make_init_grid(Spacing,size(I1));
O_grad=zeros(size(O_grid));
for zi=0:3,
    for zj=0:3,
        % The variables which will contain the controlpoints for
        % determining a central registration error gradient
        O_gradpx=O_grid; O_gradpy=O_grid;
        if(centralgrad)
            O_gradmx=O_grid; O_gradmy=O_grid;
        end
        
        %Set grid movements of every fourth grid node.
        for i=(1+zi):4:size(O_grid,1),
            for j=(1+zj):4:size(O_grid,2),
                O_gradpx(i,j,1)=O_gradpx(i,j,1)+step;
                O_gradpy(i,j,2)=O_gradpy(i,j,2)+step;
                if(centralgrad)
                    O_gradmx(i,j,1)=O_gradmx(i,j,1)-step;
                    O_gradmy(i,j,2)=O_gradmy(i,j,2)-step;
                end
            end
        end
        
        % Do the grid b-spline transformation for movement of nodes to
        % left right top and bottem.
        if(~isempty(Points1))
            [I_gradpx,B]=bspline_transform(O_gradpx,I1,Spacing,interpolation_mode);
            Bx_px=B(:,:,1); By_px=B(:,:,2);
            [I_gradpy,B]=bspline_transform(O_gradpy,I1,Spacing,interpolation_mode);
            Bx_py=B(:,:,1); By_py=B(:,:,2);
            if(centralgrad)
                [I_gradmx,B]=bspline_transform(O_gradmx,I1,Spacing,interpolation_mode);
                Bx_mx=B(:,:,1); By_mx=B(:,:,2);
                [I_gradmy,B]=bspline_transform(O_gradmy,I1,Spacing,interpolation_mode);
                Bx_my=B(:,:,1); By_my=B(:,:,2);
            end
        else
            I_gradpx=bspline_transform(O_gradpx,I1,Spacing,interpolation_mode);
            I_gradpy=bspline_transform(O_gradpy,I1,Spacing,interpolation_mode);
            if(centralgrad)
                I_gradmx=bspline_transform(O_gradmx,I1,Spacing,interpolation_mode);
                I_gradmy=bspline_transform(O_gradmy,I1,Spacing,interpolation_mode);
            end
        end
        
        for i=(1+zi):4:size(O_grid,1),
            for j=(1+zj):4:size(O_grid,2),
                
                % Calculate pixel region influenced by a grid node
                [regAx,regAy,regBx,regBy]=regioninfluenced2D(i,j,O_uniform,size(I1));
                MaskNum=numel(I1);
                % Mask some regions from the registration error
                % measure
                if(~isempty(Mask))
                    MaskRegion=Mask(regAx:regBx,regAy:regBy);
                else
                    MaskRegion=[];
                end
                
                % Determine the registration error in the region
                E_gradpx=image_difference(I_gradpx(regAx:regBx,regAy:regBy),I2(regAx:regBx,regAy:regBy),type,MaskRegion,MaskNum);
                E_gradpy=image_difference(I_gradpy(regAx:regBx,regAy:regBy),I2(regAx:regBx,regAy:regBy),type,MaskRegion,MaskNum);
                if(centralgrad)
                    E_gradmx=image_difference(I_gradmx(regAx:regBx,regAy:regBy),I2(regAx:regBx,regAy:regBy),type,MaskRegion,MaskNum);
                    E_gradmy=image_difference(I_gradmy(regAx:regBx,regAy:regBy),I2(regAx:regBx,regAy:regBy),type,MaskRegion,MaskNum);
                else
                    E_grid=image_difference(I_init(regAx:regBx,regAy:regBy),I2(regAx:regBx,regAy:regBy),type,MaskRegion,MaskNum);
                end
                
                if(~isempty(Points1))
                    E_gradpx=E_gradpx+point_error(I1,[regAx regAy],[regBx regBy],Points1,Points2,PStrength,Bx_px,By_px);
                    E_gradpy=E_gradpy+point_error(I1,[regAx regAy],[regBx regBy],Points1,Points2,PStrength,Bx_py,By_py);
                    if(centralgrad)
                        E_gradmx=E_gradmx+point_error(I1,[regAx regAy],[regBx regBy],Points1,Points2,PStrength,Bx_mx,By_mx);
                        E_gradmy=E_gradmy+point_error(I1,[regAx regAy],[regBx regBy],Points1,Points2,PStrength,Bx_my,By_my);
                    else
                        E_grid=E_grid+point_error(I1,[regAx regAy],[regBx regBy],Points1,Points2,PStrength,Bx,By);
                    end
                end
                
                
                % Calculate the error registration gradient.
                if(centralgrad)
                    O_grad(i,j,1)=(E_gradpx-E_gradmx)/(step*2);
                    O_grad(i,j,2)=(E_gradpy-E_gradmy)/(step*2);
                else
                    O_grad(i,j,1)=(E_gradpx-E_grid)/step;
                    O_grad(i,j,2)=(E_gradpy-E_grid)/step;
                end
            end
        end
    end
end

function O_grad=error_gradient3d(Spacing,I1,I2,I_init,O_grid,Mask,Points1,Points2,PStrength,Bx,By,Bz,interpolation_mode,step,centralgrad,type)
O_uniform=make_init_grid(Spacing,size(I1));
O_grad=zeros(size(O_grid));
for zi=0:(4-1),
    for zj=0:(4-1),
        for zk=0:(4-1),
            O_gradpx=O_grid; O_gradpy=O_grid; O_gradpz=O_grid;
            if(centralgrad)
                O_gradmx=O_grid; O_gradmy=O_grid; O_gradmz=O_grid;
            end
            
            %Set grid movements of every fourth grid node.
            for i=(1+zi):4:size(O_grid,1),
                for j=(1+zj):4:size(O_grid,2),
                    for k=(1+zk):4:size(O_grid,3),
                        O_gradpx(i,j,k,1)=O_gradpx(i,j,k,1)+step;
                        O_gradpy(i,j,k,2)=O_gradpy(i,j,k,2)+step;
                        O_gradpz(i,j,k,3)=O_gradpz(i,j,k,3)+step;
                        if(centralgrad)
                            O_gradmx(i,j,k,1)=O_gradmx(i,j,k,1)-step;
                            O_gradmy(i,j,k,2)=O_gradmy(i,j,k,2)-step;
                            O_gradmz(i,j,k,3)=O_gradmz(i,j,k,3)-step;
                        end
                    end
                end
            end
            
            % Do the grid b-spline transformation for movement of nodes to
            % left right top and bottem.
            if(~isempty(Points1))
                [I_gradpx,B]=bspline_transform(O_gradpx,I1,Spacing,interpolation_mode);
                Bx_px=B(:,:,:,1); By_px=B(:,:,:,2);  Bz_px=B(:,:,:,3);
                [I_gradpy,B]=bspline_transform(O_gradpy,I1,Spacing,interpolation_mode);
                Bx_py=B(:,:,:,1); By_py=B(:,:,:,2);  Bz_py=B(:,:,:,3);
                [I_gradpz,B]=bspline_transform(O_gradpz,I1,Spacing,interpolation_mode);
                Bx_pz=B(:,:,:,1); By_pz=B(:,:,:,2);  Bz_pz=B(:,:,:,3);
                if(centralgrad)
                    [I_gradmx,B]=bspline_transform(O_gradmx,I1,Spacing,interpolation_mode);
                    Bx_mx=B(:,:,:,1); By_mx=B(:,:,:,2);  Bz_mx=B(:,:,:,3);
                    [I_gradmy,B]=bspline_transform(O_gradmy,I1,Spacing,interpolation_mode);
                    Bx_my=B(:,:,:,1); By_my=B(:,:,:,2);  Bz_my=B(:,:,:,3);
                    [I_gradmz,B]=bspline_transform(O_gradmz,I1,Spacing,interpolation_mode);
                    Bx_mz=B(:,:,:,1); By_mz=B(:,:,:,2);  Bz_mz=B(:,:,:,3);
                end
            else
                I_gradpx=bspline_transform(O_gradpx,I1,Spacing,interpolation_mode);
                I_gradpy=bspline_transform(O_gradpy,I1,Spacing,interpolation_mode);
                I_gradpz=bspline_transform(O_gradpz,I1,Spacing,interpolation_mode);
                if(centralgrad)
                    I_gradmx=bspline_transform(O_gradmx,I1,Spacing,interpolation_mode);
                    I_gradmy=bspline_transform(O_gradmy,I1,Spacing,interpolation_mode);
                    I_gradmz=bspline_transform(O_gradmz,I1,Spacing,interpolation_mode);
                end
            end
            
            for i=(1+zi):4:size(O_grid,1),
                for j=(1+zj):4:size(O_grid,2),
                    for k=(1+zk):4:size(O_grid,3),
                        
                        [regAx,regAy,regAz,regBx,regBy,regBz]=regioninfluenced3D(i,j,k,O_uniform,size(I1));
                        MaskNum=numel(I1); %(regBx-regAx+1) * (regBy-regAy+1) * (regBz-regAz+1);
                        
                        % Mask some regions from the registration error
                        % measure
                        if(~isempty(Mask))
                            MaskRegion=Mask(regAx:regBx,regAy:regBy,regAz:regBz);
                        else
                            MaskRegion=[];
                        end
                        
                        % Determine the registration error in the region
                        E_gradpx=image_difference(I_gradpx(regAx:regBx,regAy:regBy,regAz:regBz),I2(regAx:regBx,regAy:regBy,regAz:regBz),type,MaskRegion,MaskNum);
                        E_gradpy=image_difference(I_gradpy(regAx:regBx,regAy:regBy,regAz:regBz),I2(regAx:regBx,regAy:regBy,regAz:regBz),type,MaskRegion,MaskNum);
                        E_gradpz=image_difference(I_gradpz(regAx:regBx,regAy:regBy,regAz:regBz),I2(regAx:regBx,regAy:regBy,regAz:regBz),type,MaskRegion,MaskNum);
                        if(centralgrad)
                            E_gradmx=image_difference(I_gradmx(regAx:regBx,regAy:regBy,regAz:regBz),I2(regAx:regBx,regAy:regBy,regAz:regBz),type,MaskRegion,MaskNum);
                            E_gradmy=image_difference(I_gradmy(regAx:regBx,regAy:regBy,regAz:regBz),I2(regAx:regBx,regAy:regBy,regAz:regBz),type,MaskRegion,MaskNum);
                            E_gradmz=image_difference(I_gradmz(regAx:regBx,regAy:regBy,regAz:regBz),I2(regAx:regBx,regAy:regBy,regAz:regBz),type,MaskRegion,MaskNum);
                        else
                            E_grid=image_difference(I_init(regAx:regBx,regAy:regBy,regAz:regBz),I2(regAx:regBx,regAy:regBy,regAz:regBz),type,MaskRegion,MaskNum);
                        end
                        
                        % Add Error for distance between
                        % corresponding points
                        if(~isempty(Points1))
                            E_gradpx=E_gradpx+point_error(I1,[regAx regAy regAz],[regBx regBy regBz],Points1,Points2,PStrength,Bx_px,By_px,Bz_px);
                            E_gradpy=E_gradpy+point_error(I1,[regAx regAy regAz],[regBx regBy regBz],Points1,Points2,PStrength,Bx_py,By_py,Bz_py);
                            E_gradpz=E_gradpz+point_error(I1,[regAx regAy regAz],[regBx regBy regBz],Points1,Points2,PStrength,Bx_pz,By_pz,Bz_pz);
                            if(centralgrad)
                                E_gradmx=E_gradmx+point_error(I1,[regAx regAy regAz],[regBx regBy regBz],Points1,Points2,PStrength,Bx_mx,By_mx,Bz_mx);
                                E_gradmy=E_gradmy+point_error(I1,[regAx regAy regAz],[regBx regBy regBz],Points1,Points2,PStrength,Bx_my,By_my,Bz_my);
                                E_gradmz=E_gradmz+point_error(I1,[regAx regAy regAz],[regBx regBy regBz],Points1,Points2,PStrength,Bx_mz,By_mz,Bz_mz);
                            else
                                E_grid=E_grid+point_error(I1,[regAx regAy regAz],[regBx regBy regBz],Points1,Points2,PStrength,Bx,By,Bz);
                            end
                        end
                                                
                        % Calculate the error registration gradient.
                        if(centralgrad)
                            O_grad(i,j,k,1)=(E_gradpx-E_gradmx)/(step*2);
                            O_grad(i,j,k,2)=(E_gradpy-E_gradmy)/(step*2);
                            O_grad(i,j,k,3)=(E_gradpz-E_gradmz)/(step*2);
                        else
                            O_grad(i,j,k,1)=(E_gradpx-E_grid)/step;
                            O_grad(i,j,k,2)=(E_gradpy-E_grid)/step;
                            O_grad(i,j,k,3)=(E_gradpz-E_grid)/step;
                        end
                    end
                end
            end
        end
    end
end


function errordis=point_error(I,reg_start,reg_end,Points1,Points2,PStrength,Bx,By,Bz)
if(size(I,3)<4)
    % Points must be inside image
    x1=round(Points2(:,1)); y1=round(Points2(:,2));
    
    % Remove points outside asked region
    check=(x1<reg_start(1))|(x1>reg_end(1))|(y1<reg_start(2))|(y1>reg_end(2));
    x1(check)=[]; y1(check)=[];
    
    % Find transformation of points in transformation field
    x_trans=Bx(sub2ind(size(I),x1,y1))+x1;
    y_trans=By(sub2ind(size(I),x1,y1))+y1;
    
    x2=Points1(:,1); y2=Points1(:,2);
    % Remove points outside asked region
    x2(check)=[]; y2(check)=[]; PStrength(check)=[];
    
    % Normalized distance between transformed and static points
    distance=((x_trans-x2).^2+(y_trans-y2).^2)/(size(I,1).^2+size(I,2).^2);
    
    % The total distance point error, weighted with point strength.
    errordis=sum(distance.*PStrength);
else
    % Points must be inside image
    x1=round(Points2(:,1)); y1=round(Points2(:,2));  z1=round(Points2(:,3));
    % Remove points outside asked region
    check=(x1<reg_start(1))|(x1>reg_end(1))|(y1<reg_start(2))|(y1>reg_end(2))|(z1<reg_start(3))|(z1>reg_end(3));
    x1(check)=[]; y1(check)=[]; z1(check)=[];
    
    % Find transformation of points in transformation field
    x_trans=Bx(sub2ind(size(I),x1,y1,z1))+x1;
    y_trans=By(sub2ind(size(I),x1,y1,z1))+y1;
    z_trans=Bz(sub2ind(size(I),x1,y1,z1))+z1;
    
    x2=Points1(:,1); y2=Points1(:,2); z2=Points1(:,3);
    % Remove points outside asked region
    x2(check)=[]; y2(check)=[]; z2(check)=[]; PStrength(check)=[];
    
    % Normalized distance between transformed and static points
    distance=((x_trans-x2).^2+(y_trans-y2).^2+(z_trans-z2).^2)/(size(I,1).^2+size(I,2).^2+size(I,3).^2);
    
    % The total distance point error, weighted with point strength.
    errordis=sum(distance.*PStrength);
end




Contact us