non-rigid b-spline grid image registration

Registration of two 2D/3D images, rigid and with smooth b-spline grid.

Description

Registration of two 2D color/grayscale images or 3D volumes both rigid and non-rigid with smooth b-spline grid. No landmarks are needed because registration is image region based, but landmarks can be added as transformation regularization.

This function can use normalized mutual information as registration error, allowing the images or volumes to be of a different type for instance a MRI T1 and T2 patient scan.

Description:
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).

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 functions register_images.m and register_volumes.m are easy to use, have some examples in the help, and will fit most applications. (If you want to write your own specialized registration code study the registration examples)

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.

note:
- In 2D b-spline registration approximately 30% of the registration time is used for calculating the grid smoothness penalty, setting this penalty to zero saves registration time. Also setting the option max-refinement steps will save a large amount of registration time.
- B-spline registration is slow and complex in comparison with demon registration see:
http://www.mathworks.fr/matlabcentral/fileexchange/loadFile.do?objectId=21451

Please report bugs, successes and questions.

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

Other Files

compile_c_files.m,
functions/backwards2forwards.m,
functions/backwards2forwards_2d_double.c,
functions/backwards2forwards_2d_double.m,
functions/backwards2forwards_3d_double.c,
functions/backwards2forwards_3d_single.c,
functions/fminlbfgs.m,
functions/image_difference.m,
functions/imgaussian.c,
functions/imgaussian.m,
functions/imresize3d.m,
functions/make_grid_image.m,
functions/make_init_grid.m,
functions/maple_penalties_smoothing_2d.mw,
functions/maple_penalties_smoothing_3d.mw,
functions/mutual_histogram_double.c,
functions/mutual_histogram_double.m,
functions/mutual_histogram_single.c,
functions/penalties_smoothness.m,
functions/penalty_matrix.mat,
functions/refine_grid.m,
functions/showcs3.fig,
functions/showcs3.m,
functions/squared_difference_double.c,
functions/squared_difference_double.m,
functions/squared_difference_single.c,
functions/strain.m,
functions_affine/affine_error_3d_double.c,
functions_affine/affine_error_3d_single.c,
functions_affine/affine_registration_error.m,
functions_affine/affine_registration_image.m,
functions_affine/affine_transform.m,
functions_affine/affine_transform_2d_double.c,
functions_affine/affine_transform_2d_double.m,
functions_affine/affine_transform_3d_double.c,
functions_affine/affine_transform_3d_single.c,
functions_affine/image_interpolation.c,
functions_affine/image_interpolation.h,
functions_affine/make_transformation_matrix.m,
functions_nonrigid/bspline_error_2d_double.c,
functions_nonrigid/bspline_error_3d_double.c,
functions_nonrigid/bspline_error_3d_single.c,
functions_nonrigid/bspline_registration_gradient.m,
functions_nonrigid/bspline_registration_image.m,
functions_nonrigid/bspline_transform.m,
functions_nonrigid/bspline_transform_2d_double.c,
functions_nonrigid/bspline_transform_2d_double.m,
functions_nonrigid/bspline_transform_3d_double.c,
functions_nonrigid/bspline_transform_3d_single.c,
functions_nonrigid/image_interpolation.c,
functions_nonrigid/image_interpolation.h,
functions_nonrigid/movepixels.m,
functions_nonrigid/movepixels_2d_double.c,
functions_nonrigid/movepixels_2d_double.m,
functions_nonrigid/movepixels_3d_double.c,
functions_nonrigid/movepixels_3d_single.c,
images/brain1.png,
images/brain2.png,
images/brain3.png,
images/landmarks.mat,
images/landmarks1.png,
images/landmarks2.png,
images/lena.jpg,
images/lenag1.png,
images/lenag2.png,
images/lenag3.png,
images/lenag4.png,
license.txt,
low_level_examples/addpaths.m,
low_level_examples/example_2d_affine.m,
low_level_examples/example_2d_nonrigid_1.m,
low_level_examples/example_2d_nonrigid_2.m,
low_level_examples/example_2d_nonrigid_3.m,
low_level_examples/example_3d_affine.m,
low_level_examples/example_3d_nonrigid.m,
low_level_examples/get_example_data.m,
register_images.m,
register_volumes.m

[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