用双线性插值做旋转缩放操作

缩放操作

缩放无非就是放大和缩小两种,不管是哪种操作,都可以得到长宽各自缩放的比值。得到的新图要么是像素增加了(拉长),要么是像素减少了(缩短)。不管是哪种情况,都需要对新图中的像素值赋值(也叫插值)。本文只介绍一下最常见的最邻近插值和双线性插值。

什么是插值

首先讲讲什么叫插值。简单地说,给缩放后得到的新图的每一个像素点赋值,就叫做插值。最常见的插值方法分两种:

一、向前映射:遍历原图中的每一个像素点,根据映射关系可以计算出它们在新图中的位置(遇到小数的情况要取整),这个新位置可以直接取原图对应点的像素值,也可以取一个接近的值。这种方法可能导致新图中存在空洞点,即原图中的所有像素点都无法映射到这个位置点,另一个缺陷是,原图中两个不同的像素点可能映射到同一个位置,这样肯定会有一个像素的信息丢失了;

二、向后映射:这里我们是反过来求出新图到原图的映射变换(按照线性代数的知识,就是将向前映射取一个逆),然后遍历新图中的每个像素,计算出它们在原图中的对应位置,再取值。同样的,这个位置可能不是整数,简单的做法可以直接取整,比较好的做法是综合一下周围点的像素值取一个均值。前一种做法称为最邻近插值,后一种做法根据综合像素的方式又分为几种方法,这里只讲最常见的双线性插值。

最邻近插值

我们就拿一个方向的缩放为例。

屏幕快照 2016-07-18 下午10.25.07
屏幕快照 2016-07-18 下午10.25.07

假设右边两个像素点的图是原图,左边四个像素点的图是新图,对于这种宽度扩大两倍,高度保持不变的情况,最邻近插值的做法是这样:我们先计算出新图回到原图的缩放因子是 0.5,所以,对新图中的每个像素,我们计算出它们乘以 0.5 后得到的位置,这个位置是它们对应于原图的位置。新图中,1、2 像素点的位置分别为 (0,0)、(0,1),缩小后得到的新位置为 (0,0),也就是它们对应原图中的 1 像素点,所以它们的像素值都取原图中 1 像素点的值,同理,3、4 像素点的值取原图中 2 像素点的值。这样,新图的像素肯定能被插满。对于原图缩小的情况道理是一样的。

这种做法的优点是简单粗暴,计算量小且容易理解,缺点是由于太过简单粗暴,导致没有很好的利用原图的信息,比方说上图,原图被拉长后,1 和 2 之间的像素应该要有一个过渡,这样人眼看的时候才不会觉得突兀(或者说变化太快),图像中常见的锯齿就是这种由突兀产生的。

上一段简单的代码,用到了 CImg 库(太简单不解释):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
CImg<unsigned char> nearest_scale(const CImg<unsigned char>& srcImg, 
int width, int height) {
int srcWidth = srcImg.width();
int srcHeight = srcImg.height();
CImg<unsigned char> outImg(width, height, 1, 3, 0);
// 缩放因子
double width_scale = srcWidth / (double)width;
double height_scale = srcHeight / (double)height;
for (int r = 0; r < height; r++) {
for (int c = 0; c < width; c++) {
outImg(c, r, 0, 0) = srcImg((int)(c*width_scale),
(int)(r*height_scale), 0, 0);
outImg(c, r, 0, 1) = srcImg((int)(c*width_scale),
(int)(r*height_scale), 0, 1);
outImg(c, r, 0, 2) = srcImg((int)(c*width_scale),
(int)(r*height_scale), 0, 2);
}
}
return outImg;
}

双线性插值

既然有双线性插值,我们不妨先看下单线性插值。

上图中,我们一样假设右图为原图,左图为缩放后的新图。我们通过计算缩放因子,发现左图中 2 这个像素点还原到新图后,对应的坐标不是整数,它恰好落在原图两个像素点之间 0.75 : 0.25 的位置(注意最邻近插值有个取整操作,把位置上的这点细节忽略了)。此时,为了更好地取值,我们不是简单粗暴地用原图中的一个像素值去取代新图,而是根据距离的比值取了一个综合的结果。线性插值的做法是,根据距离比 0.75 : 0.25,对原像素值取一个加权平均,距离近的,像素值要接近,所以权值更大,应该取 0.75,距离远的则相反,取 0.25,这样一来,原图中的两个像素值,根据距离的反比取一个加权平均,就是新图中 2 这个像素的值了。双线性插值道理同上,只不过推广到二维的情况:

除了水平方向,我们发现高度其实也不是整数,双线性插值会进一步再考虑高度的距离比,假设宽度的距离比是 u : (1-u),高度的距离比是 v : (1-v),此时,我们要综合考虑周围四个点的情况,比方说,对于原图右下角的像素值,对应的加权值为 (1-0.25)*(1-0.25),而右上角的像素值,对应的加权值为 (1-0.75)*(1-0.25),左边的类似,最后将这些值全部加起来就是新的像素值了。

示例代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
CImg<unsigned char> bilinear_scale(const CImg<unsigned char>& srcImg,
int width, int height) {
int srcWidth = srcImg.width();
int srcHeight = srcImg.height();
CImg<unsigned char> outImg(width, height, 1, 3, 0);
// 缩放因子
double width_scale = srcWidth / (double)width;
double height_scale = srcHeight / (double)height;
// 新图对应原图的坐标
double srcX, srcY, u, v;
for (int r = 0; r < height; r++) {
for (int c = 0; c < width; c++) {
srcX = c * width_scale;
srcY = r * height_scale;
u = srcX - (int)srcX;
v = srcY - (int)srcY;
for (int channel = 0; channel < 3; channel++) {
outImg(c, r, 0, channel) =
(int)((1-u)*(1-v)*srcImg(valueWidth(srcX, srcWidth),
valueHeight(srcY, srcHeight), 0, channel)
+(1-u)*v*srcImg(valueWidth(srcX, srcWidth),
valueHeight(srcY+1, srcHeight), 0, channel)
+u*(1-v)*srcImg(valueWidth(srcX+1, srcWidth),
valueHeight(srcY, srcHeight), 0, channel)
+u*v*srcImg(valueWidth(srcX+1, srcWidth),
valueHeight(srcY+1, srcHeight), 0, channel));
}
}
}
return outImg;
}

旋转操作

旋转跟缩放本质上是一样的,只是映射变换不同。

先看看映射函数怎么求。

假设点 (\(x_0\), \(y_0\)) 逆时针旋转到 (\(x_1\), \(y_1\)) 的位置,旋转角度为 A,(\(x_1\), \(y_1\)) 与 x 轴正方向的夹角为 B,且知道点 (\(x_0\), \(y_0\)) 到原点的距离为 r (这个量只是起到中间变量的作用)。现在需要求出 \(x_0\)\(y_0\)\(x_1\)\(y_1\) 之间的映射关系。

根据高中三角函数的知识,可以列出如下等式: \[x_0=r cos(B-A)=rcosBcosA+rsinBsinA \tag{1}\] \[y_0=r sin(B-A)=r sinBcosA-rsinAcosB \tag{2}\] \[x_1=r cosB \tag{3}\] \[y_1=r sinB \tag{4}\] 将 (3) (4) 式分别代入 (1) (2) 得到: \[x_0 = x_1 cosA + y_1 sinA \tag{5} \] \[y_0=y_1cosA-x_1sinA \tag{6} \] 这样,向后映射的表达式我们就求出来了。

接着,由 (5) (6) 式分别可得: \[y_1=\frac{x_0-x_1cosA}{sinA} \tag{7}\] \[y_1=\frac{y_0+x_1sinA}{cosA} \tag{8}\] 再由 (7) = (8) 可以解出: \[ x_1=x_0*cosA-y_0*sinA \] 同样的方法可以解出: \[ y_1=y_0*cosA+x_0*sinA \] 这样,向前映射也得到了。之后,按照缩放里面的思路,我们只要先算出新图的大小(根据向前映射可以求得),然后遍历新图里的像素,根据向后映射反推在原图中的像素位置,再用插值的方法就可以给每个像素点赋值了。要注意的一点是,我们这里是以原点为旋转中心计算出来的表达式,如果直接套用上面的计算结果,那么你得到的新图像是原图绕左上角的点旋转A度角后得到的图(旋转方向取决于 A 的符号)。这是因为计算机中,图片的坐标原点在左上角。如果要让图片让自己的中心旋转,必须先将图片的中心与左上角原点对齐(只需做一下平移即可),旋转完成之后,再按照之前平移的距离方向重新调整位置。

示例代码(采用双线性插值):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
/* 
* 旋转图片,旋转角为theta,并采用双线性插值减少锯齿
*/
CImg<unsigned char> rotate_biliinear(const CImg<unsigned char>& srcImg, double theta) {
int width = srcImg.width();
int height = srcImg.height();
// 原图的四个顶点坐标,这里以图片中心为坐标原点
Position lt(0-width/2, 0+height/2), lb(0-width/2, 0-height/2),
rt(0+width/2, 0+height/2), rb(0+width/2, 0-height/2);
// 获得旋转后的图片的四个顶点坐标
Position new_lt((int)(lt.x*cos(theta)+lt.y*sin(theta)), (int)(lt.y*cos(theta)-lt.x*sin(theta))),
new_lb((int)(lb.x*cos(theta)+lb.y*sin(theta)), (int)(lb.y*cos(theta)-lb.x*sin(theta))),
new_rt((int)(rt.x*cos(theta)+rt.y*sin(theta)), (int)(rt.y*cos(theta)-rt.x*sin(theta))),
new_rb((int)(rb.x*cos(theta)+rb.y*sin(theta)), (int)(rb.y*cos(theta)-rb.x*sin(theta)));
int newWidth = max(abs(new_rt.x-new_lb.x), abs(new_lt.x-new_rb.x));
int newHeight = max(abs(new_lt.y-new_rb.y), abs(new_lb.y-new_rt.y));

CImg<unsigned char> newImg(newWidth, newHeight, 1, 3, 0);
// 开始填充新图片的灰度值
bilinear_interpolation(newImg, srcImg, theta);

return newImg;
}


/**
* 采用双线性插值填充新图
*/
void bilinear_interpolation(CImg<unsigned char>& outImg,
const CImg<unsigned char>& srcImg, double theta) {
int newWidth = outImg.width(), newHeight = outImg.height();
int srcWidth = srcImg.width(), srcHeight = srcImg.height();
int halfW = newWidth / 2;
int halfH = newHeight / 2;
double srcX, srcY, u, v;
for (int r = 0; r < newHeight; r++) {
for (int c = 0; c < newWidth; c++) {
if (get_origin_pos(c-halfW, r-halfH, srcWidth, srcHeight, theta, srcX, srcY)) {
u = srcX - (int)srcX;
v = srcY - (int)srcY;
for (int channel = 0; channel < 3; channel++) {
outImg(c, r, 0, channel) =
(int)((1-u)*(1-v)*srcImg(valueWidth(srcX, srcWidth), valueHeight(srcY, srcHeight), 0, channel)
+(1-u)*v*srcImg(valueWidth(srcX, srcWidth), valueHeight(srcY+1, srcHeight), 0, channel)
+u*(1-v)*srcImg(valueWidth(srcX+1, srcWidth), valueHeight(srcY, srcHeight), 0, channel)
+u*v*srcImg(valueWidth(srcX+1, srcWidth), valueHeight(srcY+1, srcHeight), 0, channel));
}
}
}
}
}

// 获得旋转后图片的像素对应于原图的像素位置,用于双线性插值
bool get_origin_pos(int x, int y, int srcWidth, int srcHeight, double theta,
double& srcX, double& srcY) {
// 找到(x, y)在原图中对应的位置(srcX, srcY)
srcX = (double)x * cos(theta) - (double)y * sin(theta);
srcY = (double)x * sin(theta) + (double)y * cos(theta);
if (srcX >= (0-srcWidth/2-1) && srcX <= srcWidth/2+1 && srcY >= (0-srcHeight/2-1) && srcY <= srcHeight/2+1) {
srcX += srcWidth/2;
srcY += srcHeight/2;
return true;
} else {
return false;
}
}

// 检查像素位置,防止超过图片宽度
int valueWidth(double srcX, int width) {
if (srcX < 0) srcX = 0;
if (srcX >= width) srcX--;
return srcX;
}

// 检查像素位置,防止超过图片高度
int valueHeight(double srcY, int height) {
if (srcY < 0) srcY = 0;
if (srcY >= height) srcY--;
return srcY;
}

最后放上最邻近插值和双线性插值的旋转效果图,可以看到最邻近插值的锯齿略多于双线性插值

(原图)

(最邻近插值)

(双线性插值)

参考

图像旋转算法原理– 旋转矩阵

以上代码均在 https://github.com/Jermmy/ComputerVision_hw/tree/master/Ex1/src