Image rotation detection probelm

[Pradyumna A]

Hello,

I have 2 images. One is the reference image and the other is the rotated version of the reference image. I would like to find what is the angle by which the 2nd image is rotated wrt 1st one. Also in future, the second image will be scaled and translated along with this rotation !. This has to be done automatically (without human interference in selecting feature points etc)

Now to do this I am using Fourier Mellin transform. But this is not accurate and dosent work for all angles :(. Is there any other means to achive this task ?

Is there a Fourier Mellin transform readily implemented on Matlab ?

Thanks in advance,
Ayyala

[Matt]

"Pradyumna A" <a_pra...@yahoo.com> wrote in message <h1fn0q$aeg$1...@fred.mathworks.com>...

> Hello,
>
> I have 2 images. One is the reference image and the other is the rotated version of the reference image. I would like to find what is the angle by which the 2nd image is rotated wrt 1st one. Also in future, the second image will be scaled

----------

In what way will the image be scaled? Do you mean stretched geometrically or do you mean all the pixel values will be scaled?

------

>and translated along with this rotation !. This has to be done automatically (without human interference in selecting feature points etc)

-----

One thing you can try is to compute

P1=abs(fft(radon(Image1)))
P2=abs(fft(radon(Image2)))

You will see that P1 is shifted with respect to P2 by an amount corresponding to the rotation. You can use your favorite shift measurement technique (e.g. correlation) to calculate the shift, and hence the rotation. Also, translation will not affect P1 and P2.

[Yuan-Liang Tang]

Matt:
The code you provided is really good. What I did is to further sum the columns of P1 and P2 and perform cross-correlation to find the angle of rotation. To be more specific,

%%%%%%%%%%%%%%%%
P1=abs(fft(radon(Image1)));
P2=abs(fft(radon(Image2)));
a = sum(P1);
b = sum(P2);
a = a-mean(a);
b = b-mean(b);
[m n] = max(xcorr(a, b));
angle = length(a)-n+1 % The angle of rotation
%%%%%%%%%%%%%%%%

The code above worked just fine. However, what is the purpose of "fft" after radon transform? Can we skip the "fft" part? Thanks.

Y.L. Tang

"Matt " <x...@whatever.com> wrote in message <h1g8gu$mrc$1...@fred.mathworks.com>...

[Phillip]

I stumbled across this while trying to solve the same problem (determine rotation between 2 images). Thanks to those who posted here. You should only need the radon projection along one angle for one of the images. Also, the use of the FFT confused me so I have used it at the end to do the correlation in a more familiar way to me. Code snippet:

thetas = [-45 : 45]; % Enter desired range

P1 = radon(Image1,0); % Only 1 projection is needed for reference image
P2 = radon(Image2,thetas); % All projections to be considered for the other image

% Compare the projections for Image2 to the projection for Image1 using cross
% correlation. The most correlated projection gives the amount of image rotation.

for ii = 1:length(thetas)
b = P2(:,ii);
this_correlation = abs(ifft(P1.*conj(b)));
angle_scores(ii) = max(max(this_correlation));
end

"Yuan-Liang Tang" <john.do...@example.com> wrote in message <h8q5ba$fds$1...@fred.mathworks.com>...

[blue ice]

"Phillip" wrote in message <i37oqi$p65$1...@fred.mathworks.com>...

> I stumbled across this while trying to solve the same problem (determine rotation between 2 images). Thanks to those who posted here. You should only need the radon projection along one angle for one of the images. Also, the use of the FFT confused me so I have used it at the end to do the correlation in a more familiar way to me. Code snippet:
>
> thetas = [-45 : 45]; % Enter desired range
>
> P1 = radon(Image1,0); % Only 1 projection is needed for reference image
> P2 = radon(Image2,thetas); % All projections to be considered for the other image
>
> % Compare the projections for Image2 to the projection for Image1 using cross
> % correlation. The most correlated projection gives the amount of image rotation.
>

while trying this code i got this error....can you tell me WHY ?????

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

[Phillip]

> while trying this code i got this error....can you tell me WHY ?????
>
> ??? Error using ==> times
> Matrix dimensions must agree.
>
> > for ii = 1:length(thetas)
> > b = P2(:,ii);
> > this_correlation = abs(ifft(P1.*conj(b)));
> > angle_scores(ii) = max(max(this_correlation));
> > end

Make sure that P1 is a column vector and P2 has the same number of rows as P1.

However, I've changed the way I solve this problem now. Two points to note for anybody stumbling on this:

1) What I said about only needing one radon projection for the reference might be true in theory, but I found in practice with noisy data that this was not so good. I've actually abandoned the radon approach.

2) I now solve my rotation problem in an entirely different way. I'm not sure how generally applicable this solution is, but it worked well for my application (very large images, only small angles of rotation possible, very precise estimate of rotation necessary).

I find the polar transforms of the reference and acquired images, and then do a typical registration between these using cross-correlation. The shift required in the theta dimension tells you the rotation required (the other dimension has radius info which I don't use).
I got a Matlab function for creating a polar transform from this page:

http://www.csse.uwa.edu.au/~pk/Research/MatlabFns/

This function lets you specify the sampling in both the theta and radial dimensions. I upsample quite a bit in the theta dimension (better precision on the angle estimate) and downsample in the radial dimension (to keep memory costs down - my images are quite large).

Hope this approach is useful to others.

[Matt J]

"Phillip" wrote in message <ierba4$kko$1...@fred.mathworks.com>...

>
>
> 1) What I said about only needing one radon projection for the reference might be true in theory, but I found in practice with noisy data that this was not so good. I've actually abandoned the radon approach.

=====

Why not just increase the number of projection angles used for the reference image, so as to decrease the SNR?

[Phillip]

"Matt J" wrote in message <iercbt$rf9$1...@fred.mathworks.com>...

I did try that when I started having issues. I didn't have much luck with it; I'm not sure why but it's probably because the angles of rotation are so small in my data. I could probably upsample the images to make this work better, but my images are quite large and I frequently have memory issues in Matlab. I liked the polar transform approach better because I was able to upsample in just the one dimension.

[Darren g]

Hi Phillip,
I am also trying to measure the image rotation for very small angle (fractional as well) noisy images. Could you please let me know the function you used from the link you provided. Also, how did you calculate the angle of rotation from the polar transform. I have used polar transform function from the one got from the internet but do not know how to proceed for the angular measurement after polar transforming the two images.
Any help would be much appreciated.
Thanks.
Darren

"Phillip" wrote in message <ierba4$kko$1...@fred.mathworks.com>...

[Phillip]

Hi Darren,

From the link I gave above, I used the "polartrans.m" file.

Once you have the two polar images, do a plain X/Y (translation) cross-correlation. If the images are just rotated and not stretched, the peak of the cross-correlation function will be displaced by some amount (I forget in which direction, but should be obvious when you display the images). Each division in the image is worth 2*pi/ntheta angular units (where ntheta is an argument that you give to the polartrans function).

If you haven't done a cross-correlation before, have a look at xcorr2 (I think this requires the signal processing toolbox). Play around with images that just require a translation - once you have that figured out, do the same thing but using your polar images, and multiply by the amount given above to determine the angle.

NB: In R2012 onwards the Image Processing toolbox has an "imregister" function that can apparently do all sorts of registrations, including rotations. I find it works well for some kinds of images but not others; you might want to give it a shot.

"Darren g" wrote in message <k806kl$s5m$1...@newscl01ah.mathworks.com>...

[Darren g]

Hi Phillip,

Thanks for the detailed information. I am not an expert in matlab but currently learning for this image rotation measurement. I have tried the polartrans.m function. But, when calling this function what values I need to supply within the code. I tried nrad=256 and ntheta=360, are the width of my image and a full range of theta. My understanding is that the other parameters cx, cy, linlog, shape should take some default values but when running the code it show error "??? Undefined function or variable 'cx'". Here is the code I tried for polar transform.

a=1;
image = strcat('',num2str(a,'%04d'),'.png');
im = imread(image);
figure(1);
imagesc(im);
axis image;
nrad = 256; ntheta = 360;
pim = polartrans(im, nrad, ntheta, cx, cy, linlog, shape);
figure(2);
imagesc(pim);title('Polar transform of image');
axis image;

Any error in this code?

Thanks in advance.
Darren


"Phillip" wrote in message <k81qc1$app$1...@newscl01ah.mathworks.com>...

[Phillip]

Hi Darren,

When you have optional arguments like that, you have to omit them if you don't want to use them. Otherwise the function thinks you are trying to specify them, but because you haven't defined cx, cy, linlog, shape it gives an error. i.e. you need to change the call to polartrans to:

pim = polartrans(im, nrad, ntheta);


"Darren g" wrote in message <k8duk1$kjn$1...@newscl01ah.mathworks.com>...

[Darren g]

Hi Phillip,

Many thanks. Yes. It is now working with your suggession. The problem I have encountered with those optional arguments is that I tried without quote as below
pim = polartrans(im, nrad, ntheta,cx,cy, log, valid);

Now it is working with quote for the 'log' and 'valid' as below.
pim = polartrans(im, nrad, ntheta,cx,cy, 'log', 'valid');
So far got to work up to transforming the two images to polatranform and can clearly see a shift on the two transformed images.
I have not understand the xcorre2 for cross-correlation that you mentioned earlier and the multiplication to calculate the angle of rotation and would appreciate any suggestion or part of the code if you can remember now.

Thanks.
Darren


"Phillip" wrote in message <k8em3u$eut$1...@newscl01ah.mathworks.com>...

[Phillip]

Hi Darren,

Try something like

correlationMatrix = xcorr2(image1,image2);
figure; imagesc(correlationMatrix)

You will see a bright patch on a dark background. If the 2 input images were the same, the patch would be centered in the array. If the 2 input images are displaced relative to each other, the bright patch will be displaced by that amount. To quantify the amount, you can start with:

[~,maxColumn] = max(max(correlationMatrix));
[~,maxRow] = max(correlationMatrix(:,maxColumn));

This will tell you the column and row that the brightest point of the image was found in. You can then compare this to the (half) size of the whole array to quantify the displacement.

If you want a more precise (sub-pixel) estimate of the displacement, you can look into calculating the centre of mass of the bright patch instead of just finding the peak.

"Darren g" wrote in message <k8h0dd$4uj$1...@newscl01ah.mathworks.com>...

[Darren g]

Hi Phillip,

Many thanks. Yes. My code is now working after adopting your suggestion for finding the column and row in which the maximum peak intensity is found. For finding the centre of mass of the image to measure the angle of rotation smaller than one degree I am trying the regionprobs but not yet found a way on how to integrate this within my code. Is this the one you mentioned for finding to centre of mass or any other easy way of doing? Your help for this is much appreciated.

Darren

"Phillip" wrote in message <k8hjkn$69s$1...@newscl01ah.mathworks.com>...

[Phillip]

Hi Darren,

Below is a basic function for finding centre of mass, I can not recall where I stole it from. Try to find the centre of mass in a small region of interest around the peak (instead of using the whole array). Just remember to convert back to whole array co-ordinates to give the final displacement.

Phil

function [a,b] = find_com(x,y,A)

% Finds centre of mass of matrix A in x, y index notation

if nargin == 1
A = x;
[Ar Ac] = size(A);
x = 1:Ac;
y = (1:Ar)';
else
y = y(:);
x = (x(:))';
end

sumA = sum(A(:));
if sumA == 0
a = 0;
b = 0;
else
[row col] = size(A);
% Find centroid across columns:
temp = A.*(ones(row,1)*x);
a = sum(temp(:))/sumA;
% Across rows:
temp = A.*(y*ones(1,col));
b = sum(temp(:))/sumA;
end

"Darren g" wrote in message <k8v2ok$lam$1...@newscl01ah.mathworks.com>...
> Hi Phillip,
>

[Darren g]

Hi Phillip,

Thank you for the function details for finding the centre of mass. I have used this function within my code to find the centre of mass of the intensity peak region. However, the angular rotation is not resulting correctly but displayed with four decimal places. Below is the part of the code. The maximum intensity is found at column 361 and row 141 for the first two set of images and a small change for the rest.

[max_cc_Matrix, maxColumn] = max(max(cc_Matrix))
[max_cc_Matrix, maxRow] = max(cc_Matrix(:, maxColumn))
%Angle_of_rotation = maxColumn-360 % for finding rotation using the intensity peak.
cc1_Matrix = cc_Matrix(130:150,340:380); % only the region surrounding the
% maximum peak is considered.
[a,b] = find_com(cc1_Matrix) % calling function for finding centre of mass
Angle_of_rotation = 340+a-360
A1(k,:) = [k, Angle_of_rotation];
disp(A)

This will give the image no and the angle of rotation of a set of images.
Any error in calling the function or calculation?

Thanks.
Darren


"Phillip" wrote in message <k91fgr$e3d$1...@newscl01ah.mathworks.com>...

[Phillip]

Hi Darren,

Is it a coincidence that the peak of your cross-correlation appears to be at ~360, or are you getting that confused with the angular expanse of the polar image? What is the array size?

To avoid confusion you should not "hard code" these numbers. eg for one of your lines I would replace:

cc1_Matrix = cc_Matrix(130:150,340:380);

With something like:
cc1_Matrix = cc_Matrix(rowCentre - roiDiam/2 : rowCentre + roiDiam/2, colCentre - roiDiam/2 : colCentre + colDiam/2);

You also have to calculate the angle based on the ntheta input to the polartrans function. For example, if you specified ntheta = 720, each position in the array would be worth 0.5 deg.

So instead of:
Angle_of_rotation = 340+a-360

Try (I think):

Angle_of_rotation = (a + colCentre - xSize/2) * 360/ntheta

Where xSize is the diameter of the whole array.

Make sure you get an angle of zero for the test case where the images are the same.


"Darren g" wrote in message <k93t79$7d3$1...@newscl01ah.mathworks.com>...

[Darren g]

Hi Phillip,

I am completely messed up with the calculation now. Shall we forget about the centre of mass calculation for a moment and find the angle of rotation for one degree precision.

The test images I used were 283x290 pixels size and one was deliberately rotated with respect to other by 1.4 degrees.

I used the polartrans function with the following values for these images nrad = 128;

ntheta = 360; cx = (283/2)+0.5; cy = (290/2)+0.5;
pim1 = polartrans(image1, nrad, ntheta, cx, cy, 'linear', 'valid');
pim2 = polartrans(image2, nrad, ntheta, cx, cy, 'linear', 'valid');

and
cc_Matrix = xcorr2(pim1, pim2);

Now the size(cc_Matrix) gave
ans=
255 719
When I use the following for calculating the maximum intensity and its corresponding column and row it gave

maxColumn=361 and
maxRow=128

Since there is no translation and the second image is only rotated with respect to the first (by 1.4 degrees in this case), for the angular rotation calculation we only need the column value in which the maximum intensity is found which is 361 (Is it correct?).

[max_cc_Matrix, maxColumn] = max(max(cc_Matrix))
[max_cc_Matrix, maxRow] = max(cc_Matrix(:, maxColumn))

So the Angle of rotation = maxColumn-centre of the cc_matrix column
=361-360
Which resulted in the 1 degree.
Is anything wrong with this calculation?

Now move on to the centre of mass calculation.

Would you mean maxColumn and maxRow by rowCentre and colCentre in the following centre of mass (COM) calculation? Is roiDiam is the width of the square area considered for this COM calculation?

cc1_Matrix = cc_Matrix(rowCentre - roiDiam/2 : rowCentre + roiDiam/2, colCentre - roiDiam/2 : colCentre + roiDiam/2);

Is xSize is the diameter of the cc_Matrix in the following?

Angle_of_rotation = (a + colCentre - xSize/2) * 360/ntheta

Sorry for these lengthy explanation and clarification but I hope this will give an idea to find out the error in my thinking.

Thanks.
Darren

"Phillip" wrote in message <k963fc$fb5$1...@newscl01ah.mathworks.com>...

[Phillip]

Hi Darren,

In regards to the first part, the mid point for a 255 x 719 matrix is 127, 359. So that would imply some sort of error. Try subtracting the mean of both images just before the call to xcorr2, I think this is necessary for the correlation to work. i.e.

pim1 = pim1 - mean(mean(pim1));
pim2 = pim2 - mean(mean(pim2));

You might want to test your correlation routine out on a simple translation first. eg the code below will give you a 20 pixel vertical translation

xShift = 0;
yShift = 20;
translatingElement = translate(strel(1), [yShift xShift]);
image2 = imdilate(image1,translatingElement);

image1 = image1(1+yShift:end-yShift,1+xShift:end-xShift); % get rid of Infs from the imdilate function
image2 = image2(1+yShift:end-yShift,1+xShift:end-xShift);

For your other questions:

> Would you mean maxColumn and maxRow by rowCentre and colCentre in the following centre of mass (COM) calculation? Is roiDiam is the width of the square area considered for this COM calculation?

Yes.

> Is xSize is the diameter of the cc_Matrix in the following?
>
> Angle_of_rotation = (a + colCentre - xSize/2) * 360/ntheta

Yes.

Hope this helps.

Phillip

"Darren g" wrote in message <k99436$9u7$1...@newscl01ah.mathworks.com>...

[Darren g]

Hi Phillip,

I think the problem is in the Angle of rotation calculation equation but I could not figure out how to change it.
Here is the code I used to check with the code you provided to introducing the image shift.

clear all;
image1 = imread('0001.png');

xShift = 0;
yShift = 20;
translatingElement = translate(strel(1), [yShift xShift]);
image2 = imdilate(image1,translatingElement);

image1 = image1(1+yShift:end-yShift,1+xShift:end-xShift);
% get rid of Infs from the imdilate function
image2 = image2(1+yShift:end-yShift,1+xShift:end-xShift);

nrad = 128; ntheta = 360; cx = (283/2)+0.5; cy = (290/2)+0.5;

% cx & cy for a 283x290 pixels image

pim1 = polartrans(image1, nrad, ntheta, cx, cy, 'linear', 'valid');

pim2 = polartrans(image2, nrad, ntheta, cx, cy, 'linear', 'valid');

pim1 = pim1 - mean(mean(pim1));
pim2 = pim2 - mean(mean(pim2));

cc_Matrix = xcorr2(pim1, pim2);
figure (1);
imagesc(cc_Matrix);
axis image;

[max_cc_Matrix, maxColumn] = max(max(cc_Matrix))
[max_cc_Matrix, maxRow] = max(cc_Matrix(:, maxColumn))

roiDiam = 20;
cc1_Matrix = cc_Matrix(maxRow-roiDiam/2:maxRow+roiDiam/2, maxColumn-roiDiam/2:maxColumn+roiDiam/2);
figure (2);
imagesc(cc1_Matrix);
axis image;

[a,b] = find_com(cc1_Matrix) % calling function

Angle_of_rotation = (a + maxColumn - length(cc1_Matrix)/2) * 360/ntheta

Could you pls have a look of this with any of your random image?

Thanks.
Darren


"Phillip" wrote in message <k9gnhk$d5h$1...@newscl01ah.mathworks.com>...

[Phillip]

Hi Darren,

That code to do the image shift gives a pure translation. But then you are trying to take the polar transform, and determine an angle of rotation (!). So I would remove that polar image stuff and just try to extract the translation amount (i.e. you should caluclate a displacement of 20 pixels in the vertical direction).

Once you have that down pat, you can switch to rotating the image, doing the polar transform, and doing the xcorr of that. To rotate the image you can do something like

rotatedImage = imrotate(image1,angleInDegrees,'crop','bilinear');

Then finally, you can test your code by both translating and shifting the image. To solve this I usually do an initial translation correction (Cartesian images), followed by a rotation correction (polar images), then a final translation correction (Cartesian images) to mop up any residual errors.

Cheers,

Phillip

"Darren g" wrote in message <k9jeil$gc6$1...@newscl01ah.mathworks.com>...

[Darren g]

Hi Phillip,

Thanks. I have tried the following code for measuring the pure shift but it show the error

Warning: CONV2 on values of class UINT8 is obsolete.
Use CONV2(DOUBLE(A),DOUBLE(B)) or CONV2(SINGLE(A),SINGLE(B)) instead.
> In uint8.conv2 at 11
In xcorr2 at 17
In translation at 18

The code is here.

clear all;
image1 = imread('0001.jpg');

xShift = 0;
yShift = 20;
translatingElement = translate(strel(1), [yShift xShift]);
image2 = imdilate(image1,translatingElement);
image1 = image1(1+yShift:end-yShift,1+xShift:end-xShift);
% get rid of Infs from the imdilate function
image2 = image2(1+yShift:end-yShift,1+xShift:end-xShift);

figure (1);
imagesc(image1);
axis image;
figure (2);
imagesc(image2);
axis image;

d = xcorr2(image1,image2);
figure (3);
imagesc(d);
axis image;
[maxd, maxd_pos_Col] = max(max(d))
[maxd, maxd_pos_Row] = max(d(:, maxd_pos_Col))

Shift = RowCentre - maxd_pos_Row

I could not figure out the error?

Thanks.
Darren


"Phillip" wrote in message <k9k408$p5c$1...@newscl01ah.mathworks.com>...

[Phillip]

Hi Darren,

The warning appears because you can't do much useful maths on an integer array. Just follow the suggestion shown in the warning. This will convert the image from integer to double/single precision "float" (i.e. a number that can have decimal points). Single uses half the memory but you get half the precision in terms of decimal points etc... generally use single if you are worried about memory issues, otherwise, you may as well use double. The default for new Matlab variables is double.

Cheers,

Phillip

"Darren g" wrote in message <ka74i0$32b$1...@newscl01ah.mathworks.com>...

[Darren g]

Hi Phillip,

Would you mean to change the conv2 function within the corresponding folder.
I tried replacing it with double but now show another error.

Error: File: conv2.m Line: 1 Column: 26
Unbalanced or unexpected parenthesis or bracket.

Error in xcorr2 (line 17)
c = conv2(a, rot90(conj(b),2));

Error in translation (line 18)
d = xcorr2(image1,image2);

Thanks.
Darren


"Phillip" wrote in message <ka8dih$8d6$1...@newscl01ah.mathworks.com>...

[Steven_Lord]

"Darren g" <darre...@yahoo.com> wrote in message
news:kaaidk$6jc$1...@newscl01ah.mathworks.com...

> Hi Phillip,
>
> Would you mean to change the conv2 function within the corresponding
> folder.
> I tried replacing it with double but now show another error.
>
> Error: File: conv2.m Line: 1 Column: 26
> Unbalanced or unexpected parenthesis or bracket.
>
> Error in xcorr2 (line 17)
> c = conv2(a, rot90(conj(b),2));
>
> Error in translation (line 18)
> d = xcorr2(image1,image2);

Revert your changes to CONV2 then change your translation function or script
to pass double(image1) and double(image2) into XCORR2 rather than image1 and
image2.

--
Steve Lord
sl...@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

[Kev]

This solution is super useful for my application. I am basically working on an feature tracking/pattern recognition application and my method has to be rotation and intensity invariant. So one question is that will I get a higher xcorr peak if the 2 images are matching compared to none matching ones?

Cheers,
K

"Phillip" wrote in message <ierba4$kko$1...@fred.mathworks.com>...

[Phillip]

It's true that for a given pair of images the cross correlation peak will be higher when they are registered correctly. But it is possible for one pair to have a higher peak than another pair even if the latter is better aligned - it all depends on the images. In other words, the value at the peak is not directly comparable across image pairs (unless the image pairs all have quite similar spatial content).

"Kev " <piak...@gmail.com> wrote in message <kau56h$9p7$1...@newscl01ah.mathworks.com>...

[Kev]

Hi Phillip,

Thanks for replying. I do not get the part on correct registration. Say if I have an image of the letter 'A' rotated. I try to do a xcorr of the polar transform of this image with the image of the letter 'A' upright. The peak of the xcorr will be higher than xcorr of a letter 'B' with 'A' right?

K
"Phillip" wrote in message <kb03n6$n63$1...@newscl01ah.mathworks.com>...

[Kev]

Hi Phillip,

Thanks for replying. I do not get the part on correct registration. Say if I have an image of the letter 'A' rotated. I try to do a xcorr of the polar transform of this image with the image of the letter 'A' upright. The peak of the xcorr will be higher than xcorr of a letter 'B' with 'A' right?

K
"Phillip" wrote in message <kb03n6$n63$1...@newscl01ah.mathworks.com>...

[Darren g]

Dear Steve,

Thanks for your suggestion for changing to double. It took some time for me to work on the code again for some other reasons. Come back if any help needed.

Thanks.
Darren

"Steven_Lord" <sl...@mathworks.com> wrote in message <kaaj9e$9s0$1...@newscl01ah.mathworks.com>...