图像处理作业2-对图片进行任意角度旋转

Programming Assignment 2

Author：Tian YJ
(1) 2D Rotation: Write a MATLAB function(rotate.m) that takes a set of points in 2D and an angle (in degrees) and returns a new set of points which have been rotated counter-clockwise by that angle.

(2) Image Rotation: Write a MATLAB function (rotate_image.m) that takes an image and an angle (in degrees) and returns a new image which has been rotated counter-clockwise around the center of the image by that angle.

(3) User Interaction: Write a MATLAB script (stored as an ascii text file, with the name straighten.m) that loads in an image, allows the user to click on two points in the image, and then rotates the image so that the line connecting these two points is horizontal. This script will use your rotate_image.m function as a subroutine to generate the rotated image.

(4) 提交报告内容：包括实现原理，程序输入与输出对比，结果分析，及代码。上面以matlab为例进行说明，但不限编程语言。

实现原理：

如图所示

如图中所示，旋转角度为$heta$角，根据三角函数：
$X_0=Rcosalpha$
$Y_0=Rsinalpha$
旋转后：$X=Rcos(alpha- heta)$，$Y=Rsin(alpha- heta)$
根据正弦加法定理和余弦加法定理可知：
$cos(alpha- heta)=cosalpha cos heta+sinalpha sin heta$
$sin(alpha- heta)=sinalpha cos heta-cosalpha sin heta$
则$X=Rcosalpha cos heta+Rsinalpha sin heta=X_0cos heta+Y_0sin heta$
$Y=Rsinalpha cos heta-Rcosalpha sin heta=Y_0 cos heta-X_0 sin heta$
$left[ egin{array}{cccc} X& Y& 1 end{array} ight ] =left[ egin{array}{cccc} X_0& Y_0& 1 end{array} ight] left[ egin{array}{cccc} cos heta& -sin heta& 0\ sin heta& cos heta& 0\ 0& 0& 1\ end{array} ight ]$
有了上面的公式，还需要进行逆运算。因为，在进行位图旋转时，在确定旋转后的位图大小后，可以确定旋转后的位图矩形区域，需要访问该区域的每一个坐标，以获得该坐标对应于源位图的坐标，最后获得该坐标对应于源位图坐标点处的像素数据。因此需要进行逆运算。
$X=X_0cos heta+Y_0sin heta$
有$Y_0=frac{X-Xcos heta}{siin heta}$
因为$Y=Y_0cos heta-X_0sin heta$
有$frac{X-Xcos heta}{sin heta}=frac{Y+X_0sin heta}{cos heta}$
将其展开有$Xcos heta-Ysin^2 heta=Ysin heta+X_0sin^2 heta$求得$Xcos heta-Ysin heta=X_0$即$X_0=Xcos heta-Ysin heta$
同理$Y_0=Xsin heta+Ycos heta$
这样对应的逆运算矩阵为
$left[ egin{array}{cccc} X_0& Y_0& 1 end{array} ight ] =left[ egin{array}{cccc} X& Y& 1 end{array} ight] left[ egin{array}{cccc} cos heta& sin heta& 0\ -sin heta& cos heta& 0\ 0& 0& 1\ end{array} ight ]$
有了转换公式，我们还需要进行坐标转换。以一幅图为例，在计算机中的原点为图像的左上角。而在数学中，坐标原点为图像的中心点。假设位图的高度为H，宽度为W，则图像坐标$(X_0，Y_0)$与数学坐标$(X，Y)$的关系是
$X=X_0-0.5W$，$Y=-Y_0+0.5H$对应的矩阵表示为
$left[ egin{array}{cccc} X& Y& 1 end{array} ight ] =left[ egin{array}{cccc} X_0& Y_0& 1 end{array} ight] left[ egin{array}{cccc} 1& 0& 0\ 0& -1& 0\ -0.5W& 0.5H& 1\ end{array} ight ]$
逆运算为$X_0=X+0.5W$，$Y_0=-Y+0.5H$
$left[ egin{array}{cccc} X_0& Y_0& 1 end{array} ight ] =left[ egin{array}{cccc} X& Y& 1 end{array} ight] left[ egin{array}{cccc} 1& 0& 0\ 0& -1& 0\ 0.5W& 0.5H& 1\ end{array} ight ]$
在经过数学公式旋转后，还需要数学坐标转换为新位图的图像坐标。转换的公式为：$X=X_0+0.5hat{W}$，$Y=-Y_0+0.5hat{H}$
其中，$hat{W}$和$hat{H}$表示旋转位图的高度和宽度。矩阵表示为
$left[ egin{array}{cccc} X& Y& 1 end{array} ight ] =left[ egin{array}{cccc} X_0& Y_0& 1 end{array} ight] left[ egin{array}{cccc} 1& 0& 0\ 0& -1& 0\ 0.5hat{W}& 0.5hat{H}& 1\ end{array} ight ]$
逆运算公式为：$X_0=X+0.5hat{W}$，$Y_0=-Y+0.5hat{H}$
逆运算矩阵为
$left[ egin{array}{cccc} X_0& Y_0& 1 end{array} ight ] =left[ egin{array}{cccc} X& Y& 1 end{array} ight] left[ egin{array}{cccc} 1& 0& 0\ 0& -1& 0\ -0.5hat{W}& 0.5hat{H}& 1\ end{array} ight ]$
有了上述的旋转公式，下面介绍图像旋转的过程。
（1）将原始图像的坐标系转换为数学坐标系
（2）通过旋转公式对图像坐标进行旋转
（3）将旋转后的数学坐标系转换为图像坐标系
下面给出位图旋转的矩阵描述
$left[ egin{array}{cccc} X& Y& 1 end{array} ight ] =left[ egin{array}{cccc} X_0& Y_0& 1 end{array} ight] left[ egin{array}{cccc} 1& 0& 0\ 0& -1& 0\ -0.5W& 0.5H& 1\ end{array} ight ] left[ egin{array}{cccc} cos heta& -sin heta& 0\ sin heta& cos heta& 0\ 0& 0& 1\ end{array} ight ] left[ egin{array}{cccc} 1& 0& 0\ 0& -1& 0\ 0.5hat{W}& 0.5hat{H}& 1\ end{array} ight ]$
根据矩阵乘法可得：

$left[ egin{array}{cccc} cos heta& sin heta& 0\ -sin heta& cos heta& 0\ -0.5Wcos heta+0.5Hsin heta+0.5hat{W}& -0.5Wsin heta-0.5Hcos heta+0.5hat{H}& 1\ end{array} ight ]$
图像逆运算的过程如下：
（1）在目标区域中，将图像坐标系转换为数学坐标系。
（2）利用图像旋转的逆运算公式（之前导出的）确定当前坐标点在原图像中的坐标。
（3）将获得的坐标转换为源图像中的图像坐标系。
这样就获得了目标区域坐标点对应的源图像中的坐标点。
$left[ egin{array}{cccc} X_0& Y_0& 1 end{array} ight ] =left[ egin{array}{cccc} X& Y& 1 end{array} ight] left[ egin{array}{cccc} 1& 0& 0\ 0& -1& 0\ -0.5hat{W}& 0.5hat{H}& 1\ end{array} ight ] left[ egin{array}{cccc} cos heta& sin heta& 0\ -sin heta& cos heta& 0\ 0& 0& 1\ end{array} ight ] left[ egin{array}{cccc} 1& 0& 0\ 0& -1& 0\ 0.5W& 0.5H& 1\ end{array} ight ]$

根据矩阵乘法可得：

$left[ egin{array}{cccc} cos heta& -sin heta& 0\ sin heta& cos heta& 0\ -0.5hat{W}cos heta+0.5hat{H}sin heta+0.W& 0.5hat{W}sin heta-0.5hat{H}cos heta+0.5H& 1\ end{array} ight ]$
对应的函数关系式为：
$X_0=Xcos heta+Ysin heta-0.5hat{W}cos heta-0.5hat{H}sin heta+0.5W$
$Y_0=-Xsin heta+Ycos heta+0.5hat{W}sin heta-0.5hat{H}cos heta+0.5H$
设$dx=-0.5hat{W}cos heta-0.5hat{H}sin heta+0.5w$
设$dy=0.5hat{W}sin heta-0.5hat{H}cos heta+0.5H$
则有$X_0=Xcos heta+Ysin heta+dx$
$Y_0=-Xsin heta+Ysin heta+dy$

（1）尝试以前行映射的思路实现，也就是把原本进行向量的旋转，找到旋转后的向量位置，然后将原图的像素值赋值过去。
代码如下：

function [ img_new ] = rotate_forward_mapping(img,angle)         %定义旋转函数,这里尝试前向映射
angle = angle / 180 * pi;                                        %进行角度转弧度

[H,W,C] = size(img);                                             %获取输入图像的行、列和通道数

H_new = round(H*abs(cos(angle))+ W*abs(sin(angle)));             %旋转图像后得到的新高度
W_new = round(W*abs(cos(angle))+ H*abs(sin(angle)));             %旋转图像后得到的新宽度
img_new = zeros(H_new,W_new,C);                                  %定义生成目标图像的行列以及通道数

M1 = [1 0 0; 0 -1 0; -0.5*W 0.5*H 1];                            %坐标系变换矩阵M1
M2 = [cos(angle) -sin(angle) 0; sin(angle) cos(angle) 0; 0 0 1]; %角度旋转变换矩阵M2
M3 = [1 0 0; 0 -1 0; 0.5*W_new 0.5*H_new 1];                     %坐标系变换矩阵M3

for i = 1:W
    for j = 1:H
        temp = [i j 1]*M1*M2*M3;                                 %得到旋转后的矩阵temp
        x = round(temp(1));                                      %x取矩阵temp的第一行第一列,x对应i，为宽度                   
        y = round(temp(2));                                      %y取矩阵temp的第一行第二列,y对应j，为高度
        if x==0 || y==0
            x = x+1;
            y = y+1;
        end
        img_new(y, x, :) = img(j, i, :);                        
    end
end

输入图片试试看效果：

clc                                 
img = imread('crooked_horizon.jpg');
figure,imshow(img); 
title('原始图像');
img_new = rotate_forward_mapping(img,-30);    %调用rotate()函数逆时针旋转30° 
figure,imshow(uint8(img_new));
title('逆时针旋转30°');

输出结果如下：

原始图像	逆时针旋转30度

可以看到有很多瑕疵，图像中会出现很多噪声很多杂点，出现杂点的原因是从原图旋转后的像素位置在原图可能找不到，解决方法是用逆向思维，接下来我又尝试使用反向映射来实现。

（2）反向映射即从旋转后的图像出发，找到对应的原图像的点，然后将原图像中的灰度值传递过来，这样旋转后的图像的每个像素肯定可以对应到原图像中的一个点，采取不同策略可以让像素对应地更加准确，并基于上次实现的双线性插值算法在图像旋转中的应用来更好地实现图片旋转，减少噪声和杂点。

代码如下：

function [ img_new ] = rotate_reverse_mapping(img,angle)         %定义旋转函数,这里尝试反向映射
angle = angle / 180 * pi;                                        %进行角度转弧度

[H,W,C] = size(img);                                             %获取输入图像的行、列和通道数

H_new = round(H*abs(cos(angle))+ W*abs(sin(angle)));             %旋转图像后得到的新高度
W_new = round(W*abs(cos(angle))+ H*abs(sin(angle)));             %旋转图像后得到的新宽度
img_new = zeros(H_new,W_new,C);                                  %定义生成目标图像的行列以及通道数

M1 = [1 0 0; 0 -1 0; -0.5*W_new 0.5*H_new 1];                    %坐标系变换矩阵M1
M2 = [cos(angle) sin(angle) 0; -sin(angle) cos(angle) 0; 0 0 1]; %角度旋转变换矩阵M2
M3 = [1 0 0; 0 -1 0; 0.5*W 0.5*H 1];                             %坐标系变换矩阵M3

for i = 1:W_new
   for j = 1:H_new
      temp = [i j 1]*M1*M2*M3;                                   %得到旋转后的矩阵temp
      x = round(temp(1));                                        %x取矩阵temp的第一行第一列,x对应i，为宽度
      y = round(temp(2));                                        %y取矩阵temp的第一行第二列,y对应j，为高度
      
      if y < 1 || x < 1 || y > H || x > W                        %防止边界溢出
          img_new(j, i) = 0;
      else
                                                                 %双线性插值
                                                                 %计算四个角点的坐标(四邻域)
          left = floor(x);                                       %左边界
          right = ceil(x);                                       %右边界
          top = floor(y);                                        %上边界
          bottom = ceil(y);                                      %下边界
                                                
          a = x - left;                                          %取小数部分
          b = y - top;                                           %取小数部分
                                                                 %双线性插值计算
          img_new(j, i, 1) = (1-a)*(1-b)*img(top, left, 1) + a*(1-b)*img(top, right, 1) + ...
              (1-a)*b*img(bottom, left, 1) + a*b*img(bottom, right, 1);
          img_new(j, i, 2) = (1-a)*(1-b)*img(top, left, 2) + a*(1-b)*img(top, right, 2) + ...
              (1-a)*b*img(bottom, left, 2) + a*b*img(bottom, right, 2);
          img_new(j, i, 3) = (1-a)*(1-b)*img(top, left, 3) + a*(1-b)*img(top, right, 3) + ...
              (1-a)*b*img(bottom, left, 2) + a*b*img(bottom, right, 3);
      end
   end
end

输入图片试试看效果：

clc                                 
img = imread('crooked_horizon.jpg');
figure,imshow(img); 
title('原始图像');
img_new1 = rotate_reverse_mapping(img,-30);    %调用rotate()函数逆时针旋转30° 
figure,imshow(uint8(img_new1));
title('逆时针旋转30°');

输出结果如下

原始图像	逆时针旋转30度

再对比一下两种方法的结果：

前向映射	反向映射

可以清楚地看到，通过反向映射很好地解决了前向映射出现有规律噪声的问题，图片看起来更加清晰。总的来说，效果还行。

（3）接下来实现的是交互式的输入旋转角度，进行旋转，因为反向映射的效果更佳，所以接下来都会基于反向映射的方法来实现。

clc
clear
img = imread('crooked_horizon.jpg');
[H,W,C] = size(img);              %获取图像尺寸
figure,imshow(img);               % title('原始图像');
hold on;
[a, b]=ginput();                  %(a(1),b(1))为下点，(a(2),b(2))为上点
plot(a(1),b(1), 'or');
plot(a(2),b(2), 'or');            %在图上标注鼠标点击的点
                                  %进行位图坐标转换
x1 = a(1)-0.5*W;
y1 = -b(1)+0.5*H;
x2 = a(2)-0.5*W;
y2 = -b(2)+0.5*H;

r = atan((y2-y1)/(x2-x1));                    %由两点斜率求倾斜角
d = r*180/pi;                                 %弧度转角度
img_new = rotate_reverse_mapping(img, -d);    %调用函数逆时针旋转 
figure,imshow(uint8(img_new));

结果展示

鼠标点击图（1）	旋转结果（1）

鼠标点击图（2）	旋转结果（2）

我的实验分析已在上述描述中，总的来说，这次实验还算顺趟，基本达到预期结果！