Code covered by the BSD License

# multimodality non-rigid demon algorithm image registration

### Dirk-Jan Kroon (view profile)

16 Sep 2008 (Updated )

non-rigid 2D and 3D image registration with demon (fluid) algorithm, extended with modality transf.

[I1_TF,I2_TF]=MutualTransform(I1,I2,SizeKernel,SampleDistance,Bx,By,Bz)
```function [I1_TF,I2_TF]=MutualTransform(I1,I2,SizeKernel,SampleDistance,Bx,By,Bz)
% This function MutualTransform transforms two pictures with different modalities
% into the grey levels of the opposite modality.
% The function divides the images in small overlapping regions, and calculates
% mutual histograms from the regions. A (gaussian) kernel is used so that
% pixels far away from the region counts less in the histogram than nearby
% pixels. The mutual histograms are then used to find the maximum
% correlation between the intensity of the pixels in image 1 and image 2.
% This maximum correlation is then used to paint a new picture of the opposite modality.
%
% [I1_TF,I2_TF]=MutualTransform(I1,I2,SizeKernel,SampleDistance,Bx,By,Bz)
%
% inputs,
%  I1 : Input Image 1
%  I2 : Input Image 2, with a different (MRI) modality than image I1.
%  SizeKernel: The size of the gaussian kernel which is used to make a
%           local mutual histogram.
%  SampleDistance: The distance between the spacing of the
%           kernels  / image regions
%  Bx,By,Bz : Translation images from demon registration in x, y (and z)
%           which deform the images to get a better modality transformation.
%
% outputs,
%  I1_TF : Modality transformed painting of image 1
%  I2_TF : Modality transformed painting of image 2
%
% Example,
%
%
% Bx=zeros(size(Imoving))-38;
% By=zeros(size(Imoving));
%
% [Imoving_TF,Istatic_TF]=MutualTransform(Imoving,Istatic,14,7,Bx,By);
%
% figure,
% subplot(2,2,1), imshow(Imoving)
% subplot(2,2,2), imshow(Istatic)
% subplot(2,2,3), imshow(Imoving_TF)
% subplot(2,2,4), imshow(Istatic_TF)

SizeKernel=floor(SizeKernel/2)*2+1;
SampleDistance=floor(SampleDistance/2)*2+1;

% Histogram number of bins
HistogramBins=255;

% Range of input image
range=getrangefromclass(I1);

% Limit images to the range specified by the input class
I1(I1>range(2))=range(2); I1(I1<range(1))=range(1);
I2(I2>range(2))=range(2); I2(I2<range(1))=range(1);

if(ndims(I1)==2)
% Make a gaussian filter kernel
Hkernel=kernel_gaussian([SizeKernel SizeKernel], SizeKernel/6);

% Make transformed images
if(exist('Bx','var'))
% Estimate the forward transformation field from the backward
% transformation field
[Fx,Fy]=backwards2forwards(Bx,By);

% Make the transformed I1 image
I1_moved=movepixels(I1,Bx,By);

% Make the transformed I2 image
I2_moved=movepixels(I2,Fx,Fy);
else
I1_moved=I1; I2_moved=I2;
end
[I1_TF,I2_TF]=mutual_transform_2d_double(double(I1),double(I2),double(I1_moved),double(I2_moved),double(range),double(HistogramBins),double(Hkernel),double(SampleDistance));
else
% Make a gaussian filter kernel
Hkernel=kernel_gaussian([SizeKernel SizeKernel SizeKernel], SizeKernel/6);

% Make transformed images
if(exist('Bx','var'))
% Estimate the forward transformation field from the backward
% transformation field
[Fx,Fy,Fz]=backwards2forwards(Bx,By,Bz);

% Make the transformed I1 image
I1_moved=movepixels(I1,Bx,By,Bz);

% Make the transformed I2 image
I2_moved=movepixels(I2,Fx,Fy,Fz);
else
I1_moved=I1; I2_moved=I2;
end
if(isa(I1,'double'))
[I1_TF,I2_TF]=mutual_transform_3d_double(double(I1),double(I2),double(I1_moved),double(I2_moved),double(range),double(HistogramBins),double(Hkernel),double(SampleDistance));
else
keyboard;
[I1_TF,I2_TF]=mutual_transform_3d_single(single(I1),single(I2),single(I1_moved),single(I2_moved),single(range),single(HistogramBins),single(Hkernel),single(SampleDistance));
end
end

function h=kernel_gaussian(Hsize,sigma)
siz   = (Hsize-1)/2; std = sigma;
% Make gaussian kernel
if(length(Hsize)==2) % Detect if 2D or 3D kernel
[x,y] = ndgrid(-siz(1):siz(1),-siz(2):siz(2));
arg   = -(x.*x + y.*y)/(2*std*std);
else
[x,y,z] = ndgrid(-siz(1):siz(1),-siz(2):siz(2),-siz(3):siz(3));
arg   = -(x.*x + y.*y + z.*z)/(2*std*std);
end
% No negative values
h = exp(arg);
h(h<eps*max(h(:))) = 0;
% Total sum of kernel must be one
sumh = sum(h(:));
if (sumh ~= 0), h  = h/sumh; end

```