Log-polar image sampling

Resamples an image from a conventional grid to a log-polar grid, and back.

Description

Images sampled on a log-polar grid have interesting and useful properties. Rotation and scaling are trivial operations, and this can lead to efficient algorithms for straight line detection and optic flow estimation. Such sampling schemes also provide a rough model of biological foveal vision, of interest in the development of active computer vision systems.

The functions in this file carry out resampling from a conventional image to a log-polar image, and back. There is usually considerable information loss in each direction, but resampling to log-polar can still be useful for computational experiments.

The toolbox function imtransform does the main work; these functions are a front end to it.

Acknowledgements

This file inspired Affine Optic Flow.

• Required Products: Image Processing Toolbox
• MATLAB release: MATLAB 7.10 (R2010a)

[Gomathy]

Hi Dave, a quick update. The log-sample routine works fine for phase correlation too. Thanks,G

[David Young]

Hi Gomathy, I'm sorry, but I don't know very much about phase correlation, and not enough to help with your code I'm afraid.

[Gomathy]

Hi David, the correct LP parameters for the earlier post is
rmin = 5;
rmax = 125;
nw = 500;

I tried various ways to measure theta using Phase Correlation. It doesn't work. Can you please look at the test code ?

[Gomathy]

Hi David, This is Gomathy again. I am trying to use Phase correlation to measure displacements after converting the image into LP. It fails to detect displacement, whereas the same LP converted image works fine for normxcorr2. I am not sure, where am i making mistakes. Here is the test code. I would greatly appreciate, if you can you can please help.

Iref = imread('cameraman.tif');
deg = 73;
I = imrotate(Iref, deg, 'bicubic', 'crop');

rmin = 5;
rmax = 300;
nw = 1500;
xc = size(Iref, 2)/2;
yc = size(Iref, 1)/2;

lp_Iref = logsample(Iref, rmin, rmax, xc, yc, [], nw);
lp_I = logsample(I, rmin, rmax, xc, yc, [], nw);

THETA_F1 =fftshift(fft2(lp_Iref));
THETA_F2 =fftshift(fft2(lp_I));
a1 = angle(THETA_F1); a2 = angle(THETA_F2);
THETA_CROSS = exp(-1* (a1 - a2)); % Phase Difference
THETA_PHASE = abs(fftshift(ifft2(THETA_CROSS)));
[peak, idx] = max(THETA_PHASE (:));
[ypeak,xpeak] = ind2sub(size(THETA_PHASE ),idx);
theta=xpeak-(size(lp_I,2)+1)/2;

[Michael Meyers]

Hello David

How would your code be modified to include calculating the scale from the xpeak? I can't find any examples of that calculation on the web. I understand the angle as the y axis covers 2pi but what does the xaxis cover in terms of scale?

[Bharath Ramesh]

Hello David,

Do you know how to find the image point corresponding a point in the log-polar transformed space? Is there an easy way to do it?

[Janaka]

Hi David
Thanks for your submission. I have been following your answers and submissions and they all are very useful. Greatly appreciate your work.

I am trying to use imrotate for checking the rotation between two image frames (1024x1024 int16) coming from HH and HV polarized antennas of a marine radar. They are 16-bit data (comes as Matlab matrices) and the HV is rotated w.r.t HH. I tried to use the above piece of code but first I tried with the rotation of HH w.r.t. itself but drot always comes as 0. Is it due to the data? Should this valid for int8 only? Here is the modified code I used. Greatly appreciate your reply.

HHdble=double(HH);
HVdble=double(HV);

Iref = HHdble;
deg = 10;

I = imrotate(Iref, deg, 'bilinear', 'crop');
rmin = 5;
rmax = 512; %%125;
nw = 512;
xc = size(Iref, 2)/2;
yc = size(Iref, 1)/2;
lp_Iref = logsample(Iref, rmin, rmax, xc, yc, [], nw);
lp_I = logsample(I, rmin, rmax, xc, yc, [], nw);
cc = normxcorr2(lp_I, lp_Iref);
[max_cc, imax] = max(cc(:));
[ypeak, xpeak] = ind2sub(size(cc), imax)
drot = 360 * (ypeak - nw) / nw %% comes as '0' .

[Gomathy]

Thanks David.

[David Young]

Hi Gomathy. Measuring small rotations is difficult, but you could try my affine optic flow estimator at http://www.mathworks.co.uk/matlabcentral/fileexchange/27093-affine-optic-flow If you have difficulty with it, please let me know.

[Gomathy]

Hi David, Thanks for your reply and explaining about estimation of deg from ypeak. The medical image sequences, what i am dealing with, have rotations less than 1degree. Even, 0.05,.. Though it is a small value, it turn up as a swing in registered video seq. I tried couple of other methods. But, Nothing works for less than a degree. Now, Log-polar also fails. Do you have any other recommendation for < 1deg? Thanks agian.

[David Young]

Hi Gomathy. Your code is almost OK. One problem is that you have chosen very small rotations to test, so it is hard to see if the results are right. Instead of fractions of a degree, try using many degrees, so that you can tell if there's a big error or just a small inaccuracy.

The second problem is that your conversion from the y-coordinate of the correlation peak to an angle is not correct. Remember that the number of wedges is just a full circle, so in your case a change of 500 in the y-coordinate corresponds to 360 degrees or 2*pi radians. If you use this, then the code works.

Here is my test code, based on yours, but simplified to make it easier to analyse.

Iref = imread('cameraman.tif');
deg = 73;
I = imrotate(Iref, deg, 'bilinear', 'crop');

rmin = 5;
rmax = 125;
nw = 500;
xc = size(Iref, 2)/2;
yc = size(Iref, 1)/2;

lp_Iref = logsample(Iref, rmin, rmax, xc, yc, [], nw);
lp_I = logsample(I, rmin, rmax, xc, yc, [], nw);

cc = normxcorr2(lp_I, lp_Iref);
[max_cc, imax] = max(cc(:));
[ypeak, xpeak] = ind2sub(size(cc), imax)

drot = 360 * (ypeak - nw) / nw % prints 72.72 for initial deg=73

The choice of the parameters is up to you - you seem to have chosen good ones for this problem.

[Gomathy]

Hi, i am using Log-Polar for measuring and correcting the image rotation. For test purpose, cameraman.tif is rotated to known value. And using log-polar , we like to measure the rotation and de-rotate to end. It is too difficult to find the suitable parameters for rmin/rmax/nr/nw. Below code doesn't measure the rotation. Can you please help to fix this issue.
***************
I(:,:,1)=imread('cameraman.tif');
deg=[0 0.1,-0.2,0.3,-0.4,0.5,-0.6,0.7,-0.8,0.9];
rad=deg*(pi/180);
for i=2:10 I(:,:,i)=imrotate(I(:,:,1),-deg(i),'bilinear','crop');
end
[h w]=size(I(:,:,1));
lp_Iref = logsample(I(:,:,1), 5, 125, w/2, h/2, [], 500); %mri
for i=2:10
clear lp_I ;
lp_I = logsample(I(:,:,i),5, 125, w/2, h/2, [], 500);
clear cc;
cc = normxcorr2(lp_I,lp_Iref);
[max_cc, imax] = max((cc(:)));
[ypeak(i),xpeak(i)] = ind2sub(size(cc), imax(1));

xpeak(i)=xpeak(i)+xx1;
ypeak(i)=ypeak(i)+yy1;
drot_r(i)=size(lp_Iref,1)-ypeak(i);
dscale(i)=size(lp_Iref,2)-xpeak(i);
end
drot_d=drot_r*(180/pi);

*************
The measured rotation in drot_d is not same as deg. It would be great if you can help to fix this.

Regards

[David Young]

Yeoh cs: There are no general rules for setting those parameters. It depends on the original image and the purpose of doing the transform.

rmax determines the "field of view" of the transform. That is, it gives the radius of the region of the original image that is sampled.

nw determines the resolution at the outside ring. If you want to capture the full detail of the original image, with one transform pixel for each original pixel, you need to make nw equal to 2*pi*rmax. However, this results in massive oversampling nearer the centre, so you may want a smaller value.

I normally set only one of rmin and nr, and allow the other to default, so that the pixels in the sampled image have an aspect ratio of 1. There is always a "hole" at the centre of a log-polar image, and rmin determines how big this hole is. Near the centre of the log-polar image the original image is grossly oversampled, and the smaller rmin is, the worse this is.

In the end though, you have to look at what you are using the log-polar image for, and use either theoretical arguments or empirical tests to decide the parameters.

[Yeoh cs]

Young: thx for ur advice, 1 more question .... how to defined the rmin /rmax /nr /nw value of different pattern ?

[David Young]

Yeoh cs: I'm sorry, but I have not written a tutorial on the use of this transform, and there isn't room to write one here. There are many papers on applications of log-polar transforms, and it would be best to start with those and then ask specific questions on MATLAB Answers.

[Yeoh cs]

can this code use to do the recognition? how do i compare the 2 image using function of logsample? i m new to log polar transform. i wish to use log polar trans to do the object recognition.
can u giv me an example that how to compare 2 logimage after going through the logsample function.. Thanks

[David Young]

Masthan Zayid: Sorry, I don't know what the log polar gabor transform is.

[masthan zayid]

hi...i need code for 2D log polar gabor transform for my project....

[David Young]

Aravind: The words and equations in the documentation give the meanings of the arguments as clearly and precisely as I can - but it's true there isn't a tutorial there on the log-polar transform. Perhaps you need to look at the literature first to make sure you're familiar with the main ideas - perhaps start with Weiman & Chaikin's 1979 paper, which Google will find for you.

Here's an example of applying the function - you can change parameters to see what they do:

im = imread('cameraman.tif');
imshow(im);
logim = logsample(im, 5, 200, 150, 150, [], 300);
figure; imshow(logim);

There is a much more complex example in the code for my FEX contribution 27093, but unless you're already familiar with the basic ideas this is unlikely to shed a lot of light.

[Aravind Natarajan]

Since I'm new to log polar transform I'm not able to understand the arguments for the function logsample.Can you please post an example of how to use the function.Thank you.

[David Young]

Thanks Vlatko. Yes, the tilde as a marker for an usused result only came in with release 2009b. However, if you just omit it, the function will no longer work properly for RGB images. A better change is to use
[U, V, unused] = size(arr);
which is compatible with older versions and preserves the functionality.

[Vlatko Becanovic]

At line 65 in the function logsample.m it says:
[U, V, ~] = size(arr);
If I change this to be:
[U, V] = size(arr);
it all works great.
Thanks for a nice program!