小白学习图像处理3——图像旋转原理

一、图像旋转的原理

  图像的旋转类似坐标平面中XOY点的旋转，如下图，点P以坐标原点O为旋转中心，逆时针旋转角度 β 后得到点Q：

  设P、Q点的坐标分别为 (x, y) 和 (x’, y’)，假设P点离坐标原点的距离为D，则通过三角函数有：

$$x^{\prime }=Dcos(\alpha +\beta )=D(cos\alpha cos\beta -sin\alpha sin\beta )=xcos\beta -ysin\beta$$
$$y^{\prime }=Dsin(\alpha +\beta )=D(cos\alpha sin\beta +sin\alpha cos\beta )=xsin\beta +ycos\beta$$

  在数字图像中，由于图像是用二维数组存储的，因此上面的公式可以转化为矩阵运算：

$$R=\left[\begin{array}{ccc}cos\beta & sin\beta & 0\\ -sin\beta & cos\beta & 0\\ 0& 0& 1\end{array}\right]$$

$$\left[\begin{array}{ccc}{x}^{\prime }& y^{\prime }& 1\end{array}\right]=\left[\begin{array}{ccc}x& y& 1\end{array}\right]\ast R$$

---

二、使用matlab实现

1、思路

步骤：

1. 构建新的矩阵存储旋转得到的图像
2. 选择旋转中心
3. 计算原图像每一个像素点经旋转得到的坐标
4. 将原图像的像素点值映射到新图像中

问题1：如何确定像素点旋转后的坐标：假设原坐标为P: (x, y)，旋转中心为C: (cx, cy)，新坐标为Q: (x’, y’)，则

C = [cx, cy];
R = [cos(theta) sin(theta); -sin(theta) cos(theta)];
Q = round(R * (P - C) + C);

问题2：如何确定旋转后得到的图像的大小（可以存放下整个旋转后图像的矩形框按大小）：
思路：计算原图像四个角上的像素点旋转后的位置，计算行列的覆盖范围：

nw = round(([0, 0]-C)*R + C);
ne = round(([0, width]-C)*R + C);
sw = round(([height, 0]-C)*R + C);
se = round(([height, width]-C)*R + C);
h = max([nw(1), ne(1), sw(1), se(1)]) - min([nw(1), ne(1), sw(1), se(1)]);
w = max([nw(2), ne(2), sw(2), se(2)]) - min([nw(2), ne(2), sw(2), se(2)]);

---

2、实现代码

根据上述分析，代码如下：

img = imread('lena_color_512.tif');
[height, width,  color] =  size(img);

theta = 1/3*pi;
C = round([height/2, width/2]);     % 原图像中心
R = [cos(theta) sin(theta); -sin(theta) cos(theta)];

% 计算新图像大小
nw = round(([0, 0]-C)*R + C);
ne = round(([0, width]-C)*R + C);
sw = round(([height, 0]-C)*R + C);
se = round(([height, width]-C)*R + C);
h = max([nw(1), ne(1), sw(1), se(1)]) - min([nw(1), ne(1), sw(1), se(1)]);
w = max([nw(2), ne(2), sw(2), se(2)]) - min([nw(2), ne(2), sw(2), se(2)]);

target = zeros(h, w, 3, 'uint8') + 255;
c = round([h/2, w/2]);      % 新图像中心

% 旋转
for x = 1 : height
    for y = 1 : width
        P = [x, y];
        Q = round((P-C)*R + c);
        if (Q(1)>0 && Q(1)<=h) && (Q(2)>0 && Q(2)<=w)   % 越界判断
            target(Q(1), Q(2), :) = img(x, y, :);
        end
    end
end
imshow(target)

得到的结果和我们预料的有些不太一样，有很多小白点，主要原因是根据原图像坐标旋转计算得到的新坐标经过四舍五入后可能一样，导致一些坐标点未赋值，出现白点

---

三、优化

1、思路

  由于小白点是像素点未赋值造成的，我们可以遍历新图像的所有像素点，通过逆向寻找原图像对应像素点坐标，从而实现对所有像素点赋值。
问题1：如何根据新图像的坐标点找到原图中对应的坐标点
思路：考虑逆矩阵，从原图到旋转后的图像使用的是变换R，则根据新图像的坐标找到原图像对应坐标可以用R的逆矩阵

R = inv(R);	% 求逆
% 对新图像的坐标 Q(x, y)，求原图像坐标P
for x = 1 : h
    for y = 1 : w
        Q = [x, y];
        % 根据目标像素点计算原图像像素点
        P = round((Q-c)*R + C);
        if (P(1)>0 && P(1)<=height) && (P(2)>0 && P(2)<=width)   % 越界判断
            target(x, y, :) = img(P(1), P(2), :);
        end
    end
end

---

2、代码实现

img = imread('lena_color_512.tif');
[height, width,  color] =  size(img);

theta = 1/4*pi;
C = round([height/2, width/2]);     % 原图像中心
R = [cos(theta) sin(theta); -sin(theta) cos(theta)];
R = inv(R);     % 逆矩阵

% 原图像四个角旋转后的坐标
nw = round(([0, 0]-C)*R + C);
ne = round(([0, width]-C)*R + C);
sw = round(([height, 0]-C)*R + C);
se = round(([height, width]-C)*R + C);

% 新图像
h = max([nw(1), ne(1), sw(1), se(1)]) - min([nw(1), ne(1), sw(1), se(1)]);
w = max([nw(2), ne(2), sw(2), se(2)]) - min([nw(2), ne(2), sw(2), se(2)]);
target = zeros(h, w, 3, 'uint8') + 255;
c = round([h/2, w/2]);      % 新图像中心

for x = 1 : h
    for y = 1 : w
        Q = [x, y];
        % 根据目标像素点计算原图像像素点
        P = round((Q-c)*R + C);
        if (P(1)>0 && P(1)<=height) && (P(2)>0 && P(2)<=width)   % 越界判断
            target(x, y, :) = img(P(1), P(2), :);
        end
    end
end
imshow(target)

结果很棒！

---

3、使用双线性插值

  使用双线性插值是图像过度更加平滑，原来程序的for循环中，是通过 round 函数四舍五入得到原图像坐标，现在使用双线性插值（双线性插值的原理可以参考：图像插值）
  我们修改for循环的内容：

for x = 1 : h
    for y = 1 : w
        Q = [x, y];
        % 使用双线性插值
        P = (Q-c)*R + C;
        p1 = [ceil(P(1)), ceil(P(2))];
        p2 = [ceil(P(1)), ceil(P(2))+1];
        p3 = [ceil(P(1))+1, ceil(P(2))];
        p4 = [ceil(P(1))+1, ceil(P(2))];
        dx = ceil(P(1)) - P(1);
        dy = ceil(P(2)) - P(2);
        if (P(1)>=1 && P(1)<height-1) && (P(2)>=0 && P(2)<width-1)   % 越界判断
            target(x, y, :) = (1-dx)*(1-dy)*img(p1(1), p1(2), :) + (1-dx)*dy*img(p3(1), p3(2), :)...
                + dx*(1-dy)*img(p2(1), p2(2), :) + dx*dy*img(p4(1), p4(2), :);
        end
    end
end

---

四、matlab函数实现图像旋转

  其实，matlab有自带的图像旋转函数 imrotate（会不会觉得前面全部白学了，hhh），举个例子：

1、imrotate函数

%% imrotate
img=imread('lena_color_512.tif');
figure
imshow(img)
figure
% 对I旋转逆时针30°，使用最邻近插值（也可以用bilinear）
imshow(imrotate(img, 30, 'nearest'))

---

2、imtransform函数

%% imtransform
I = imread('lena_color_512.tif');
theta = -pi/6;
T = [cos(theta) sin(theta) 0; -sin(theta) cos(theta) 0; 0 0 1];	% 构造旋转变换矩阵
tform = maketform('affine',T);	% 生成tform结构体
J = imtransform(I,tform);	% 调用tform进行旋转
imshow(J)

这两个函数执行结果都是一样的：

---

3、图像块的旋转

拓展 ：我们使用blockproc函数实现对图像块的操作，blockproc使用方法如下：

% A 代表图像, [m, n]代表块的大小, fun是函数句柄，表示对每个块的处理函数
B = blockproc(A,[m n],fun);
% processes the image A by applying the function fun to each distinct block of size [m n] and concatenating the results into the output matrix, B.

拓展 ：函数句柄相当于C++里面的指针，可以自定义函数句柄：

add = @(a, b) a+b;
class(add)  % function_handle
add(2, 9)	% 相当于 a+b

实现图像块旋转的代码很简洁：

img = imread('lena_color_512.tif');
img = im2double(img);

fun = @(block_struct)...	% 换行符
    imrotate(block_struct.data, 30);	
img1 = blockproc(img, [64, 64], fun);
figure, imshow(img1)

效果如下：

---

  以上就是全部的内容啦，介绍了旋转图像的原理和matlab的imrotate函数的使用，通过类比，我们可以知道一些其他图像空间变换的原理，比如平移、放大、缩小、水平/垂直偏移等…

完结 cheers 🍻~