I imrotate() an image, draw two lines, rotate back the lines and draw them in orginal image, but don't get the expected result in MATLAB?

Question

What I want to do :

Suppose I have an image I that I rotate -45° using imrotate (I get I_R). Then I draw two lines AB and CD (parallels). Finally, I rotate the two lines back (45°) and draw them in the original image I.

How I do that ##

I rotate I using the MATLAB function imrotate():

I_R = imrotate(I,-45);

From Matlab help, I get : B = imrotate(A,angle) rotates image A by angle degrees in a counterclockwise direction around its center point.

But it seems that imrotate add a translation to the image ! I have read the code of built-in matlab function, it seems that it uses a function called getOutputBound to check if the rotated image will fit in the figure. This translation is what I am looking for !!

The four points A,B,C,D form a two parallels lines AB & CD.

A = [x_A; u];
B = [x_B; u];

C = [x_A; d];
D = [x_B; d];

So now , I rotate the two lines, I use my function rotateTwoPoints(), by simply calling the following two lines :

[Af,Bf] = rotateTwoPoints(A,B,-45,O,true);
[Cf,Df] = rotateTwoPoints(C,D,-45,O,true);

Where O is the origin around which the rotation will be.

• I have tried O = [0;0] I mean it is the origin of the plot. No success !
• So I choose O the centroid of the image I using regionprops(I,"Centroid"). It was wrong because the centroid is not a center.
• Now, I use the center of the image O = floor(size(I)/2+0.5)' or using ceil !

But when I draw the resulting lines AfBf & CfDf in image I like this :

plot([Af(1) Bf(1)],[Af(2) Bf(2)],'k');
plot([Cf(1) Df(1)],[Cf(2) Df(2)],'k'); 

I get a result that is not correct !

Problem : In I_R, the AB & CD contain what I call the BlueZone (see image 3). But the rotated back lines AfBf & CfDf do not cover it !

The results in Images

Here is the rotated image I_R and the two lines drawn (The two middle red lines correspond to AB and CD) :

I then draw the rotated lines AfBf & CfDf in original image I (the black bold point correspond to the center on which I have done the rotation) :

IMAGE UPDATED

Problem : As you can see the BlueZone was inside the two lines AB and CD. But when rotated back it become outside, as shown in the following image (red arrows point to BlueZone) :

---

UPDATE ADDED a SNIPPET

Since my problem is not yet resolved , I selected the code that causes the problem and I add it as the following snippet (there is a variable stored in a file that you can download here) :

function Question()

% load image in I, the image is available online in the below link  
load I ;

% rotate I with -45° using imrotate
I_R = imrotate(I,-45);

% some data
x_A = 3 ;
x_B = 79;

u = 24;
d = 44;

% some meaningful Points : A,B,C and D that form two lines AB and CD
% parallels
A = [x_A; u];
B = [x_B; u];

C = [x_A; d];
D = [x_B; d];

% figure 1 contain two subplots 
figure(1);
% draw rotated image I_R
subplot(1,2,1), axis image, imagesc(I_R), hold on;
% draw two lines AB and CD in red in rotated image 
plot([A(1) B(1)],[A(2) B(2)],'r');
plot([C(1) D(1)],[C(2) D(2)],'r');
title('I_R the rotated image with the two lines AB and CD');

% draw original image I
subplot(1,2,2), axis image, imagesc(I)  , hold on;

% compute the middle of image I
axises=axis;
center = [mean(axises(1:2)),mean(axises(3:4))]';
% draw the center in red and as a point
plot(center(1),center(2),'ro');

% rotate the two lines, the result is the two lines AfBf and CfDf
[Af,Bf] = rotateTwoPoints(A,B,-45,center,true);
[Cf,Df] = rotateTwoPoints(C,D,-45,center,true);

% draw the rotated back lines in original image I
figure(1);
subplot(1,2,2);
plot([Af(1) Bf(1)],[Af(2) Bf(2)],'k');
plot([Cf(1) Df(1)],[Cf(2) Df(2)],'k');
title('the original image I with the two lines AfBf and CfDf');

function [Af,Bf] = rotateTwoPoints (A,B,t,Origin,isPlot)

% Definition of the rotation matrix (rotation around origin)
R=[ ...
    cosd(t) -sind(t)
    sind(t) cosd(t)
    ];

% translation 
At = A - Origin;
Bt = B - Origin;

% rotation of the points A and B
Ar = R*At;
Br = R*Bt;

% translation 
Af = Ar + Origin;
Bf = Br + Origin;

if isPlot == true

    figure(100)

    % Plot of the original line
    plot(A(1),A(2),'k*', B(1),B(2),'b*');
    line([A(1) B(1)],[A(2) B(2)], 'Color','r');

    grid on
    hold on

    % Plot the Origin around which the rotation will be
    plot(Origin(1),Origin(2),'k*','LineWidth',3);

    % Plot of the rotated line
    plot(Af(1),Af(2),'g*', Bf(1),Bf(2),'r*');
    line([Af(1) Bf(1)],[Af(2) Bf(2)], 'Color','b');
    legend('A','B','line AB','Origin','Af','Bf','line AfBf',['angle: ',num2str(t)],'Location','northeastoutside');
    daspect([1 1 1])

end

^{PS: I am using MATLAB R2012b}

Answer 1

As you say the rotation is being computed correctly, as shown by the first plot. The problem then is just with the final display of the results. When you do this

plot([Af(1) Bf(1)],[Af(2) Bf(2)],'k');
plot([Cf(1) Df(1)],[Cf(2) Bf(2)],'k'); 

you have a typo in the second line (second element of second argument) - you're plotting Bf(2) as the end of the second line instead of Df(2). When I replaced this with Df(2) it plotted parallel lines as expected.

UPDATE:

In the comments I suggested refactoring this code to use transform matrices for everything, such that the same set of transforms can be applied to both the image and the overlays. Here is a generic outline of how you might set that up.

There are two main points to be aware of.

1. Matlab's convention with images is that y=0 is the top. imshow and related then put y=0 at the top and y increasing down the plot. This is the opposite convention from plot, and so when you overlay plots and images one of them has to be inverted. This has consequences for the use of transform matrices, the main one being that the rotation angle has to be inverted.

2. Images in matlab can be thought of as the z-data for a regular grid of pixels. However those pixel coordinates aren't stored, they're just inferred when you display the image. Consequently, when you translate an image and then imshow it again the translation is not evident, because it re-infers a new set of pixel coordinates which happen to be the same as the ones for the untransformed image. To get round this we must fix spatial referencing to an image.

In addition to this the 2d transforms are carried out using 3space matrices so we need to define out points in 3space (but then truncate them for display in 2d plots).

So let's put it all together.

% get an image
I = imread('cameraman.tif');
% get spatial referencing information
Rcb = imref2d(size(I));

% define some test points in 3space
pta = [10; 10; 0];
ptb = [ 50;  10; 0];
% construct a test line in 2space from our test points
testline = [pta(1:2) ptb(1:2)];

% overlay line on plot
figure(90); clf
imshow(I,Rcb); truesize; hold on;
plot(testline(1, :), testline(2, :), 'y', 'linewidth', 4)

% define our transforms
% rotation angle (deg)
t = 45;

% a transform suitable for images
Rimage=[ ...
    cosd(t) -sind(t) 0
    sind(t) cosd(t) 0
    0 0 1
    ];
% the same transform suitable for points
Rpoints=[ ...
    cosd(-t) -sind(-t) 0
    sind(-t) cosd(-t) 0
    0 0 1
    ];

% make tform object suitable for imwarp
tform = affine2d(Rimage);

% transform image and spatial referencing with tform
[Ir, Rr] = imwarp(I, tform);

% transform points directly using matrix multplication
ptar = Rpoints*pta;
ptbr = Rpoints*ptb;

% construct the rotated line for plotting
newline = [ptar(1:2) ptbr(1:2)];

% the results
figure(91); clf
imshow(Ir, Rr); truesize; 
hold on;
plot(testline(1, :), testline(2, :), 'y', 'linewidth', 4)
plot(newline(1, :), newline(2, :), 'g', 'linewidth', 4)

In the second plot I've plotted the transormed line in green and the original line in yellow, for comparison.

Points to note:

1. The imwarp function doesn't use the rotation matrix directly, you have to construct a tform object from it first and then supply that. You can of course do all this in one line.

2. Having two matrices is a bit inelegant. It's okay here where we only have a single rotation to worry about, but the whole point of using transform matrices is that you can chain a sequence of transforms together just by multiplying the matrices. If you have to do this twice it would spoil an otherwise elegant piece of code and that would never do, so probably the cleanest way of doing it would be to flip your images right at the start of the whole process, then flip them back at the very end if necessary (for exporting, say).

3. The bookkeeping for the point data is somewhat tedious. There are lots of ways to do this, depending on the conventions you choose to adopt with respect to which coordinate is being held in which position in the column vectors. I've never found a set that plays nicely with both 3space transforms and plot, and what works best changes depending on the application. Helper functions can save some headaches for you, if not the next maintainer, and wrapping it all up in a class is the neatest, if you can justify the time.

UPDATE2:

To minimise the use of 3space coordinates I would simply use 2space vectors to describe all your points, as normal, and then add the dummy z-coordinate only when you need to perform the transform.

So you might have

testpoint = [1; 5]; % x and y xoordinates only
trpoint = Rpoints*[testpoint; 1];
trpoint = trpoint(1:2);

You could put this in a wrapper but there's no getting round the fact that you need to use a 3x3 matrix for the imwarp section, which means you need to specify your coordinates in 3space too.

Alternatively you could truncate the rotation matrix for the coordinate transform:

trpoint = Rpoints(1:2, 1:2)*testpoint;

but one way or another you'll have to do some index bookkeeping.

UPDATE3:

The OP doesn't have access to imref2d, so here's a hacky way to achieve the same result. The important line is this one

imshow(Ir, Rr);

where Rr is the spstial referencing object output by the imwarp function. If you don't have this you can supply spatial referencing to imshow manually with the 'XData' and 'YData' arguments. You need to work out the extents first, and decide on a convention. imref2d uses the convention of (0,0) as top left and when the image is rotated this corner of the original image remains as coordinate (0,0), such that the rotated image now extends from y=[-181:181] and x=[0:363]. To get these values you need to transform all four corners of the image and then work out the maxima and minima.

xmax = size(I,2);
ymax = size(I,2);
corner1 = [xmax; 0; 0]+0.5;
corner2 = [xmax; ymax; 0]+0.5;
corner3 = [0; ymax; 0]+0.5;
corner4 = [0; 0; 0]+0.5;
tc1 = Rpoints*corner1;
tc2 = Rpoints*corner2;
tc3 = Rpoints*corner3;
tc4 = Rpoints*corner4;
corners = [tc1 tc2 tc3 tc4];
xlims = minmax(corners(1,:));
ylims = minmax(corners(2,:));

you can then replace the imshow line with this one

imshow(Ir, 'xdata', xlims, 'ydata', ylims);

everything else should be the same.

Note that, when I did all this, there was a slight difference between the values I got above and the ones produced by imwarp in the imref2d object. That function returned

XWorldLimits: [0.2264 363.2264]
YWorldLimits: [-181.5000 181.5000]

Whereas mine gives

xlims: [0.7071  362.7458];
ylims: [-181.0193  181.0193];

I can't account for this without looking in the source of imwarp and, even if it isn't MEX at that level, I'm not sure how far you could go before you could be accused of IP infringement.

These numbers account for the placement of the image on the axes, so if you're using the image to choose points they may be off by a few tenths of a pixel. If you're just using the image for reference then it should be close enough.

Answer 2

1) you have image I

2)you rotated image I about its center (-45 Degree)=> I_R is rotated image

3)you draw two parallel lines in I_R (line1 & line2)

4)you rotate lines in respect to center of imge I

step number 4 is wrong you must choose center of image I_R for rotation center.

From Matlab Docs: By defaul imrotate creates an output image large enough to include the entire original image. pixels that fall outside the boundary of original image are set to 0 and apear as black background in the output image.

>

from comments: 1) you rotated the image about center point of "I" (-45degree)

2) when you do this using imrotate , the function rotate and translate the main picture and generates "I_R" ,

3) so the center point of I_R is the center point of "I" + some translation!.

4) so the main problem is calculating amount of translation,

5) to do this you can use this: Calculating translation value and rotation angle of a rotated 2D image or this: http://angeljohnsy.blogspot.com/2011/06/image-rotation.html or

6) simply you can use some markings on your original image (for example marking on the center of original image ) to see how much it translated ( or to see center of I_R directly) and the use this center for rotation back of lines.

Hope you find some good way to calculate translation , even you can write your own "image_rotation" script it is not so hard .

Answer 3

I found the solution !

Well, what happens is that the built-in function imrotate does not only rotate by default the image, but it does also a translation so that the rotated-image fit in the figure.

the solution is to use :

I_R = imrotate(I,-45,'nearest','crop');

the important parameter here is 'crop' (bbox) :

B = imrotate(A,angle,method,bbox) rotates image A, where bbox specifies the size of the returned image. bbox is a text string that can have one of the following values. The default value is enclosed in braces ({}).

'crop' Make output image B the same size as the input image A, cropping the rotated image to fit

{'loose'} Make output image B large enough to contain the entire rotated image. B is generally larger than A.

from the help of the function imrotate here.

Answer 4

Might not fix it, but try defining your matrix as:

R=[ ...
    cosd(t) -sind(t)
    -sind(t) -cosd(t)
    ];

When plotting over an image you need to take care of the y-axis, which automatically goes from top to bottom.

Answer 5

I may be misunderstanding something, but what you want to do is to turn a picture, draw two lines and then turn everything including the two lines back again, right?

If so I guess that the easiest way to turn the two lines would be to translate them into polar coordinates using car2pol. http://se.mathworks.com/help/matlab/ref/cart2pol.html

Then it should be easy to turn the lines. Let me know if I misunderstood anything.

edit: I tried to make some example code, apparently I am too tired to think straight and there is some error, but I guess the idea is clear. Also please do not use structures as I did, I just thought it would be the easiest way to name the variables.

%% Load picture
close all
[pic] = dicomread('C0011866_00094.DCM');

%% Get four points
    figure, imshow(pic,[]); 
    title('Select four points'); hold on; 
    [x,y] = ginput;
    
%% Plot the lines on top of the figure
hold on
plot(x(1:2),y(1:2),'b');
hold on
plot(x(3:4),y(3:4),'r');

%% Define the lines
%In cartesian coordinates
line1.start.cart = [x(1),y(1)];
line1.end.cart = [x(2),y(2)];

line2.start.cart = [x(3),y(3)];
line2.end.cart = [x(4),y(4)];

%Middle of picture
axises=axis;
center = [mean(axises(1:2)),mean(axises(3:4))];
plot(center(1),center(2),'ro')

%In polar coordinates from center [angle,radius]
[line1.start.pol(1),line1.start.pol(2)] = cart2pol(line1.start.cart(1)-center(1),line1.start.cart(2)-center(2));
[line1.end.pol(1),line1.end.pol(2)] = cart2pol(line1.end.cart(1)-center(1),line1.end.cart(2)-center(2));

[line2.start.pol(1),line2.start.pol(2)] = cart2pol(line2.start.cart(1)-center(1),line2.start.cart(2)-center(2));
[line2.end.pol(1),line2.end.pol(2)] = cart2pol(line2.end.cart(1)-center(1),line2.end.cart(2)-center(2));

%% Rotate the image a degrees
a = 10;
newpic = imrotate(pic,a);

%% Rotate the lines
line1.start.pol(1) = line1.start.pol(1) + a;
line1.end.pol(1) = line1.end.pol(1) + a;

line2.start.pol(1) = line2.start.pol(1) + a;
line2.end.pol(1) = line2.end.pol(1) + a;

%% Transform the coordinates back to cartesian
[line1.start.cart(1), line1.start.cart(2)] = pol2cart(line1.start.pol(1),line1.start.pol(2));
line1.start.cart = line1.start.cart+center;
line1.end.cart = line1.end.cart+center;
[line1.end.cart(1), line1.end.cart(2)] = pol2cart(line1.end.pol(1),line1.end.pol(2));
line2.start.cart = line2.start.cart+center;
line2.end.cart = line2.end.cart+center;

[line2.start.cart(1), line2.start.cart(2)] = pol2cart(line2.start.pol(1),line2.start.pol(2));
[line2.end.cart(1), line2.end.cart(2)] = pol2cart(line2.end.pol(1),line2.end.pol(2));

%% Plot it again
figure
imshow(newpic,[]);
hold on
plot([line1.start.cart(1),line1.end.cart(1)],[line1.start.cart(2),line1.end.cart(2)])
hold on
plot([line2.start.cart(1),line2.end.cart(1)],[line2.start.cart(2),line2.end.cart(2)])
hold on
plot(center(2),center(1),'ro')