B-spline Grid, Image and Point based Registration

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

Description

Affine and B-spline grid based registration and data-fitting of two 2D color/grayscale images or 3D volumes or point-data. Registration can be done intensity / pixel based, or landmark / corresponding points based (See OpenSurf), or a combination.

Description Pixel-based registration:
This function is an (enhanced) implementation of the b-spline registration algorithm in D. Rueckert et al. "Nonrigid Registration Using Free-Form Deformations: Application to Breast MR Images". Including the smoothness penalty of Rueckert (thin sheet of metal bending energy), and Jacobian (diffeomorphic) function. Also including, localized normalized mutual information as registration error, allowing the images or volumes to be of a different type/modality for instance a MRI T1 and T2 patient scan.

How it works:
A grid of b-spline control points is constructed which controls the transformation of an input image. An error measure is used to measure the registration error between the moving and static image. The quasi newton Matlab optimizer fminlbfgs (also on Mathworks) is used to move the control points to achieve the optimal registration between both images with minimal registration error.

Usage:
- The function image_registration.m is easy to use, and contains examples in the help, and will fit most applications. (If you want to write your own specialized registration code study the registration examples)
- The function point_registration is fast fitting of a b-spline grid to 2D/3D corresponding points, for landmark based registration.
- There is also the function manually_warp_images which allow control grid changes with the mouse, to get better registration.

First, you need to compile mex / C code with compile_c_files.m. (2D registration also works without mex files but will be slower)

The multi-threaded mex code supports Windows, Linux (and Mac OS?)

Some Features:
- 2-D / 3-D Eulerian strain tensor images can be made from the transformations fields. For example to describe cardiac motion in the images.
- Landmarks can be used for already known corresponding points (for example from Sift). Influence of every landmark on the registration process can be adjust.
- It is possible to register a number of movie frames by using the registration grid of previous two images as initial registration grid of the next two images.
- It is possible to mask parts of the images, to decrease or increase the influence of an image structure on the registration result.

Literature:
- D. Rueckert et al. "Nonrigid Registration Using Free-Form Deformations: Application to Breast MR Images".
- Seungyong Lee, George Wolberg, and Sung Yong Shing, "Scattered Data interpolation with Multilevel B-splines"

note:
- B-spline registration is slower and more complex than demon registration see:
http://www.mathworks.fr/matlabcentral/fileexchange/loadFile.do?objectId=21451
- Why still use B-spine registration?, because the resulting transformation-field corresponds better to real-live deformation than transformation fields from fluid registration.

Please report bugs, successes and questions.

• Required Products: Image Processing Toolbox
• MATLAB release: MATLAB 7.7 (R2008b)

[Dirk-Jan Kroon]

Fast version available May 28 2008.

[Dirk-Jan Kroon]

nonrigid3dfast.m also uploaded soon available...

[fabrizio fasano]

great job. very useful for the MR community

[Ramin Mafi]

Thanks for the codes

[Xu Han]

Very brilliant work! I will try to program it and learn it. Thank you very much

[Perico Palotes]

There must be something wrong on the C code because some error are poped out when I try to compile the c files.

Error bspline_transform_3d_double.c: 37 undeclared identifier `mwSize'
Warning bspline_transform_3d_double.c: 37 Statement has no effect
Error bspline_transform_3d_double.c: 37 syntax error; found `ZOsizex' expecting `;'
Error bspline_transform_3d_double.c: 37 undeclared identifier `ZOsizex'
Error bspline_transform_3d_double.c: 37 undeclared identifier `ZOsizey'
Error bspline_transform_3d_double.c: 37 undeclared identifier `ZOsizez'
Warning bspline_transform_3d_double.c: 37 Statement has no effect
Error bspline_transform_3d_double.c: 38 illegal statement termination
Error bspline_transform_3d_double.c: 38 skipping `const'
Error bspline_transform_3d_double.c: 38 undeclared identifier `dims'
Warning bspline_transform_3d_double.c: 38 Statement has no effect
Warning bspline_transform_3d_double.c: 40 Statement has no effect
Error bspline_transform_3d_double.c: 40 syntax error; found `Osizex' expecting `;'
Error bspline_transform_3d_double.c: 40 undeclared identifier `Osizex'
Error bspline_transform_3d_double.c: 40 undeclared identifier `Osizey'
Error bspline_transform_3d_double.c: 40 undeclared identifier `Osizez'
Warning bspline_transform_3d_double.c: 40 Statement has no effect
Error bspline_transform_3d_double.c: 43 illegal statement termination
Error bspline_transform_3d_double.c: 43 skipping `double'
Error bspline_transform_3d_double.c: 43 undeclared identifier `u'
Error bspline_transform_3d_double.c: 43 undeclared identifier `v'
Error bspline_transform_3d_double.c: 43 undeclared identifier `w'
Warning bspline_transform_3d_double.c: 43 Statement has no effect
Error bspline_transform_3d_double.c: 44 illegal statement termination
Error bspline_transform_3d_double.c: 44 skipping `int'
Error bspline_transform_3d_double.c: 44 undeclared identifier `i'
Error bspline_transform_3d_double.c: 44 too many errors
 
  C:\MATLAB701\BIN\WIN32\MEX.PL: Error: Compile of 'bspline_transform_3d_double.c' failed.
 
??? Error using ==> mex
Unable to complete successfully

[Dirk-Jan Kroon]

Matlab 2007b and newer releases require "mwSize" as type to store data dimensions.
You can simply replace "mwSize" by (a normal) "int", to make the code work with your Matlab version.

Dirk-Jan

[N M]

Thank you for posting this - it looks like it could be very useful for a lot of applications. I'm wondering how difficult it would be to modify it to register two 3D data sets with different voxel dimensions, numbers of slices and prescriptions (i.e., different, x,y,z vectors for the different datasets)? If so, any chance you could you post an example?

Thanks,
N

[deborah alley]

good, although processing time is long.

[Dirk-Jan Kroon]

*deborah alley
This week or the next I will release multi-threaded code, with more speed optimizations and control grid smooth optimization.

[siwo]

seems ur codes just work with Windows:(, could u tell me which C-compiler u used under windows?

[rania hassen]

i got this error when i tried to compile one file
mex rigid_transform_2d_double.c
Error rigid_transform_2d_double.c: 132 syntax error; found `*' expecting `;'
Error rigid_transform_2d_double.c: 132 undeclared identifier `dims'
Error rigid_transform_2d_double.c: 132 type error: pointer expected
Warning rigid_transform_2d_double.c: 132 Statement has no effect
Error rigid_transform_2d_double.c: 134 illegal statement termination
Error rigid_transform_2d_double.c: 134 skipping `double'
Error rigid_transform_2d_double.c: 134 undeclared identifier `mean'
Error rigid_transform_2d_double.c: 134 type error: pointer expected
Error rigid_transform_2d_double.c: 134 illegal expression
Error rigid_transform_2d_double.c: 134 syntax error; found `0' expecting `;'
Warning rigid_transform_2d_double.c: 134 Statement has no effect
Error rigid_transform_2d_double.c: 134 syntax error; found `}' expecting `;'
Warning rigid_transform_2d_double.c: 134 possible usage of mean before definition
Warning rigid_transform_2d_double.c: 132 possible usage of dims before definition
Warning rigid_transform_2d_double.c: 132 local `const int mwSize' is not referenced
Warning rigid_transform_2d_double.c: 125 local `pointer to pointer to void ThreadList' is not referenced
Warning rigid_transform_2d_double.c: 120 local `pointer to pointer to double ThreadArgs2' is not referenced
Warning rigid_transform_2d_double.c: 120 local `pointer to pointer to double ThreadArgs1' is not referenced
Warning rigid_transform_2d_double.c: 117 local `pointer to double M' is not referenced
Warning rigid_transform_2d_double.c: 117 local `pointer to double Iout' is not referenced
Warning rigid_transform_2d_double.c: 117 local `pointer to double Iin' is not referenced
Warning rigid_transform_2d_double.c: 134 empty declaration
Error rigid_transform_2d_double.c: 137 invalid initialization type; found `pointer to pointer to void' expected `int'
Error rigid_transform_2d_double.c: 137 initializer must be constant
Error rigid_transform_2d_double.c: 140 invalid initialization type; found `pointer to pointer to double' expected `int'
Error rigid_transform_2d_double.c: 140 initializer must be constant
Error rigid_transform_2d_double.c: 141 invalid initialization type; found `pointer to pointer to double' expected `int'
Error rigid_transform_2d_double.c: 141 initializer must be constant
Error rigid_transform_2d_double.c: 145 unrecognized declaration
Error rigid_transform_2d_double.c: 145 syntax error; found `!=' expecting `)'
Error rigid_transform_2d_double.c: 145 syntax error; found `!=' expecting `;'
Error rigid_transform_2d_double.c: 145 skipping `!='
Error rigid_transform_2d_double.c: 145 too many errors
 
  C:\PROGRAM FILES\MATLAB\R2006A\BIN\MEX.PL: Error: Compile of 'rigid_transform_2d_double.c' failed.
 
??? Error using ==> mex
Unable to complete successfully

Help Me plz

[Hu Jun]

 Error: File: C:\MATLAB6p5\work\nonrigid_version10b\example_2d_nonrigid_1.m Line: 38 Column: 19
"identifier" expected, "(" found.

[Maria Garcia]

Error using all example_3d :
??? Seven inputs are required.

Error in ==> bspline_transform at 46
             I=bspline_transform_3d_single(single(O(:,:,:,1)),single(O(:,:,:,2)),single(O(:,:,:,3)),single(I),single(Spacing(1)),single(Spacing(2)),single(Spacing(3)),single(black));

Error in ==> get_example_data at 17
V1=bspline_transform(O_trans,V2,Spacing);

Error in ==> example_3d_affine at 10
[I1,I2]=get_example_data;

[Dirk-Jan Kroon]

Thanks for your error report, I forgot to update the single threaded version of bspline_transform_3d_single.c , and will update the code today.

[Am VR]

Dirk-Jan Kroon ,

Thanks a lot for sharing this great piece of work with us. This is a very nice contribution to the Mathworks file exchange repository.

I have one question about your code though. I am looking at your submission on Mar 4. When I run one of your examples, e.g. example_2d_nonrigid_1.m, the results figure shows four images, and I am not sure I understand what is being displayed.

Shouldn't 'transformed image 1' (Icor in your code) be the result of registering Image 1 (I1) onto image 2 (I2)? If so, 'transformed image 1' should look similar to 'image 2', right? However, in your figure 'transformed image 1' resembles more image 1 than image 2. Am I missing something? Is the registration really working in this example?

Thanks a lot!
Best,
AmR

[Dirk-Jan Kroon]

Am VR,

Image 1 is a spherical distored version of Image 2. If you look at the registration result "Icor" you can notice that for instance the hat of Lena has shrinked to the size of the hat in image 2.

I know the nose in the image registration result still looks like the same as in image 1. This is because the grid is to coarse, and that is why image_registration.m uses grid refinement.

[Am VR]

Thanks a lot DIrk-Jan Kroon.

Your script works really great on our machines. One more thing, my research group has access to a cluster and we'd love to run your alignment program on it. I was wondering how difficult would it be to change the .c files so that they can also be compiled on Linux. Any thoughts on this issue?

[Dirk-Jan Kroon]

Am VR,

It is not that difficult to change the c-code files to work multi-threaded on Linux, it is simply replacing the windows functions such as _beginthreadex with pthread_create.

The code is now optimized for "shared caches": one thread does the even and another the odd image lines. With separate caches/processors it will be better to divide the image in several image regions.

When you add a large number of cpu-cores, the memory and the single-threaded Matlab code will be come the bottleneck instead of the cpu power.

---

I hope to release version 12 which uses not 2 but all available cores, and uses a new optimizer FMINLBFGS (also on mathworks) in a few days...

[Brett]

This does not work for arbitrary sized images. I haven't tried above 256x256 yet, but it won't work for smaller ROIs, I get the following error:
??? Attempted to access O_uniform(3,18,1); index out of bounds because size(O_uniform)=[17,17,2].

Error in ==> bspline_registration_gradient at 164
                    regBx=O_uniform(irp,jrp,1); regBy=O_uniform(irp,jrp,2);

Error in ==>
register_images>@(x)bspline_registration_gradient(x,sizes,Spacing/resize_per,imresize(ISmoving,1/resize_per),imresize(ISstatic,1/resize_per),options,MASKmovingsmall,MASKstaticsmall)
at 300
        O_trans = resize_per*
        fminlbfgs(@(x)bspline_registration_gradient(x,sizes,Spacing/resize_per,imresize(ISmoving,1/resize_per),imresize(ISstatic,1/resize_per),options,MASKmovingsmall,MASKstatic
        Error in ==> fminlbfgs>gradient_function at 729
            [fval, grad]=feval(funfcn,reshape(x,data.xsizes));

Error in ==> fminlbfgs at 154
[data,fval,grad]=gradient_function(data.xInitial,funfcn, data, optim);

Error in ==> register_images at 300
        O_trans = resize_per*
        fminlbfgs(@(x)bspline_registration_gradient(x,sizes,Spacing/resize_per,imresize(ISmoving,1/resize_per),imresize(ISstatic,1/resize_per),options,MASKmovingsmall,MASKstatic
        Error in ==> example at 9
  [Ireg,O_trans,Spacing,Bx,By,Fx,Fy] = register_images(Imoving,Istatic);

[Dirk-Jan Kroon]

Dear Brett,

Thank you for the comment, this bug was unfortunately introduced by fixing a bug in the previous update, and is fixed in the update which will be on Mathworks.com today or tomorrow

[MUNIER Mederic]

Hi,

Is the example_3d_affine.m supposed to work on a Matlab 7.1 version ? I encounter some problems while executing the script (segmentation fault).

Thks.

[Dirk-Jan Kroon]

Dear MUNIER Mederic,

I never tested my scripts with Matlab version 7.1, but it should work. You can replace the multi-threaded code files with the single threaded files, maybe that will work?

Does the "segmentation fault", happen every time, you call example_3d_affine?, I have checked and tested the code again but couldn't find a bug...

[Dirk-Jan Kroon]

I've tested the code with Matlab version 7.1.
With this old Matlab version, you must first replace the word "mwsize" (unsupported type) with "int" in all c-code files, and use a newer compiler than the default LCC.

[Abdoul-Nasser DJIDO]

great code, really useful for both starter and advanced ...
working very well on MATLAB 7.6 R2008a

[MUNIER Mederic]

Hi,

Indeedn each time I run example_3d_affine, I got a Segmentation violation error. I can copy the last part of the report error here if it can help someone telling me what's wrong with my version or my manipulation :

Error in ==> gui_mainfcn>local_openfig at 216
    gui_hFigure = openfig(name, singleton, 'auto');

Error in ==> gui_mainfcn at 94
        gui_hFigure = local_openfig(gui_State.gui_Name, gui_SingletonOpt);

Error in ==> showcs3 at 29
    gui_mainfcn(gui_State, varargin{:});

Error in ==> example_3d_affine at 19
showcs3(I1);

I will try to get a newer version of Matlab from my university, hope it will solve the problem.

Thks a lot for trying to help me.

[MUNIER Mederic]

Yeah, in fact it seems to crash almost at the very beginning, when trying to show I1 and I2. Which probably leads to a problem with my version of Matlab.

[charu hans]

Whenever I try to run register_volumes.m on my rgb images.It gives me following error:

??? Out of memory. Type HELP MEMORY for your options.

Error in ==> fminlbfgs at 140
    data.deltaX=zeros(data.numberOfVariables,optim.StoreN);

Error in ==> register_volumes at 272
    O_trans =
    resize_per*fminlbfgs(@(x)bspline_registration_gradient(x,sizes,Spacing/resize_per,imresize3d(ISmoving,1/resize_per,[],'linear','bound'),imresize3d(ISstatic,1/resize_per,[],'linear','bound')

I tried to resize my image to smaller size like 50x50, it still gives me error

[Dirk-Jan Kroon]

*MUNIER Mederic

Your segmentation fault when the affine transformation example shows the 3D data is caused by the fact that GUI's created by new Matlab versions are not 100% compatible with older Matlab versions, even if you set your save settings to an older setting.

Only solution is to not use show3d and use a simple imshow to show a slice, or upgrade to a recent Matlab version

[Ali B]

Dirk-Jan, thanks for sharing your code.
I have a pretty basic question.
I am trying to use the B-spline FFD part (bspline_transform.m) of your code for image morphing. More specifically, I am trying to find out how I can change transformation grid to manipulate a picture.
Transformation grid determines where each of the control points of the new image should be. is that right?
You have mentioned that you've used Fumihiko Ino et al. paper. Just wondering if you have used the same x and y directions mentioned in that paper (meaning (x,y)=(0,0) as the lower left corner of image).
thanks a lot.

[MUNIER Mederic]

Thanks, it perfectly works on a Matlab R2007a version.

Med.

[Dirk-Jan Kroon]

* alib B,

I use the left-top of the image as (0,0) coordinate, and the transformation is done backwards.

[Gerald Lodron Lodron]

is it possible to add a funcitonality that crop the transformed image in movepixels? If I rotate an image via the ffd the images becomes smaller...

[Gerald Lodron Lodron]

is it possible to add a funcitonality that crop the transformed image in movepixels? If I rotate an image via the ffd the images becomes smaller...and I think the result becomes completely wrong at rotations >40 degrees

[Gerald Lodron Lodron]

and I think the result becomes completely wrong at rotations >40 degrees

[Luca]

Dear Dirk
your code is so usefull for my M.s thesis and I,m enjoying it.
but I was wondering if I could just optimize the objective functions using some optimization methods like PSO,...
but when I try this specially in non-rigid transformation I cant do ie well.
infact, I dont know how to use these objective functions(i.e affine_registration_image in affine transform and bspline_registration_gradient in non-rigid)?
would you please post an example of their optimization using PSO?

[aya b]

Hi
I encountered a compilation problem, which seems different than the ones already posted:

Writing library for bspline_transform_3d_single.mexw32
c:\docume~1\admini~1\locals~1\temp\mex_vrmztx\bspline_transform_3d_single.obj .text: undefined reference to '_floorf'
 
  C:\PROGRA~1\MATLAB\R2008A\BIN\MEX.PL: Error: Link of 'bspline_transform_3d_single.mexw32' failed.

I'm using matlab 2008a

Thanks!
Aya

[Dirk-Jan Kroon]

* aya b
Thank you for your comments, I fixed this LCC compiler error in release of today.
(you can also install the last lcc compiler or use another compiler)

[aya b]

It worked - thanks!
Is there help somewhere? or should I use the examples?

[aya b]

Hi
I'm sorry for all these questions, but I encountered another problem - I tried to run the example_nonrigid_3d example and got the following output:
     Iteration Func-count Grad-count f(x) Step-size
         0 1 1 0.048856
         1 4 2 0.0488542 5
         2 7 3 0.0488539 1
         3 10 4 0.0488535 1
         4 13 5 0.0488532 1
    Optimizer Results
        Algorithm Used: limited memory BFGS (L-BFGS)
        Exit message : Number of iterations exceeded options.MaxIter or number of function evaluations exceeded options.FunEvals.
        iterations : 5
        Function Count : 16
        Minimum found : 0.048853
        Intern Time : 0.1239 seconds
        Total Time : 2884.6883 seconds
??? Index exceeds matrix dimensions.

Error in ==> refine_grid at 111
    O_new=O_new(1:size(X,1),1:size(X,2),1:size(X,3),1:3);

Error in ==> example_3d_nonrigid at 51
[O_trans,Spacing]=refine_grid(O_trans,Spacing,size(I1));

Do I need to change any params?
Thanks again,
Aya

[christos]

very nice code... there are some issues though

a)
regarding the masks, could you please confirm what their structure should be? Is it just an image that contains 0s and 1s (where 0s correspond to areas to be masked out)

b)
should the images to be registered be the same size? I also see that at some point in register_images.m you resize the moving image to have the same size with static... fair enough... you don't do this for the masks and I think the system complains about this

c)
In image_difference.m

RE: (lines 243-252)
% Remove unmasked pixels
if(~isempty(Mask)),
    sizes=size(V);
    V=V(:);
    V(~Mask)=[];
    V=reshape(V,sizes);
    U=U(:);
    U(~Mask)=[];
    U=reshape(U,sizes);
end

I don't understand how is it possible to expect, if you remove the masked pixels, to be able to reshape U or V to the same size prior to the pixel removal? A 32x32 image (vectorized 1024x1). Your remove lets say 5 pixels from the vector - now 1019x1, and then you ask it to be reshaped back to 32x32
(my version at least complains)

could you please elaborate on the usages of masks and if they have been tested?

----

thanks alot for your time and effort and for sharing this with us

[Jan Van den Bulcke]

Is it possible to do a straightforward rigid-body registration or can the non-rigid code cope with that?

Jan.

[Jan Van den Bulcke]

Sorry

Yes it can, of course, by reading the help. I apologize.

Jan.

[Dirk-Jan Kroon]

*christos
Sorry a little late, but yesterday I had time to work on my b-spline registration method.

Your 1th question, must masks be zero and ones?:
No the mask is just multiplied with the error image from the similarity measure, before summing up all errors. There is one exception mutual information, in which elements which are zero are simply discarded before calculation of the mutual information value.

2th, must images be the same size? :
Yes, they must be the same size also the masks. I fixed the mask resize bug you mentioned in version 17a which will probably be accepted by Mathworks this evening.

3th, Code image_difference reshape?:
This is an bug, I fixed it in version 17

4th , Mask are they already been tested?:
Only the 2D masks. I developed the masks to register the backgrounds of an angiogram, in which I masked out the coronary artery in the images. After registration of the backgrounds I subtracted them from the original movie to make the coronary artery better visible.

Thank you for your comments, and bug report.

[lajulajay]

Hi, whenever I try using landmarks I get a matrix multiplication error:

*******************************************************************************
??? Error using ==> times
Matrix dimensions must agree.

Error in ==> affine_registration_error>affine_point_error at 255
    errordis=sum(distance.*PStrength);

Error in ==> affine_registration_error>affine_registration_error_2d at 161
    e=e+affine_point_error(M,I1,Points1, Points2, PStrength);

Error in ==> affine_registration_error at 107
    e=affine_registration_error_2d(par,I1,I2,type,O_trans,Spacing,MaskI1,MaskI2,Points1,Points2,PStrength);

Error in ==>
register_images>@(x)affine_registration_error(x,scale,ITmoving,ITstatic,type,O_trans,Spacing,MASKmoving,MASKstatic,Points1t,Points2t,PStrengtht)

Error in ==> fminlbfgs>gradient_function at 729
            [fval, grad]=feval(funfcn,reshape(x,data.xsizes));

Error in ==> fminlbfgs at 154
[data,fval,grad]=gradient_function(data.xInitial,funfcn, data, optim);

Error in ==> register_images at 230
x=fminlbfgs(@(x)affine_registration_error(x,scale,ITmoving,ITstatic,type,O_trans,Spacing,MASKmoving,MASKstatic,Points1t,Points2t,PStrengtht),x,optim);
*******************************************************************************

I am using three landmark points; i.e. Points1 and Points2 are 3 by 2 matrices and PStrength is a 3 by 1 matrix.

Thanks.

[siamak yousefi]

Hi Dirk,
I want to nonrigidly register two .nii volumes, but register_volumes cammand doesn't accept this format. Could you please tell me what formats it accepts.

Siamak

[Dirk-Jan Kroon]

*siamak,
You can find "Tools for NIfTI" on the Mathwork file-exchange. Those tools can be used to load your data in Matlab, and save the registration data afterward ...

[siamak yousefi]

Hi Drik,

I want to use some options in the regsitration. How can I put the options in the command line:
image_r= register_volumes(image1, image2, options);
How can I set options in the command line?
Could you please write for me an example of using Affine registration as an option.

Thanks,
Siamak

[siamak yousefi]

Hi Drik,

I srossed out my previous mistake. It was a mistake in creating Options structure. However, I did a simple rigid registration between to MRI images. When I see the result, Ireg, it is completely corrupted, half part is white, half part is a bright MRI image. Is there any this that I am not setting. I checked the input images, both of them are int16, the same as result. Any help is highly appreciated.

Siamak

[siamak yousefi]

Hi Drik,

I encountered another problem in registration_volumes. When I run this command:
Ireg= register_volumes(vol1, vol1, Options) That:
Options.Registration= 'Affine', all of the voxels of the output volume are zero.

Siamak
;
I am using a 3-D volume and try to

[Dirk-Jan Kroon]

Dear Siamak,

The problem is that you use int16 as datatype, if you convert it to double or single with im2single or im2double, it will probably just work fine.

[siamak yousefi]

Hi Drik,

Thanks for your comment. There is a typo in the help part of image_difference.m file. When I run "help image_difference", in the example section I see:
[t,I]=registration_error(I1, I2, 'sd');
It shold be changed to:
[t,I]=image_difference(I1, I2, 'sd');

Siamak

[siamak yousefi]

Hi Drik,

I have some questions related to register_volumes program:
1. Is it possible to generate the affine matrix or warp file of transformation when registering two images?
2. is it possible to reduce the run time of nonrigid registration using some parameters?
3. For nonrigid registration using a mask, How I can set the "Options.maskMoving". I tried several times using any sort of mask file type with the same dimention, but didn't work.

Siamak Yousefi

[Dirk-Jan Kroon]

* Siamak Yousefi
1 : The (inverse) Affine Matrix is returned as variable "M" if you use the latest version (18)
2 : It is faster to already set:
- the 'similarity' option
- not use all outputs but only [Ireg,Grid,Spacing,M]
- do less refinements option 'MaxRef'
- Gaussian filter your images with imfilter and fspecial to remove noise, before starting the registration.

The fastest: is using no Mask, no Landmarks, in combination with 'sd' as similarity measure, because than it uses optimized gradient mex functions named bspline_error.

3 : Do you use version 18?, otherwise I'll try to find the bug.

[siamak yousefi]

Thanks for you comments. I am using version 18 still have a problem using mask.

Siamak

[Dirk-Jan Kroon]

* siamak yousefi
I have uploaded a new version which I've tested with masks, will probably be online tomorrow.

[siamak yousefi]

hank you very much Drik.

[jichao zhao]

When I used the function "showcs3", my input images are pure black
and white, but when I used matlab built-in fucntion imagesc, I can see my images are gray color instead of pure black and white. Also when I use colormap, I can see the color ranges from 0 to 60. So my question is that why my images are pure black-white when plotted using showcs3 and whether it influences results of registration or not. Many thanks. If you need, I can send you my image data. thanks, jichao

[Dirk-Jan Kroon]

*jichao zhao
showcs3 always uses the default Matlab image ranges, thus if a double or single it sets black to 0 and white to 1, and if uint8 black to 0 and white to 255.
That your image range is 0 to 60 instead of 0 to 1 will influence registration, because the pixel based penalty will now be larger in comparison to the smoothness penalty.
I almost always normalize my data to range [0..1] by V=V-min(V(:)); V=V./max(V(:));

[jichao zhao]

many thanks for your reply.

[zhang xu]

I'm currently working on this subject, can you point some related resource or important papers to me.
Thank you!

[ff dfds]

yeha,it's a excellent work! thanks,Dirk-Jan Kroon

[James Ford]

I get a different compile error under Linux -- redefinition of a function. I think it's a typo in the code... in bspline_error_2d_double.c, you have

(line 150)
#ifdef _WIN32
  unsigned __stdcall transformvolume_error_color(double **Args) {
#else
  void transformvolume_error_color(double **Args) {
#endif

and later
(line 485)
#ifdef _WIN32
  unsigned __stdcall transformvolume_gradient_color(double **Args) {
#else
  void transformvolume_error_color(double **Args) {
#endif

Those function names should match on both?

After I change this, it compiles but doesn't link. :-(

>> compile_c_files
compiling : bspline_error_2d_double.c
image_interpolation.o: In function `mindex2':
image_interpolation.c:(.text+0x0): multiple definition of `mindex2'
bspline_error_2d_double.o:bspline_error_2d_double.c:(.text+0x0): first defined here

and 8-10 more like that.

[Dirk-Jan Kroon]

*James Ford,

Thank you, for the comments.

I've made some improvements, fixed the bugs you pointed out and tested the code in Ubuntu, new version is probably online tomorow.

Dirk-Jan Kroon

[Xi Liang]

Dear Dirk,

The volume I want to show is very noisy, and by using the function of "showcs3", the images appear to be much less noisy. I assume there are certain kind of filter used. Would you please roughly point out where it is in the code?

[Hui]

Hi, Dirk,

Thanks for sharing the excellent codes.
I have two questions about your version18e codes:

1. In image_difference.m, line 34, if(exist('type','var')==0), type='p'; end. should here type = 'd' as default? Because there is no type 'p'

2. In example_2d_affine.m, can you explain a little bit about why the scale = [1 1 1 0.1 0.1] needed here? I didn't figure out the purpose of scale here.

Thanks very much.

Hui

[Dirk-Jan Kroon]

Dear Hui,

(1) Indeed 'p' doesn't exist, I will replace it with 'd' in the next update.

(2)
The scaling [1 1 1 0.1 0.1] is because the parameters are [translateX translateY rotate resizeX resizeY]. If you change the size of the image only a little bit it will produce large pixel movements. An optimizer produces the best result if all parameters have approximately the same influence on the image transformation, thus we introduced the scaling.

Dirk-Jan Kroon

[Chao Li]

Hi Dirk,

excellent code.
One trivial thing I'd like to report is that:
in deciding max refinements steps (e.g., line 255 in 'register_image.m, Version18e'), its better to define the lower bound of MaxItt to account for small images.

[christos]

you need to put t=0 if the expression in the if does not stand

if( isempty(Mask) | (~isempty(Mask) & ~isempty(find(Mask~=0) ) ) )
 
   %your code
   if(size(U,3)==3),
   .....
   end
else
    t = 0;
end

[christos]

I ve sent two posts but the first did not show up.
RE: my post on 23 NOV 09

If one of the sub-regions in b-spline registration corresponds to a totally masked out part of the image (hence the mask is not empty but populated only by zeros), then the image_difference in the MI returns an NaN value.

As a workaround you can probably check for this case and simply set the image similarity to 0 for these regions. See previous post for the exact if statement

[Dirk-Jan Kroon]

*Christos

Thank you for your bug report, I will try to fix them in the next code update.

Best regards, Dirk-Jan

[Alex Kararg]

I am getting this error when i try to run example3d_non_rigid.m : "Undefined function or variable 'maxNumCompThreads' " I run matlab 2007a. I compiled successfully.

[Dirk-Jan Kroon]

*Alex Kararg,
The function maxNumCompThreads is a build-in function which returns the number of cores in the computer. If it doesn't exist in your Matlab you simply create the function your self like:

function a=maxNumCompThreads()
a=2;

[Alex Kararg]

Thank you for your reply. I upgraded to the newest version of Matlab and it works.

However, I would like to ask you this simple question. I have two volumes V1 (sphere) and V2 (cube). V2 is smaller than V1. I am using the "registervolumes.m" to register V2 over V1 . The output shows the final registered volumes and that is only the sphere. My question is how I could access volume V2 as it deforms to V1 on each iteration.

Thank you!

[Dirk-Jan Kroon]

Note: bug in manually_warp_images (will be fixed in next update).

Replace line 279 with:
[j,i]=ind2sub([size(data.Grid,2) size(data.Grid,1)],data.pointselected);

[winnie]

Hi Dirk,

I have just been trying to register an image to a template using your script, but unfortunately it did not work. I used all the default options. Here is a summary of the output:

Start Affine part of Non-Rigid registration
...
    Optimizer Results
        Algorithm Used: limited memory BFGS (L-BFGS)
        Exit message : Change in the objective function value was less than the specified tolerance TolFun.
... (repeat twice)

Start non-rigid b-spline grid registration
Current Grid size : 7x7
...
Optimizer Results
        Algorithm Used: limited memory BFGS (L-BFGS)
        Exit message : Change in the objective function value was less than the specified tolerance TolFun.
...
Registration Refinement
Current Grid size : 11x11
...
Optimizer Results
        Algorithm Used: limited memory BFGS (L-BFGS)
        Exit message : Line search cannot find an acceptable point along the current search
...
Registration Refinement
Current Grid size : 19x19
...
Optimizer Results
        Algorithm Used: limited memory BFGS (L-BFGS)
        Exit message : Line search cannot find an acceptable point along the current search

Then the program quits without doing any registration, as far as I can tell. I was just wondering whether you had any suggestions that I might try. For example, could I increase the specified tolerance TolFun somehow?

Thanks, Winnie

[medical imaging Ibrahim]

I have used this implementation as a step in one of my analysis systems, could you please tell me the paper that I should reference for this work ?

[Jhimli Mitra]

Hi Dirk,
   Thanks for your code. The 2D B-spline works slower without compiling the C-files and I am getting a strange error while compiling the C-files. Could you please help me with this.

compiling : bspline_error_2d_double.c
Error bspline_error_2d_double.c: 742 syntax error; found `*' expecting `;'
Error bspline_error_2d_double.c: 742 undeclared identifier `dims'
Error bspline_error_2d_double.c: 742 type error: pointer expected
Warning bspline_error_2d_double.c: 742 Statement has no effect
Warning bspline_error_2d_double.c: 745 Statement has no effect
Error bspline_error_2d_double.c: 745 syntax error; found `Osizex' expecting `;'
Error bspline_error_2d_double.c: 745 undeclared identifier `Osizex'
Error bspline_error_2d_double.c: 745 undeclared identifier `Osizey'
Warning bspline_error_2d_double.c: 745 Statement has no effect
Error bspline_error_2d_double.c: 746 illegal statement termination
Error bspline_error_2d_double.c: 746 skipping `double'
Error bspline_error_2d_double.c: 746 undeclared identifier `Osize_d'
Error bspline_error_2d_double.c: 746 type error: pointer expected
Error bspline_error_2d_double.c: 746 illegal expression
Error bspline_error_2d_double.c: 746 syntax error; found `0' expecting `;'
Warning bspline_error_2d_double.c: 746 Statement has no effect
Error bspline_error_2d_double.c: 746 syntax error; found `}' expecting `;'
Warning bspline_error_2d_double.c: 746 possible usage of Osize_d before definition
Warning bspline_error_2d_double.c: 745 possible usage of Osizey before definition
Warning bspline_error_2d_double.c: 745 possible usage of Osizex before definition
Warning bspline_error_2d_double.c: 742 possible usage of dims before definition
Warning bspline_error_2d_double.c: 742 possible usage of mwSize before definition
Warning bspline_error_2d_double.c: 734 local `pointer to double ThreadID1' is not referenced
Warning bspline_error_2d_double.c: 733 local `pointer to pointer to double ThreadID' is not referenced
Warning bspline_error_2d_double.c: 727 local `pointer to pointer to void ThreadList' is not referenced
Warning bspline_error_2d_double.c: 723 local `pointer to pointer to double ThreadArgs1' is not referenced
Warning bspline_error_2d_double.c: 722 local `pointer to pointer to pointer to double ThreadArgs' is not referenced
Warning bspline_error_2d_double.c: 719 local `int offset2' is not referenced
Warning bspline_error_2d_double.c: 719 local `int offset1' is not referenced
Warning bspline_error_2d_double.c: 715 local `int Nthreads' is not referenced
Warning bspline_error_2d_double.c: 714 local `pointer to double Nthreadsd' is not referenced
Warning bspline_error_2d_double.c: 711 local `pointer to double ThreadOut' is not referenced
Warning bspline_error_2d_double.c: 711 local `pointer to double Egradient' is not referenced
Warning bspline_error_2d_double.c: 711 local `pointer to double E' is not referenced
Warning bspline_error_2d_double.c: 711 local `pointer to double dya' is not referenced
Warning bspline_error_2d_double.c: 711 local `pointer to double dxa' is not referenced
Warning bspline_error_2d_double.c: 711 local `pointer to double I2' is not referenced
Warning bspline_error_2d_double.c: 711 local `pointer to double I1' is not referenced
Warning bspline_error_2d_double.c: 711 local `pointer to double Oy' is not referenced
Warning bspline_error_2d_double.c: 711 local `pointer to double Ox' is not referenced
Warning bspline_error_2d_double.c: 746 empty declaration
Error bspline_error_2d_double.c: 763 unrecognized declaration
Error bspline_error_2d_double.c: 763 syntax error; found `!=' expecting `)'
Error bspline_error_2d_double.c: 763 syntax error; found `!=' expecting `;'
Error bspline_error_2d_double.c: 763 skipping `!='
Error bspline_error_2d_double.c: 763 unrecognized declaration
Error bspline_error_2d_double.c: 763 unrecognized declaration
Error bspline_error_2d_double.c: 763 unrecognized declaration
Error bspline_error_2d_double.c: 764 too many errors
 
  C:\PROGRAM FILES\MATLAB71\BIN\MEX.PL: Error: Compile of 'bspline_error_2d_double.c' failed.

[Carl]

Hi Dirk

I did not find any bicubic interpolation in affine_transform_2d.m but bilinear instead.

How to use b-spline interpolation for affine_tranform?

Is there any way to calculate the affine transform matrix first rather than its inverse matrix?

Thanks,

Carl

[Dirk-Jan Kroon]

*Carl

bicubic-interpolation is available if you compile the c-code with compile_c_files. (not implemented in .m file because it would be really slow)

If you want to keep track of the forward affine matrix during registration, the best way is to put: M = inv(M); as first line in affine_transform.m and as first line in make_init_grid.m

[Tan Chye Cheah]

Hi Dirk,

When i try to run manually_warp_images, I received an error below. Thanks for your help.

Error using ==> manually_warp_images>manually_warp_images_OpeningFcn at 60
no images set

Error in ==> gui_mainfcn at 221
    feval(gui_State.gui_OpeningFcn, gui_hFigure, [], guidata(gui_hFigure), varargin{:});

Error in ==> manually_warp_images at 38
    gui_mainfcn(gui_State, varargin{:});

[Dirk-Jan Kroon]

* Tan Chye Cheah
See the example in the help op manually_warp_images.
(Note!, tomorrow a whole new version of code will be online)

[fpuspitasari puspitasari]

hi, I encountered error while compiling the C files.

compiling : bspline_error_2d_double.c
   Creating library C:\DOCUME~1\user~3\LOCALS~1\TEMP\MEX_FO~1\templib.x and object C:\DOCUME~1\user~3\LOCALS~1\TEMP\MEX_FO~1\templib.exp
image_interpolation.obj : error LNK2019: unresolved external symbol _floorf referenced in function _floorfloat
bspline_error_2d_double.mexw32 : fatal error LNK1120: 1 unresolved externals
 
  C:\PROGRA~1\MATLAB\R2008B\BIN\MEX.PL: Error: Link of 'bspline_error_2d_double.mexw32' failed.
 
??? Error using ==> mex at 213
Unable to complete successfully.

Error in ==> compile_c_files at 8
        mex(filename,'image_interpolation.c');

I wonder why this error occurs? I'm using matlab R2008b by the way

[Harrison Parks]

I can't make this code work so perhaps I'm applying it incorrectly. I want to register a smaller image to a larger one where the scale of the images differ and the smaller image will match only a portion of the larger. For example, crop the rotated Lena image and try matching it with the full unrotated image using either rigid or affine transformations. It doesn't work. In particular, the code first resizes the smaller image to the size of the larger image. A rigid transformation now can't possibly work. If I put a border around the cropped image, it still doesn't help. Suggestions will be welcomed.

[Borislav Oreshkin]

Hello Dirk,

When I tried to compile under R2007a, the maxNumCompThreads.c returned the error:

lcc preprocessor error: C:\PROGRA~1\MATLAB\R2007a\sys\lcc\include\win.h:13510 C:\PROGRA~1\MATLAB\R2007a\sys\lcc\include\windows.h:4 maxNumCompThreads.c:4 Macro redefinition of EXTERN_C

The solution for this is to exchange the order of header inclusions so that windows.h goes first:

....
#include <windows.h>
...

#include "mex.h"
#include "math.h"

May be you find this useful for the next release.

[Dirk-Jan Kroon]

*Boris Oreshkin
Thanks, I can use your solution in almost all my files on Mathworks.

[Borislav Oreshkin]

Hello Dirk,

Error in the imgaussian.m:

Hz=reshape(H,[1 1 length(H)]);
        function I=imgaussian(I,sigma,siz)

Should be
function I=imgaussian(I,sigma,siz)

[Johnson Gad Elkarim]

the 3D matlab functions are absent

[Jamie]

An incredible contribution, thank you.

[goodman]

Hi Dirk,

can you give more information about the penalty_smooth in your code?
thanks

[Carl]

Hi Dirk,

Could you please tell me how to 'b-spline' option in 'image_interpolation.m'?

Thanks,

[Chao Li]

Hi Dirk,
In function 'make_transformation_matrix', r=r*(pi/180) is used in 3D but not in 2D version. Why different representations are used in 2D/3D? Is it to compromise the scale difference between them?

Thanks.

[xia wei]

Hi Dirk,
In bspline_registration_gradient.m,I find that there is a function bspline_error_3d_double to caculate the O_error and O_grad,but it also has other two functions image_difference and error_gradient3d to caculate the O_error and O_grad, why we need the functions image_difference and error_gradient3d?And in example_3d_norigid.m,which functions give us the results of O_error and O_grad?

[ucd puri]

Dear Dirk, I want to know if this algo can register a 2d image to a 3d Volume? According to my limited understanding, such registrations are done using projection based methods.. Any suggestion on this. Many thanks.

[Xi]

Hi Dirk, very good work and thank you. I think there is a minor bug here.
In the file "affine_registration_error.m You should have if(~exist('Spacing','var')),Spacing=[]; end

Or else, when you run rigid or affine registration, you will always ave to specify the space which is only necessary when bspline registration is used.

[zhang zhen]

useful!

[Tripp]

Great submission!

I believe there is a slight bug in the 3D registration if one does non-rigid registration only. I get the following error:

??? Attempted to access M(1,1); index out of bounds because size(M)=[0,0].

Error in ==> make_init_grid at 63
        O_trans(:,:,:,1) = mean(1) + M(1,1) * xd + M(1,2) *yd + M(1,3) *zd + M(1,4)* 1;

Error in ==> image_registration>Make_Initial_Grid at 475
        O_trans=make_init_grid(Spacing,size(Imoving),M);

Error in ==> image_registration at 219
[O_trans,Spacing,MaxItt]=Make_Initial_Grid(O_trans,Spacing,Options,Imoving,M,IS3D);

If I had to guess, it seems like M (inverse transformation from rigid registration) isn't initialized in the 3D, non-rigid case. I do not get this error when running 2D non-rigid. I also do not get this error when running 3D rigid + non-rigid.

[Alberto Gomez]

@Tripp, I have found that bug too, I think it suffice to replace M=[] with M=eye(4).

There is however another bug I am afraid. I have not gone deep into the code yet but it seems that, when doing 3D NonRigid registration only 1 direction out of 3 is computed, i.e. the resulting grid only has the estimated deformation in the third component, and components 1 and 2 remain unchanged. Ani hint on how to fix this?