图像旋转的原理
图像旋转的原理其实很简单，为了简化公式的推导，这里我们假设绕原点 ( 0 , 0 ) (0,0) (0,0)旋转。
本文规定逆时针旋转为正，当然你也可以规定顺时针为正（此时，正角度就表示顺时针旋转）
很容易得到三个方程：
t
a
n
(
θ
+
α
)
=
y
′
x
′
(
1
)
tan(\theta + \alpha)={y' \over x'} \quad \quad \quad (1)
tan(θ+α)=x′y′(1)
t
a
n
α
=
y
x
(
2
)
tan\alpha={y \over x} \quad \quad \quad (2)
tanα=xy(2)
x
2
+
y
2
=
x
′
2
+
y
′
2
(
3
)
x^2+y^2={x'}^2+{y'}^2 \quad \quad \quad (3)
x2+y2=x′2+y′2(3)
由和角公式可知
t
a
n
(
θ
+
α
)
=
t
a
n
θ
+
t
a
n
α
1
−
t
a
n
θ
t
a
n
α
tan(\theta + \alpha)={{tan\theta+tan\alpha} \over {1-tan\theta tan\alpha}}
tan(θ+α)=1−tanθtanαtanθ+tanα
代入（1）和（2）可得
y
′
=
x
′
(
y
+
x
t
a
n
θ
)
x
−
y
t
a
n
θ
(
4
)
y'={{x'(y+xtan\theta)} \over {x-ytan\theta}} \quad \quad \quad (4)
y′=x−ytanθx′(y+xtanθ)(4)
(4)代入(3)得到：
x
′
2
=
(
x
c
o
s
θ
−
y
s
i
n
θ
)
2
(
5
)
x'^2=(xcos\theta-ysin\theta)^2 \quad \quad \quad (5)
x′2=(xcosθ−ysinθ)2(5)
（5）代入（3）
y
′
2
=
(
x
s
i
n
θ
+
y
c
o
s
θ
)
2
y'^2=(xsin\theta+ycos\theta)^2
y′2=(xsinθ+ycosθ)2
所以
{
x
′
=
x
c
o
s
θ
+
y
(
−
s
i
n
θ
)
y
′
=
x
s
i
n
θ
+
y
c
o
s
θ
(
6
)
\begin{cases} x '=xcos\theta +y(-sin\theta) \\ y'=xsin\theta + ycos \theta \end{cases} \quad \quad \quad (6)
{x′=xcosθ+y(−sinθ)y′=xsinθ+ycosθ(6)
我们知道仿射变换矩阵可以表示为：
A
=
[
c
o
s
θ
−
s
i
n
θ
s
i
n
θ
c
o
s
θ
]
(7)
A=\left[ \begin{matrix} cos\theta & -sin\theta \\ sin\theta & cos \theta \\ \end{matrix} \right] \tag{7}
A=[cosθsinθ−sinθcosθ](7)
公式7与公式6在形式上是一致的(这里没有加入平移和缩放)
注意公式6中
x
,
y
,
x
′
,
y
′
x,y,x',y'
x,y,x′,y′实际表示
x
−
0
,
y
−
0
,
x
′
−
0
,
y
′
−
0
x-0,y-0,x'-0,y'-0
x−0,y−0,x′−0,y′−0即与旋转中心横坐标以及纵坐标之间的距离，现在我们将旋转中心变为更具有一般性的
(
x
0
,
y
0
)
(x_0,y_0)
(x0,y0),得旋转公式为：
{
x
′
=
(
x
−
x
0
)
c
o
s
θ
+
(
y
−
y
0
)
(
−
s
i
n
θ
)
+
x
0
y
′
=
(
x
−
x
0
)
s
i
n
θ
+
(
y
−
y
0
)
c
o
s
θ
+
y
0
(
8
)
\begin{cases} x '=(x-x_0)cos\theta +(y-y_0)(-sin\theta)+x_0 \\ y'=(x-x_0)sin\theta + (y-y_0)cos \theta+y_0 \end{cases} \quad \quad \quad (8)
{x′=(x−x0)cosθ+(y−y0)(−sinθ)+x0y′=(x−x0)sinθ+(y−y0)cosθ+y0(8)
对于图像来说，规定顺时针旋转角度为正，则旋转公式与上述公式一致，由于任意角度的旋转后会出现像素坐标为负的情况，如果不将旋转后的图像坐标平移，会缺失部分图像，所以，对于旋转后的图像，我们通常会加入一个平移
{
x
′
=
(
x
−
x
0
)
c
o
s
θ
+
(
y
−
y
0
)
(
−
s
i
n
θ
)
+
x
0
+
Δ
x
y
′
=
(
x
−
x
0
)
s
i
n
θ
+
(
y
−
y
0
)
c
o
s
θ
+
y
0
+
Δ
y
(
9
)
\begin{cases} x '=(x-x_0)cos\theta +(y-y_0)(-sin\theta)+x_0+\Delta x \\ y'=(x-x_0)sin\theta + (y-y_0)cos \theta+y_0+\Delta y \end{cases} \quad \quad \quad (9)
{x′=(x−x0)cosθ+(y−y0)(−sinθ)+x0+Δxy′=(x−x0)sinθ+(y−y0)cosθ+y0+Δy(9)
图像中使用后向映射，使用
x
′
,
y
′
x',y'
x′,y′表示
x
,
y
x,y
x,y,则上述公式变为
{
x
=
(
(
x
′
−
x
0
−
Δ
x
)
c
o
s
θ
+
(
y
′
−
y
0
−
Δ
y
)
s
i
n
θ
)
/
s
c
a
l
e
+
x
0
y
=
(
(
x
′
−
x
0
−
Δ
x
)
(
−
s
i
n
θ
)
+
(
y
′
−
y
0
−
Δ
y
)
c
o
s
θ
)
/
s
c
a
l
e
+
y
0
(
10
)
\begin{cases} x=((x'-x_0-\Delta x)cos\theta +(y'-y_0-\Delta y)sin\theta)/scale+x_0 \\ y=((x'-x_0-\Delta x)(-sin\theta) + (y'-y_0-\Delta y)cos \theta)/scale+y_0 \end{cases} \quad \quad \quad (10)
{x=((x′−x0−Δx)cosθ+(y′−y0−Δy)sinθ)/scale+x0y=((x′−x0−Δx)(−sinθ)+(y′−y0−Δy)cosθ)/scale+y0(10)
这里加入了缩放，其中：scale>1表示原图放大 <1表示原图缩小
有了上面的公式，下面就可以写出相应的代码了
注：为什么公式推导的时候选用右手系？因为右手系与图像坐标系推导出来的公式在形式上是统一的，而右手系又是我们熟悉的坐标系。
图像旋转的实现
最近邻插值
先给出一个结构最简洁的版本，采用最近邻插值
// 宏定义
#define DEGREE2RADIAN(x) (x*CV_PI/180)//角度转弧度
#define RADIAN2DEGREE(x) (x*180/CV_PI)//弧度转角度
#define SHIFT 10
#define DESCALE(x,n) (((x)+(1 << ((n)-1))) >> (n))
/* center:原图像的旋转中心
dstSize:旋转后图像的大小
theta:旋转角度，单位弧度，顺时针为正
scale:缩放，scale>1表示放大 <1表示缩小
*/
void Rotate_Nearest(const Mat &srcImage, Mat &dstImage, Point center, Size dstSize, double theta, double scale)
{
CV_Assert(srcImage.depth() == CV_8U);
dstImage.create(dstSize, srcImage.type());
int x0 = center.x;
int y0 = center.y;
theta = DEGREE2RADIAN(theta);
// dx,dy就是dst与src图像中心的距离
int dx = dstImage.cols/2 - srcImage.cols/2;
int dy = dstImage.rows/2 - srcImage.rows/2;
int numberOfChannels = srcImage.channels();
int widthOfDst = dstImage.cols;
int heightOfDst = dstImage.rows;
for (int y = 0; y <= heightOfDst - 1; ++y)
{
for (int x = 0; x <= widthOfDst - 1; ++x)
{
float srcX = ((x - x0 - dx)*cos(theta) + (y - y0 - dy)*sin(theta))/scale + x0;
float srcY = ((x0 + dx - x)*sin(theta) + (y - y0 - dy)*cos(theta))/scale + y0;
// get the nearest coordinate of src
int x1 = (int)srcX;
int y1 = (int)srcY;
if (numberOfChannels == 1)
{
if ((x1 >= 0 && x1 <= srcImage.cols - 1) && (y1 >= 0 && y1 <= srcImage.rows - 1))
{
dstImage.at<uchar>(y, x) = srcImage.at<uchar>(y1, x1);
}
else
{
// 越界赋值0
dstImage.at<uchar>(y, x) = 0;
}
}
else
{
if ((x1 >= 0 && x1 <= srcImage.cols - 1) && (y1 >= 0 && y1 <= srcImage.rows - 1))
{
dstImage.at<cv::Vec3b>(y, x) = srcImage.at<cv::Vec3b>(y1, x1);
}
else
{
dstImage.at<cv::Vec3b>(y, x) = cv::Vec3b(0,0,0);
}
}
}
}
}
测试使用的原始图像大小为1000 x 580 (三通道彩色图)
测试代码：
Rotate_Nearest(srcImage, dstImage, Point(srcImage.cols / 2, srcImage.rows / 2), Size(2500, 2500), 30.0, 2);
效果：
原图：
结果：
无缩放测试：
Rotate_Nearest(srcImage, dstImage, Point(srcImage.cols / 2, srcImage.rows / 2), Size(2500, 2500), 30.0, 1);
结果：
双线性插值
下面我们对他进行优化，首先我们将最近邻插值改用更为通用的双线性插值，代码如下
void Rotate_Bilinear(const Mat &srcImage, Mat &dstImage, Point center, Size dstSize, double theta, double scale)
{
CV_Assert(srcImage.depth() == CV_8U);
dstImage.create(dstSize, srcImage.type());
dstImage.setTo(Scalar(0, 0, 0));
int x0 = center.x;
int y0 = center.y;
theta = DEGREE2RADIAN(theta);
//
Mat extendedImage;
copyMakeBorder(srcImage, extendedImage, 1, 1, 1, 1, BORDER_CONSTANT,Scalar(0,0,0)); // 使用0填充边界
// dx,dy就是dst与src图像中心的距离
int dx = dstImage.cols / 2 - srcImage.cols / 2;
int dy = dstImage.rows / 2 - srcImage.rows / 2;
int numberOfChannels = srcImage.channels();
int widthOfDst = dstImage.cols;
int heightOfDst = dstImage.rows;
for (int y = 0; y <= heightOfDst - 1; ++y)
{
for (int x = 0; x <= widthOfDst - 1; ++x)
{
// 按照原来的方式计算原图坐标
float srcX = ((x - x0 - dx)*cos(theta) + (y - y0 - dy)*sin(theta)) / scale + x0;
float srcY = ((x0 + dx - x)*sin(theta) + (y - y0 - dy)*cos(theta)) / scale + y0;
// 加1,得到在extendedImage中的坐标
srcX++;
srcY++;
// get the nearest coordinate of src
int x1 = (int)(srcX);
int y1 = (int)(srcY);
// 浮点转化为整数
int dx1 = (srcX - x1)*(1<< SHIFT);
int dy1 = (srcY - y1)*(1<< SHIFT);
if (numberOfChannels == 1)
{
// !！！注意这里的范围，在extendedImage中，原图的范围就是1~cols - 2了
if ((x1 >= 1 && x1 <= extendedImage.cols - 2) && (y1 >= 1 && y1 <= extendedImage.rows - 2))
{
//双线性插值
//周围4个点
//a就是最近邻像素
//a b
// p
//c d
uchar a = extendedImage.at<uchar>(y1, x1);
uchar b = extendedImage.at<uchar>(y1, x1 + 1);
uchar c = extendedImage.at<uchar>(y1 + 1, x1);
uchar d = extendedImage.at<uchar>(y1 + 1, x1 + 1);
//int p = (a*((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b*dx1*((1 << SHIFT) - dy1) + c*((1 << SHIFT) - dx1)*dy1 + d*dx1*dy1)/(1<<(2* SHIFT));
int p = a*((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b*dx1*((1 << SHIFT) - dy1) + c*((1 << SHIFT) - dx1)*dy1 + d*dx1*dy1;
p = DESCALE(p, 2*SHIFT);
dstImage.at<uchar>(y, x) = p;
}
else
{
// 越界赋值0
dstImage.at<uchar>(y, x) = 0;
}
}
else
{
if ((x1 >= 1 && x1 <= extendedImage.cols - 2) && (y1 >= 1 && y1 <= extendedImage.rows - 2))
{
//双线性插值
//周围4个点
//a就是最近邻像素
//a b
// p
//c d
Vec3b a = extendedImage.at<Vec3b>(y1, x1);
Vec3b b = extendedImage.at<Vec3b>(y1, x1 + 1);
Vec3b c = extendedImage.at<Vec3b>(y1 + 1, x1);
Vec3b d = extendedImage.at<Vec3b>(y1 + 1, x1 + 1);
/*int p1 = (a[0] * ((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[0] * dx1*((1 << SHIFT) - dy1) + c[0] * ((1 << SHIFT) - dx1)*dy1 + d[0] * dx1*dy1)/(1<<(2*SHIFT));
int p2 = (a[1] * ((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[1] * dx1*((1 << SHIFT) - dy1) + c[1] * ((1 << SHIFT) - dx1)*dy1 + d[1] * dx1*dy1)/ (1 << (2 * SHIFT));
int p3 = (a[2] * ((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[2] * dx1*((1 << SHIFT) - dy1) + c[2] * ((1 << SHIFT) - dx1)*dy1 + d[2] * dx1*dy1)/ (1 << (2 * SHIFT));*/
int p1 = a[0]*((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[0]*dx1*((1 << SHIFT) - dy1) + c[0]*((1 << SHIFT) - dx1)*dy1 + d[0]*dx1*dy1;
p1 = DESCALE(p1, 2 * SHIFT);
int p2 = a[1]*((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[1]*dx1*((1 << SHIFT) - dy1) + c[1]*((1 << SHIFT) - dx1)*dy1 + d[1] *dx1*dy1;
p2 = DESCALE(p2, 2 * SHIFT);
int p3 = a[2]*((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[2]*dx1*((1 << SHIFT) - dy1) + c[2]*((1 << SHIFT) - dx1)*dy1 + d[2] *dx1*dy1;
p3 = DESCALE(p3, 2 * SHIFT);
dstImage.at<cv::Vec3b>(y, x) = Vec3b(p1,p2,p3);
}
else
{
dstImage.at<cv::Vec3b>(y, x) = cv::Vec3b(0, 0, 0);
}
}
}
}
}
测试方法与上述相同，测试代码如下：
Rotate_Bilinear(srcImage, dstImage, Point(srcImage.cols / 2, srcImage.rows / 2), Size(2500, 2500),30.0, 2);
基本旋转效果都是一样的，下面我们看局部放大图
这是最近邻的眼部区域
双线性在保持图像细节方面要好于最近邻，而且不容易产生锯齿
双线性的优化
上面的双线性插值速度比较慢，我们对其进行优化，主要优化有：
1.将循环内部不变量提取出来
由于sin和cos的计算是比较慢的函数，所以可以将他们提前计算好，而不是每次循环都要计算
如：
double sinTheta = sin(theta);
double cosTheta = cos(theta);
2.改变循环体内循环变量自增的方式
采用加法自增的方式，而不是每次都要计算乘法，详见代码
注：浮点数转化为整数的优化，上面已经涉及，这里不再说明
优化后的代码如下：
void Rotate_Bilinear2(const Mat &srcImage, Mat &dstImage, Point center, Size dstSize, double theta, double scale)
{
CV_Assert(srcImage.depth() == CV_8U);
dstImage.create(dstSize, srcImage.type());
dstImage.setTo(Scalar(0, 0, 0));
int x0 = center.x;
int y0 = center.y;
theta = DEGREE2RADIAN(theta);
//
Mat extendedImage;
copyMakeBorder(srcImage, extendedImage, 1, 1, 1, 1, BORDER_CONSTANT, Scalar(0, 0, 0)); // 使用0填充边界
// dx,dy就是dst与src图像中心的距离
int dx = dstImage.cols / 2 - srcImage.cols / 2;
int dy = dstImage.rows / 2 - srcImage.rows / 2;
int numberOfChannels = srcImage.channels();
int widthOfDst = dstImage.cols;
int heightOfDst = dstImage.rows;
// 优化部分/
// 将循环内的不变量提取出来
double sinTheta = sin(theta);
double cosTheta = cos(theta);
scale = 1.0 / scale;
// 改变了循环内部增量的方式
double temp1= (0 - y0 - dy)*sinTheta;
double temp2 = (0 - y0 - dy)*cosTheta;
double dtemp1 = sinTheta;
double dtemp2 = cosTheta;
for (int y = 0; y <= heightOfDst - 1; ++y,temp1+=dtemp1,temp2+=dtemp2)
{
// 改变了循环内部增量的方式
double temp3= ((0 - x0 - dx)*cosTheta + temp1)*scale + x0;
double temp4= (-(0 - x0 - dx)*sinTheta + temp2)*scale + y0;
double dtemp3 = (cosTheta)*scale;
double dtemp4= (-sinTheta)*scale;
for (int x = 0; x <= widthOfDst - 1; ++x,temp3+=dtemp3,temp4+=dtemp4)
{
// 计算原图坐标
double srcX = temp3;
double srcY = temp4;
// 加1,得到在extendedImage中的坐标
srcX++;
srcY++;
// get the nearest coordinate of src
int x1 = (int)(srcX);
int y1 = (int)(srcY);
// 浮点转化为整数
int dx1 = (srcX - x1)*(1 << SHIFT);
int dy1 = (srcY - y1)*(1 << SHIFT);
if (numberOfChannels == 1)
{
// !！！注意这里的范围，在extendedImage中，原图的范围就是1~cols - 2了
if ((x1 >= 1 && x1 <= extendedImage.cols - 2) && (y1 >= 1 && y1 <= extendedImage.rows - 2))
{
//双线性插值
//周围4个点
//a就是最近邻像素
//a b
// p
//c d
uchar a = extendedImage.at<uchar>(y1, x1);
uchar b = extendedImage.at<uchar>(y1, x1 + 1);
uchar c = extendedImage.at<uchar>(y1 + 1, x1);
uchar d = extendedImage.at<uchar>(y1 + 1, x1 + 1);
//int p = (a*((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b*dx1*((1 << SHIFT) - dy1) + c*((1 << SHIFT) - dx1)*dy1 + d*dx1*dy1)/(1<<(2* SHIFT));
int p = a*((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b*dx1*((1 << SHIFT) - dy1) + c*((1 << SHIFT) - dx1)*dy1 + d*dx1*dy1;
p = DESCALE(p, 2 * SHIFT);
dstImage.at<uchar>(y, x) = p;
}
else
{
// 越界赋值0
dstImage.at<uchar>(y, x) = 0;
}
}
else
{
if ((x1 >= 1 && x1 <= extendedImage.cols - 2) && (y1 >= 1 && y1 <= extendedImage.rows - 2))
{
//双线性插值
//周围4个点
//a就是最近邻像素
//a b
// p
//c d
Vec3b a = extendedImage.at<Vec3b>(y1, x1);
Vec3b b = extendedImage.at<Vec3b>(y1, x1 + 1);
Vec3b c = extendedImage.at<Vec3b>(y1 + 1, x1);
Vec3b d = extendedImage.at<Vec3b>(y1 + 1, x1 + 1);
/*int p1 = (a[0] * ((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[0] * dx1*((1 << SHIFT) - dy1) + c[0] * ((1 << SHIFT) - dx1)*dy1 + d[0] * dx1*dy1)/(1<<(2*SHIFT));
int p2 = (a[1] * ((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[1] * dx1*((1 << SHIFT) - dy1) + c[1] * ((1 << SHIFT) - dx1)*dy1 + d[1] * dx1*dy1)/ (1 << (2 * SHIFT));
int p3 = (a[2] * ((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[2] * dx1*((1 << SHIFT) - dy1) + c[2] * ((1 << SHIFT) - dx1)*dy1 + d[2] * dx1*dy1)/ (1 << (2 * SHIFT));*/
int p1 = a[0] * ((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[0] * dx1*((1 << SHIFT) - dy1) + c[0] * ((1 << SHIFT) - dx1)*dy1 + d[0] * dx1*dy1;
p1 = DESCALE(p1, 2 * SHIFT);
int p2 = a[1] * ((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[1] * dx1*((1 << SHIFT) - dy1) + c[1] * ((1 << SHIFT) - dx1)*dy1 + d[1] * dx1*dy1;
p2 = DESCALE(p2, 2 * SHIFT);
int p3 = a[2] * ((1 << SHIFT) - dx1)*((1 << SHIFT) - dy1) + b[2] * dx1*((1 << SHIFT) - dy1) + c[2] * ((1 << SHIFT) - dx1)*dy1 + d[2] * dx1*dy1;
p3 = DESCALE(p3, 2 * SHIFT);
dstImage.at<cv::Vec3b>(y, x) = Vec3b(p1, p2, p3);
}
else
{
dstImage.at<cv::Vec3b>(y, x) = cv::Vec3b(0, 0, 0);
}
}
}
}
}
下面我们测试一下他们的速度
测试环境：Intel core i5-6200U,12G，放大2倍
| 测试方法 | Rotate_Nearest | Rotate_Bilinear | Rotate_Bilinear2 |
|---|---|---|---|
| 速度(单位ms) | 58.2 | 144.9 | 78.7 |
可以看出，优化后的代码速度提升还是很明显的。
双线性代码还可以进一步的优化，有兴趣的朋友可以自己实现。
完整工程见github项目：QQImageProcess_OpenCV
其中图像旋转优化的实现在 Src/ImageProcess/GeometryTransformation.h中的Rotate_系列函数中。
2016-9-11 15:12:24
Last Updated:2017-3-11 23:22:16
非常感谢您的阅读，如果您觉得这篇文章对您有帮助，欢迎扫码进行赞赏。
5639
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
