multimodality non-rigid demon algorithm image registration

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

Description

This function will perform demon registration which is an type of fast non-rigid fluid like registration between two 2D or 3D images. Registration between different (MRI) modalities is also supported, through a function which transform one image modality so it looks likes the modality of the second image.

The demon registration is described by the paper of Thirion 1998 and extended by Cachier 1999 and He Wang 2005.

Basic algorithm: On each pixel a velocity (movement) is defined with use of the intensity differences and gradient information. This velocity field is smoothed by an Gaussian, and iteratively used to transform the moving image, and register on to the static image. (Easy to understand code example in file basic_demon_example.m)

Instead of using the basic equations for the "demonregistration" function, we have rewritten it to be used by an limit memory BFGS optimizer in an iterative and multi-resolution way, with also support of diffusion regularization. (see also Tom Vercauteren et al. "Non-parametric Diffeomorphic Image..." )

Transforming one modality into the fake modality of the other image is done with use of 2D mutual histograms between regions of both images, and choosing the grey values which have the highest correlation. D. Kroon et al. "MRI Modality Transformation in Demon Registration" (MutualTransform.m)

See the screenshot for an example result.

Usage:
The functions register_images.m and register_volumes.m are easy to use, and will fit most applications.
First compile the c-code : compile_c_files.m.

Notes:
- All the mex-code is multi-threaded and is tested on both Windows and Linux.
- Please leave comments and report bugs

• Required Products: Image Processing Toolbox
• MATLAB release: MATLAB 7.6 (R2008a)

[Dirk-Jan Kroon]

There is a small bug that images with odd valued dimensions now gives an error. Will be updated tomorrow...

[K K]

Hi, seems it is good. But I cannot complie it..is it due to the Matlab version? I have R2006a and Matlab 6.5. both cannot complie. Thanks a lot!

[K K]

I try to complie it with other machine, thanks. Would like to know if you have version for color image~ Thanks~

[Dirk-Jan Kroon]

* K K,
Color image registration example:
C1=im2double(imread('lena1.jpg'));
    C2=im2double(imread('lena2.jpg'));
    [I3,Tx,Ty]=demonregistration(rgb2gray(C1),rgb2gray(C2),[],false);
    C3 = zeros(size(C1));
    C3(:,:,1)=movepixels(C1(:,:,1),Tx,Ty);
    C3(:,:,2)=movepixels(C1(:,:,2),Tx,Ty);
    C3(:,:,3)=movepixels(C1(:,:,3),Tx,Ty);
    imshow(C3,[]); title('Registered image');

[K K]

thanks a lot. I know what you mean. But what I mean is change the Idiff as something like this:
Idiff_color=M_color-S_color;
Idiff = Idiff_color(:,1)+Idiff_color(:,2)+Idiff_color(:,3) ; Idiff = Idiff/3 ;
where M_color is color version of M previous obtained.

[Franz Rudolf]

I am a student and want to learn something about nonrigid registration techniques. Can anybody tell me which transformation model it uses (is the Thin Plate Transformation the same as the B Spline?)

[K K]

I test the code by adding the following:

% Register image 1 to image 2
[I3,Tx,Ty]=demonregistration(I1,I2,[],true);

[h,w] = size(Tx) ;
Ty = zeros(h,w) ;
Tx = 5*ones(h,w) ;

I3=movepixels(I1,Tx,Ty);

but the result is not what I expected.
The image is shifted up instead of left~~~
Indeed what is Tx,Ty doing? is it a deformation field? forward or backward warping?

[K K]

It seems that the Tx and Ty are indeed swapped???
I have tried the code on a very simple line registration, the "registered image" is correct, however, Tx and Ty seem being swapped...

[Dirk-Jan Kroon]

*K K,
Thanks for your comments. The method uses backward warping. Today, I will look at my code if I accidentally swapped Tx and Ty. This week or the next I will upload a new version, which can be used with any Matlab optimizer.

* Franz Rudolf
This transformation algorithm is called: non-parametric, iconic feature based, using point features.

[Assaf cohen]

Will this be good to register a degraded image due to Atmospheric Turbulence to a reference frame (avarage of say 10 degraded images) and this how to improve the quality of the video?

[Dirk-Jan Kroon]

Please upgrade to version 10b if you experience Matlab crashes with affine transformation in version 10.

[Michel Kocher]

Hello, I am using Matlab from a Mac and the files system.h and process.h do not seem to be available. How can I compile movepixels_2d_double.c

[Dirk-Jan Kroon]

*Michel Kocher...
After this weekend i will add a (previous) single threaded version, which works with the Mac computer.

[yichen fan]

Having an error when compile the program:

??? Error using ==> mex
Unable to complete successfully

Error in ==> compile_c_files at 4
    mex(files(i).name,'-v');

Error in ==> basic_demon_example at 7
compile_c_files

Is that a common error?
I am using Matlab ver 7.1.0.246(R14) SP3

[Dirk-Jan Kroon]

*yichen fan
I have never seen this error.
Please have a look at the help of Mex : http://www.mathworks.com/access/helpdesk/help/techdoc/index.html?/access/helpdesk/help/techdoc/ref/mex.htm

Maybee Mex or the LCC compiler is not installed with Matlab.
You can also trie the single threaded c-code files...

[nor ki]

works pretty fine!

[Alper Yaman]

Can you add single thread codes of the latest version. I am using linux and it gives error.
In addition, as using in linux, compilation of mex files gives error just because of the comment lines. You can avoid such errors by replacing "//" with "/*" and of course add "*/" at the end of the line.

[Cheng Lu]

good job!

[Geoff]

 

[jichao zhao]

Could you please make good examples for 3D rigid and nonrigid registration which demonstrate how to set parameters for registration? thanks, jichao

[jichao zhao]

Hi, Sir, your codes are definitely good. When I use your 3D example in register_volumes, the results shown by showcs3
looks weird. Also the values range from -0.2 to 0.3 for example.
Can you explain why it. On the other hand, when I used Ireg2=movepixels(Istatic,Fx,Fy,Fz); showcs3(Ireg2);
Ireg2, however, turns out to be surprising good. thanks.

[Adeola]

Hi Sir, your work is very good & relevant. Im still a beginner in Matlab environment. Im working on a project that is related to your work. I tried running the codes but for the basic_demon_example, it gave me an error message at line 50, that is, [M=movepixels(I1,Tx,Ty);]. Can you pls help me? Hoping to hear from you soonest.

[gx]

Good work. Thank you.

One question: When you determine the bin size in function "registration_error_mutual_info", you use:

bins=round(numel(V)^(1/ndims(V)));

May I know what this formulation is obtained?

[gx]

I meant how the equation "bins=round(numel(V)^(1/ndims(V)));" is obtained? Is it based on some existing publication?

[Dirk-Jan Kroon]

*gx,

The equation "bins=round(numel(V)^(1/ndims(V)));" is not from literature. You can also use a fix bin size of 255, or so.

- Mutual information only works if there are overlaps in the histogram bins.
- More histogram bins will give better registration accuracy.
Thus this equation gives a good trade-off based on the number of available pixels ...

[Rob Campbell]

 

[Rob Campbell]

 

[Rob Campbell]

oops, two ratings. But I do like the program!

[Yin Yang]

Hello, I know it may be a stupid, but I would like to find a function which could compute a pixel say, (x1, y1) in the moving image 's location in the static image, say (x2, y2). It seems to me that the (x1, y1) and (x2, y2) are related through the so-called transformation . Could you explain the detailed relationship among them and what exactly the direction of x, y and z in your code? (from left to right, up to down and top tp buttom?) Thanks

[LI guang]

I ran the example of multi-registration.but the result is not correct.The pixels of floating image turned to be choas when I ran that multi-registration.So is there something wrong about the multi-registration codes?

[LI guang]

sorry i made a mistake .The result in multi-registration is also good

[Toni]

Just a quick question please!

So the pixel value from the moving image is transformed into the deformed moving one and the new position is e.g. x + Tx(x) , y + Ty(y), right? Does this mean that the old image is just updated with these values? what happens to the old pixel values where the pixel was transformed from and is the 'transformed to' pixel value the new value + the old one which was there already? And where does the interpolation take place?

Too much confusion....

[Vijay]

This is an excellent algorithm... I have a quick question. Is this algorithm affected by noise in the image? Will it make sense to denoise image before doing image registration?

Thanks,
VIjay

[Dirk-Jan Kroon]

*Vijay
Denoising the data before registration is not really needed, but when done with an edge preserving filter, it will slightly improve the final results.

[yang xiaomei]

good job.thanks.
I want to try the example in the demons energy,but there is an erro message for the fminsd,I can not find the information of fminsd. look forward to receiving your reply.thanks
                    xiaomei

[emilie]

The file ''movepixels.m'' does not exist in the zipped folder that I downloaded. So, I searched it in the net but the version that I found does not fit to the existing 3D image registration :(
Can you please tell me where I can find the 3D version of the file ''movepixels.m''.
Thanks.

[Dirk-Jan Kroon]

*emilie

The functionmovepixels.m is present in the sub-folder functions_nonrigid. To make 3D work, you have to run compile_c_files ..., or compile the c-code with the mex function.

[sajan]

hi Dirk,
Ive been using demon's off late. Had a question on the forward and backward transformation functions and its use in conjuction with movepixels.m . Here is what I have: A deformation map (D) derived by registering a moving object (M) to a static object (S). I would now like to use D^-1 and M and solve for S.

Ive tried using movepixels.m and the bacward2forward.m, but as I see it, movepixels.m only considers static images as its inputs, which is not what I desire. Your comments and suggestions on this would be highly appreciated.

Thanks

[Cheng Lu]

What a perfect work~~~ I try it and foud that the gradient of the image is calculated as [My,Mx]=gradient(M); But why you want to inverse the x and y? Since normally we use [Mx,My]=gradient(M); Maybe it is related to your 'Movepixels' function?

Your implemntation is based on the additive Demons, do you have any idea to add a compositive Demons? How do we calculate the composition of fields in Matlab? Thanks.

[arghavan]

Hi Dirk
i have 3d( face ) laser scan database .
by registration i want to establish correspondenc for all points of the face, in order to after
 bringing them into correspondence linear combinations of scans,are again faces. Could you please make good examples for 3D nonrigid registration which demonstrate how to set parameters for registration? thanks, arghavan.

[fan]

so, I think that it can be use less memory if you break the large arry into many little arrays.

[fan]

I tried to do 3D nonrigid registration use this program, but I found that it use very large of memory. I tried to run the program in a 64byte OS, memory is 3G , the size of Moving image and Static Image are both 256*256*128. but it run to "out of memory". I found that there a very large of memory not released, I don't know how much memory it should use,so you should do much more in memory management.

[Dirk-Jan Kroon]

*fan,
When you optimize with the quasi newton method, you need memory space for 10 x the amount of unknowns.
In your case 256 * 256 * 128 * 3 (x y z) *10 = 251 658 240 values in the Estimated Hessian, and to store the transformation fields of x,y,z and both images 256 * 256 * 128 * 5 = 41 943 040 values.
In total you store 293 601 280 values, if all those values are in double we have to store 8x this amount in bytes, is 2.2GB ...
Thus be sure to also use Matlab 64bit ..

[fan]

Thank you!
so it must use a 1920M continuous memory space to create a huge array, I think ihis is the reason of "out of memory".
I tested you 3D Affine registration, but the registered Image(Ireg) is just identical to the Imoving image, and the deformation field Fx,Fy,Fz is all zeros. is there any mistake in your program?

[Luca]

 

[yang liu]

??? Undefined function or method 'movepixels' for input arguments of type 'double'.

Error in ==> basic_demon_example at 50
        M=movepixels(I1,Tx,Ty);

Error in ==> run at 57
          evalin('caller', [s ';']);

After compling, above error was presented. Did I make some mistake?

[yang liu]

Sorry,it's my fault. With the step of "Set Path", the program run normally.
Thank you for excellent job.

[Heng]

Hi. I'm using your program to register MRI T1 and T2 images. I used [Ireg,Bx,By] = register_images(Imoving,Istatic);
it ran well at first, then it gave me "Undefined function or variable "resizepercentage".

Error in register_images (line 323)
    if(resizepercentage~=1)"
this error.

Could you tell me what is wrong?
Thanks.