Homographic image transformation distortion issue

Question

I am trying to transform an image using a 3D transformation matrix and assuming my camera is orthonormal.

I am defining my homography using the plane-induced homography formula H=R-t*n'/d (with d=Inf so H=R) as given in Hartley and Zisserman Chapter 13.

What I am confused about is when I use a rather modest rotation, the image seems to be distorting much more than I expect (I'm sure I'm not confounding radians and degrees).

What could be going wrong here?

I've attached my code and example output.

n = [0;0;-1];
d = Inf;

im = imread('cameraman.tif');

 rotations = [0 0.01 0.1 1 10];

 for ind = 1:length(rotations)
       theta = rotations(ind)*pi/180;

       R = [ 1     0           0 ;
           0  cos(theta) -sin(theta);
           0  sin(theta)  cos(theta)];

       t = [0;0;0];

       H = R-t*n'/d;

      tform = maketform('projective',H');
      imT = imtransform(im,tform);

      subplot(1,5,ind) ;
      imshow(imT)
      title(['Rot=' num2str(rotations(ind)) 'deg']);
      axis square
 end

Answer 1

The formula H = R-t*n'/d has one assumption which is not met in your case:

This formula implies that you are using pinhole camera model with focal length=1

But in your case, for your camera to be more real and for your code to work, you should set the focal length to some positive number much greater than 1. (focal length is the distance from your camera center to the image plane)

To do this you can define a calibration matrix K which handles the focal length. You just need to change your formula to H=K R inv(K) - 1/d K t n' inv(K) in which K is a 3-by-3 identity matrix whose two first elements along the diagonal are set to the focal length (e.g. f=300). The formula can be easily derived if you assume a projective camera.

Below is the corrected version of your code, in which the angles make sense.

n = [0;0;-1];
d = Inf;

im = imread('cameraman.tif');

rotations = [0 0.01 0.1 30 60];

 for ind = 1:length(rotations)
   theta = rotations(ind)*pi/180;

   R = [ 1     0           0 ;
       0  cos(theta) -sin(theta);
       0  sin(theta)  cos(theta)];

   t = [0;0;0];

  K=[300 0    0;
        0    300 0;
        0    0    1];

  H=K*R/K-1/d*K*t*n'/K;

  tform = maketform('projective',H');
  imT = imtransform(im,tform);

  subplot(1,5,ind) ;
  imshow(imT)
  title(['Rot=' num2str(rotations(ind)) 'deg']);
  axis square
end

You can see the result in the image below:

You can also rotate the image around its center. For it to happen you should set the image plane origin to the center of the image which I think is not possible with that method of matlab (maketform). You can use the method below instead.

imT=imagehomog(im,H','c');

Note that if you use this method, you'll have to change some settings in n, d, t and R to get the appropriate result. That method can be found at: https://github.com/covarep/covarep/blob/master/external/voicebox/imagehomog.m

The result of the program with imagehomog and some changes in n, d, t , and R is shown below which seems more real.

New settings are:

n = [0 0 1]';
d = 2;
t = [1 0 0]';
R = [cos(theta),  0, sin(theta);
     0,           1,          0;
     -sin(theta), 0, cos(theta)];

Answer 2

Hmm... I'm not 100% percent on this stuff, but it was an interesting question and relevant to my work, so I thought I'd play around and give it a shot.

EDIT: I tried this once using no built-ins. That was my original answer. Then I realized that you could do it your way pretty easily:

The easy answer to your question is to use the correct rotation matrix about the z-axis:

R = [cos(theta) -sin(theta)   0;
    sin(theta)  cos(theta)    0;
    0             0           1];

Here's another way to do it (my original answer):

I'm going to share what I did; hopefully this is useful to you. I only did it in 2D (though that should be easy to expand to 3D). Note that if you want to rotate the image in plane, you will need to use a different rotation matrix that you have currently coded. You need to rotate about the Z-axis.
I did not use those matlab built-ins.

I referred to http://en.wikipedia.org/wiki/Rotation_matrix for some info.

im = double(imread('cameraman.tif')); % must be double for interpn

[x y] = ndgrid(1:size(im,1), 1:size(im,2)); 

rotation = 10;
theta = rotation*pi/180; 

% calculate rotation matrix
R = [ cos(theta) -sin(theta);
     sin(theta)  cos(theta)]; % just 2D case

% calculate new positions of image indicies

tmp = R*[x(:)' ; y(:)']; % 2 by numel(im)
xi = reshape(tmp(1,:),size(x)); % new x-indicies
yi = reshape(tmp(2,:),size(y)); % new y-indicies

imrot = interpn(x,y,im,xi,yi); % interpolate from old->new indicies

imagesc(imrot); 

My own question now is: "How do you change the origin about which you are rotating the image? Clearly, I'm rotating about (0,0), the top left corner.

EDIT 2 In response to the asker's comment, I've tried again.
This time I fixed a couple of things. Now I'm using the same transformation matrix (about x) as in the original question.
I rotated about the center of the image by redoing the way i do the ndgrids (put 0,0,0) in the center of the image. I also decided to show 3 planes of the image. This was not in the original question. The middle plane is the plane of interest. To get just the middle plane, you can leave out the zero-padding and redefine the 3rd ndgrid option to be just 1 instead of -1:1.

im = double(imread('cameraman.tif')); % must be double for interpn
im = padarray(im, [0 0 1],'both');

[x y z] = ndgrid(-floor(size(im,1)/2):floor(size(im,1)/2)-1, ...
    -floor(size(im,2)/2):floor(size(im,2)/2)-1,...
    -1:1); 

rotation = 1;
theta = rotation*pi/180; 

% calculate rotation matrix
R = [ 1     0           0 ;
      0  cos(theta) -sin(theta);
      0  sin(theta)  cos(theta)];

% calculate new positions of image indicies
tmp = R*[x(:)'; y(:)'; z(:)']; % 2 by numel(im)

xi = reshape(tmp(1,:),size(x)); % new x-indicies
yi = reshape(tmp(2,:),size(y)); % new y-indicies
zi = reshape(tmp(3,:),size(z));

imrot = interpn(x,y,z,im,xi,yi,zi); % interpolate from old->new indicies

figure;
subplot(3,1,1);imagesc(imrot(:,:,1)); axis image; axis off; 
subplot(3,1,2);imagesc(imrot(:,:,2)); axis image; axis off; 
subplot(3,1,3);imagesc(imrot(:,:,3)); axis image; axis off;