Thread Subject: Image Registration gradient descent

Hello,

I am trying to implement 2d image registration (translation X Y et rotation). At the moment I am working with 'fminsearch' but I would like to test the gradient descent methods. I have calculated the gradient in the following way:

g(x) = (f(x + step) - f(x - step)) / )2*step);
g(y) = (f(y + step) - f(y - step)) / )2*step);
g(Rot) = (f(Rot + step) - f(Rot - step)) / )2*step);

where f(x) is the objective function

Is it right? When I use the option with the gradient (optimset('gradObj','on')) the image move on the other image with a triangular trajectoire (like a loop) and doesn't stop.

Another problem is related with th error of the objective function. I am using - ans(corr2). Are ther others better methods?

Do you have some ideas? Thank you very much in advance for your help.

"Fabiano" <fabiano.riva@unil.ch> wrote in message <ghla0v$bar$1@fred.mathworks.com>...
> Hello,
>
> I am trying to implement 2d image registration (translation X Y et rotation). At the moment I am working with 'fminsearch' but I would like to test the gradient descent methods. I have calculated the gradient in the following way:
>
> g(x) = (f(x + step) - f(x - step)) / )2*step);
> g(y) = (f(y + step) - f(y - step)) / )2*step);
> g(Rot) = (f(Rot + step) - f(Rot - step)) / )2*step);
>
> where f(x) is the objective function

Isn't the objective funtion f(Rot)? If so, the chain rule says you have to pre-multiply
g(Rot) with the gradient of Rot.

This might be clearer if you described how you're parametrizing Rot.

"Matt" <mjacobson.removethis@xorantech.com> wrote in message <ghm4ve$dle$1@fred.mathworks.com>...
> "Fabiano" <fabiano.riva@unil.ch> wrote in message <ghla0v$bar$1@fred.mathworks.com>...
> > Hello,
> >
> > I am trying to implement 2d image registration (translation X Y et rotation). At the moment I am working with 'fminsearch' but I would like to test the gradient descent methods. I have calculated the gradient in the following way:
> >
> > g(x) = (f(x + step) - f(x - step)) / )2*step);
> > g(y) = (f(y + step) - f(y - step)) / )2*step);
> > g(Rot) = (f(Rot + step) - f(Rot - step)) / )2*step);
> >
> > where f(x) is the objective function
>
> Isn't the objective funtion f(Rot)? If so, the chain rule says you have to pre-multiply
> g(Rot) with the gradient of Rot.
>
> This might be clearer if you described how you're parametrizing Rot.
>
>

Hello,

the objective function is the following:

Res = function_obj_min(Input)

Rotation = Input(1);
X = Input(2);
Y = Input(3);

global Structure;

% The model image

A = Structure.A;

% The image that have to be registrated

B = Structure.B;

% Rotation

T_Rot = [cos(Rot) sin(Rot) 0; -sin(Rot) cos(Rot) 0; 0 0 1];

% Translation

T_Trans = [1 0 0; 0 1 0; X Y 1];

% Transformation

Transformation = T_Rot*T_Trans;

[M N] = size(B);

% Application of the transformation

Tform = maketform('affine',Transformation);
Image_Transformed = imtransform(B,Tform,'Xdata',[1 N],'Ydata',[1 M]);

Res = -abs(corr2(A,Image_Transformed));


To calculate the gradient for each variable I calculate the derivative respect at the variable and the others variables don't change. This procedure has done for each variable (Rot,Translation X and Translation Y).

Could you help me?

Thank you

Fabiano

"Fabiano" <fabiano.riva@unil.ch> wrote in message <ghtn13$1dd$1@fred.mathworks.com>...
> "Matt" <mjacobson.removethis@xorantech.com> wrote in message <ghm4ve$dle$1@fred.mathworks.com>...
> > "Fabiano" <fabiano.riva@unil.ch> wrote in message <ghla0v$bar$1@fred.mathworks.com>...
> > > Hello,
> > >
> > > I am trying to implement 2d image registration (translation X Y et rotation). At the moment I am working with 'fminsearch' but I would like to test the gradient descent methods. I have calculated the gradient in the following way:
> > >
> > > g(x) = (f(x + step) - f(x - step)) / )2*step);
> > > g(y) = (f(y + step) - f(y - step)) / )2*step);
> > > g(Rot) = (f(Rot + step) - f(Rot - step)) / )2*step);
> > >
> > > where f(x) is the objective function
> >
> > Isn't the objective funtion f(Rot)? If so, the chain rule says you have to pre-multiply
> > g(Rot) with the gradient of Rot.
> >
> > This might be clearer if you described how you're parametrizing Rot.
> >
> >
>
> Hello,
>
> the objective function is the following:
>
> Res = function_obj_min(Input)
>
> Rotation = Input(1);
> X = Input(2);
> Y = Input(3);
>
> global Structure;
>
> % The model image
>
> A = Structure.A;
>
> % The image that have to be registrated
>
> B = Structure.B;
>
> % Rotation
>
> T_Rot = [cos(Rot) sin(Rot) 0; -sin(Rot) cos(Rot) 0; 0 0 1];
>
> % Translation
>
> T_Trans = [1 0 0; 0 1 0; X Y 1];
>
> % Transformation
>
> Transformation = T_Rot*T_Trans;
>
> [M N] = size(B);
>
> % Application of the transformation
>
> Tform = maketform('affine',Transformation);
> Image_Transformed = imtransform(B,Tform,'Xdata',[1 N],'Ydata',[1 M]);
>
> Res = -abs(corr2(A,Image_Transformed));
>
>
> To calculate the gradient for each variable I calculate the derivative respect at the variable and the others variables don't change. This procedure has done for each variable (Rot,Translation X and Translation Y).
>
> Could you help me?
>
> Thank you


You should probably turn gradObj 'off' and forget about computing the gradient.

Your objective function is not differentiable. This is partly due to the use of the abs() function in

Res = -abs(corr2(A,Image_Transformed));


which will introduce non-differentiability. Furthermore, imtransform uses linear interpolation by default, so you will have lots of non-differentiable points because of that as well.

Alternatively, you could make your objective differentiable by telling imtransform to use bicubic interpolation. I think you can also get rid of the abs() and use

Res = -(corr2(A,Image_Transformed));

In the neighbourhood of the solution, the correlation will be positive anyway...

hi I am also working on the same problem. i would like to register my image with respect to translation and rotation with gradient descent. for translation part i used following function
[gy gx]=gradient(img2);
dx=sum(2*(gx(:)).*(img1(:)-img2(:)));
dy=sum(2*(gy(:)).*(img1(:)-img2(:)));
and then in the main program i update my translation matrix like this
lenght=sqrt(dx^2+dy^2);
trans(1,3)=trans(1,3)+currentStepSize*dx/lenght; trans(2,3)=trans(2,3)+currentStepSize*dy/lenght;

where currentstepsize is learning rate which is 0.5 in my case and if the answer in each iteration go far away or approach from solution from SSD then it will decrease or increase by factor of 2. but i have problem with rotation i dont how can i find gradient decent with respect to rotation and also i dont know how to update my transform matrix. i would be nice if give me some clues what did you do for rotation part

tnx

hi I am also working on the same problem. i would like to register my image with respect to translation and rotation with gradient descent. for translation part i used following function
[gy gx]=gradient(img2);
dx=sum(2*(gx(:)).*(img1(:)-img2(:)));
dy=sum(2*(gy(:)).*(img1(:)-img2(:)));
and then in the main program i update my translation matrix like this
lenght=sqrt(dx^2+dy^2);
trans(1,3)=trans(1,3)+currentStepSize*dx/lenght; trans(2,3)=trans(2,3)+currentStepSize*dy/lenght;

where currentstepsize is learning rate which is 0.5 in my case and if the answer in each iteration go far away or approach from solution from SSD then it will decrease or increase by factor of 2. but i have problem with rotation i dont how can i find gradient decent with respect to rotation and also i dont know how to update my transform matrix. i would be nice if give me some clues what did you do for rotation part

tnx

"amir yavari" <amir.yavari@yahoo.com> wrote in message <iqgi0c$fgl$1@newscl01ah.mathworks.com>...
> hi I am also working on the same problem. i would like to register my image with respect to translation and rotation with gradient descent. for translation part i used following function
> [gy gx]=gradient(img2);
> dx=sum(2*(gx(:)).*(img1(:)-img2(:)));
> dy=sum(2*(gy(:)).*(img1(:)-img2(:)));
> and then in the main program i update my translation matrix like this
> lenght=sqrt(dx^2+dy^2);
> trans(1,3)=trans(1,3)+currentStepSize*dx/lenght; trans(2,3)=trans(2,3)+currentStepSize*dy/lenght;
>
> where currentstepsize is learning rate which is 0.5 in my case and if the answer in each iteration go far away or approach from solution from SSD then it will decrease or increase by factor of 2. but i have problem with rotation i dont how can i find gradient decent with respect to rotation and also i dont know how to update my transform matrix. i would be nice if give me some clues what did you do for rotation part
=========

Is this all really necessary? Why not use

regionprops(img~=0, img, 'Orientation','WeightedCentroid')

to get the position/orientation of img1 and img2.

"Matt J" wrote in message <iqgkbc$m37$1@newscl01ah.mathworks.com>...
> "amir yavari" <amir.yavari@yahoo.com> wrote in message <iqgi0c$fgl$1@newscl01ah.mathworks.com>...
> > hi I am also working on the same problem. i would like to register my image with respect to translation and rotation with gradient descent. for translation part i used following function
> > [gy gx]=gradient(img2);
> > dx=sum(2*(gx(:)).*(img1(:)-img2(:)));
> > dy=sum(2*(gy(:)).*(img1(:)-img2(:)));
> > and then in the main program i update my translation matrix like this
> > lenght=sqrt(dx^2+dy^2);
> > trans(1,3)=trans(1,3)+currentStepSize*dx/lenght; trans(2,3)=trans(2,3)+currentStepSize*dy/lenght;
> >
> > where currentstepsize is learning rate which is 0.5 in my case and if the answer in each iteration go far away or approach from solution from SSD then it will decrease or increase by factor of 2. but i have problem with rotation i dont how can i find gradient decent with respect to rotation and also i dont know how to update my transform matrix. i would be nice if give me some clues what did you do for rotation part
> =========
>
> Is this all really necessary? Why not use
>
> regionprops(img~=0, img, 'Orientation','WeightedCentroid')
>
> to get the position/orientation of img1 and img2.