OpenCV算法学习笔记之几何变换

发布于 2019-05-03  123 次阅读



此系列的其他文章:
OpenCV算法学习笔记之初识OpenCV
OpenCV算法学习笔记之对比度增强
OpenCV算法学习笔记之平滑算法

对于一张图片的放大、缩小、旋转等操作我们统称为几何变换。几何变换是图像最基本也是最成用的操作,常见的几何变换有仿射变换、投影变换、极坐标变换。

仿射变换

二维空间的仿射变换公式为:
\left(\begin{matrix}\bar{x}\\\bar{y}\end{matrix}\right) = \left( \begin{matrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{matrix}\right) \left(\begin{matrix}x\\y\end{matrix}\right)+\left(\begin{matrix}a_{13}\\a_{23}\end{matrix}\right)
此公式也可以表示为
\left(\begin{matrix}\bar{x}\\\bar{y}\\1\end{matrix}\right)=A\left(\begin{matrix}x\\y\\1\end{matrix}\right)
其中
A=\left(\begin{matrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\0&0&1\end{matrix}\right)
A通常称为仿射变换矩阵,由于最后一行为(0, 0, 1),所以在之后的讨论中会省略最后一行。常见的仿射变换类型有平移、缩放、旋转。

平移

在图像中,通常是取左上角为原点坐标,向右和向下为正方向。假设图像中的任一坐标为(x,y),假设图像向右平移t_x个单位,向下平移t_y个单位,则平移后的坐标(\bar{x},\bar{y})=(x+t_x,y+t_y),用矩阵表示就是
\left(\begin{matrix}\bar{x}\\\bar{y}\\\bar{1}\end{matrix}\right)=\left(\begin{matrix}1&0&t_x\\0&1&t_y\\0&0&1\end{matrix}\right)\left(\begin{matrix}x\\y\\1\end{matrix}\right)
t_x>0 ,则表示向右移动;若 t_y>0 ,则表示向下移动。

缩放

这里的缩放与我们平常的认知有所不同:(x,y)(0,0)为中心在水平方向上缩放s_x倍,在方向上垂直缩放上s_y倍,是指缩放后的坐标距离缩放中心(0,0)的水平垂直距离分别变为了原来的s_xs_y倍。若缩放中心为原点,则缩放公式为 (\bar{x},\bar{y})=(x* s_x,y* s_y) ,对应的矩阵表示则为:
\left(\begin{matrix}\bar{x}\\\bar{y}\\\bar{1}\end{matrix}\right)=\left(\begin{matrix}s_x&0&0\\0&s_y&0\\0&0&1\end{matrix}\right)\left(\begin{matrix}x\\y\\1\end{matrix}\right)
如果是以(x_0,y_0)为中心进行缩放变换,相当于先把原点平移到(x_0,y_0),然后以原点为中心进行变换,最后将原点再移回去。对应公式为 (\bar{x},\bar{y})=((x- x_0) * s_x + x_0,(y- y_0) *s_y + y_0) ,用矩阵表示为:
\left(\begin{matrix}\bar{x}\\\bar{y}\\\bar{1}\end{matrix}\right)=\left(\begin{matrix}1&0&x_0\\0&1&y_0\\0&0&1\end{matrix}\right)\left(\begin{matrix}s_x&0&0\\0&s_y&0\\0&0&1\end{matrix}\right)\left(\begin{matrix}1&0&-x_0\\0&1&-y_0\\0&0&1\end{matrix}\right)\left(\begin{matrix}x\\y\\1\end{matrix}\right)
以任意一点为中心的缩放变换矩阵是平移矩阵和以(0,0)为中心的缩放矩阵组合相乘得到的。

旋转

设坐标(x,y)(0,0)顺时针旋转到 (\bar{x},\bar{y}) ,角度从\alpha变为\alpha+\theta,cos\alpha=\frac{x}{p}sin\alpha=\frac{y}{p},其中p=\sqrt{x^2+y^2},则
cos(\alpha+\theta)=cos\alpha cos\theta-sin\alpha sin\theta=\frac{\bar{x}}{p}\\sin(\alpha+\theta)=sin\alpha cos\theta+sin\theta cos\alpha=\frac{\bar{y}}{p}
化简可得\bar{x}=xcos\theta-ysin\theta\bar{y}=xsin\theta+ycos\theta;相反,若左边(x,y)逆时针旋转到(\bar{x},\bar{y}),则取\theta-\theta即可。

矩阵表示为(顺时针):
\left(\begin{matrix}\bar{x}\\\bar{y}\\\bar{1}\end{matrix}\right)=\left(\begin{matrix}cos\theta&-sin\theta&0\\sin\theta&cos\theta&0\\0&0&1\end{matrix}\right)\left(\begin{matrix}x\\y\\1\end{matrix}\right)
若以任意一点(x_0,y_0)为中心旋转,相当于先将原点移动到旋转中心,然后绕原点旋转,最后移回坐标原点,用矩阵表示为:
\left(\begin{matrix}\bar{x}\\\bar{y}\\\bar{1}\end{matrix}\right)=\left(\begin{matrix}1&0&x_0\\0&1&y_0\\0&0&1\end{matrix}\right)\left(\begin{matrix}cos\theta&-sin\theta&0\\sin\theta&cos\theta&0\\0&0&1\end{matrix}\right)\left(\begin{matrix}1&0&-x_0\\0&1&-y_0\\0&0&1\end{matrix}\right)\left(\begin{matrix}x\\y\\1\end{matrix}\right)
上面的运算顺序是从右向左的。

OpenCV提供函数rotate(InputArray src, Output dst, int rotateCode)实现顺时针90°、180°、270°的旋转,rotateCode有以下取值:

ROTATE_90_CLOCKWISE    //顺时针旋转90度
ROTATE_180             //顺时针旋转180度
ROTATE_270_COUNTERCLOCKWISE  //顺时针旋转270度

需要注意的是OpenCV还有一个函数为flip(src, dst, int flipCode)实现了图像的水平镜像、垂直镜像和逆时针旋转180°,不过并不是通过仿射变换实现的,而是通过行列互换,它与rotate()transpose()函数一样都在core.hpp头文件中。

求解仿射变换矩阵

以上都是知道变换前坐标求变换后的坐标,如果我们已经知道了变换前的坐标和变换后的坐标,想求出仿射变换矩阵,可以通过解方程法或矩阵法。

由于仿射变换矩阵
A=\left(\begin{matrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\0&0&1\end{matrix}\right)
有6个未知数,所以我们只需三组坐标列出六个方程即可。

OpenCV提供函数getAffineTransform(src, dst)通过方程法求解,其中src和dst分别为前后坐标,函数声明在imgproc.hpp头文件,在Python中,可以用以下方式求解:

import cv2 as cv
import numpy as np
src = np.array([[0, 0], [200, 0], [0, 200]], np.float32)
dst = np.array([[0, 0], [100, 0], [0, 100]], np.float32)
A = cv.getAffineTransform(src, dst)

对于C++来说,一种方式是将坐标存在Point2f数组中,另一种方法是保存在Mat中:

// 第一种方法
Point2f src1[] = {Pointy2f(0, 0), Point2f(200, 0), Point2f(0, 200)};
Point2f dst1[] = {Pointy2f(0, 0), Point2f(100, 0), Point2f(0, 100)};
// 第二种方法
Mat src2 = (Mat_<float>(3, 2) << 0, 0, 200, 0, 0, 200);
Mat dst2 = (Mat_<float>(3, 2) << 0, 0, 100, 0, 0, 100);

Mat A = getAffineTransform(src1, dst1);

对于矩阵法求解,仿射变换矩阵是平移仿射矩阵乘以缩放仿射矩阵
\left(\begin{matrix}\bar{x}\\\bar{y}\\\bar{1}\end{matrix}\right)=\left(\begin{matrix}1&0&t_x\\0&1&t_y\\0&0&1\end{matrix}\right)\left(\begin{matrix}s_x&0&0\\0&s_y&0\\0&0&1\end{matrix}\right)\left(\begin{matrix}x\\y\\1\end{matrix}\right)
即运算的顺序是从右往左。对于矩阵的乘法,Numpy提供函数dot()实现,假设某个图像先等比例缩放2倍,然后水平向右移动100,垂直向下移动200,则

import numpy as np
# 缩放矩阵
s = np.array([[0.5, 0, 0],
              [0, 0.5, 0],
              [0, 0, 1.0]])
# 平移矩阵
t = np.array([[1, 0, 100],
              [0, 1, 200],
              [0, 0, 1]])
A = np.dot(s, t)

C++中OpenCV通过“*”或gemm()函数实现矩阵乘法,缩放矩阵和平移矩阵可以用Mat表示。

若是缩放后以(x_0,y_0)旋转,则通过以下公式求:
\left(\begin{matrix}1&0&x_0\\0&1&y_0\\0&0&1\end{matrix}\right)\left(\begin{matrix}cos\theta&-sin\theta&0\\sin\theta&cos\theta&0\\0&0&1\end{matrix}\right)\left(\begin{matrix}s_x&0&0\\0&s_y&0\\0&0&1\end{matrix}\right)\left(\begin{matrix}1&0&-x_0\\0&1&-y_0\\0&0&1\end{matrix}\right)
运算顺序仍是从右向左,如还需平移,则左乘平移仿射矩阵。

对于等比例缩放的仿射变换,OpenCV提供函数getRotationMatrix2D(center, angle, scale)来计算矩阵,center是变换中心;angle是逆时针旋转的角度,如果是负数就代表顺时针了;scale是等比例缩放的系数。

插值算法

在运算中,我们可能会遇到目标坐标有小数的情况,比如将坐标(3,3)缩放2倍变为了(1.5,1.5),但是对于图像来说并没有这个点,这时候我们就要用周围坐标的值来估算此位置的颜色,也就是插值。

最近邻插值

最近邻插值就是从(x,y)的四个相邻坐标中找到最近的那个来当作它的值,如(2.3,4.7),它的相邻坐标分别为(2,4)(3,4)(2,5)(3,5),计算(x_i,y_i)(i=1,2,3,4)(x,y)的距离,最近的为(2,5),则取(2,5)的值为(2.3,4.7)的值。

此种方法得到的图像会出现锯齿状外观,对于放大图像则更明显。

双线性插值

[x]表示不大于x的最大整数

第一步:|x-[x]|表示点(x,y)([x],[y])的水平距离,|[x]+1-x|表示点(x,y)([x]+1,[y])的水平距离,对于处于(x,[y])的值f_I(x,[y])用以下公式计算:
f_I (x,[y]) = |x-[x]| * f_I ([x+1],y) + (1-|x-[x]|) * f_I ([x],[y])
第二步:|x-[x]|表示点(x,y)([x],[y]+1)的水平距离,|[x]+1-x|表示点(x,y)([x]+1,[y]+1)的水平距离,对于处于(x,[y]+1)的值f_I(x,[y]+1)用以下公式计算:
f_I (x,[y]+1) = |x-[x]| * f_I ([x+1],[y]+1) + (1-|x-[x]|) * f_I ([x],[y]+1)
第三步:|y-[y]|表示点(x,y)(x,[y])的水平距离,|[y]+1-y|表示点(x,y)(x,[y]+1)的水平距离,结合一二步中求得的f_I(x,[y]+1)f_I(x,[y])计算f_I(x,y)的值
f_I (x,y) = |y-[y]| * f_I(x,[y]+1) + (1-|y-[y]|) * f_I (x,[y])

实现

在已知仿射矩阵的基础上,OpenCV提供函数warpAffine(src, M, dsize[, dst[, flags[, borderMode[, borderValue ]]]])实现图像的仿射变换,其中,src是输入图像矩阵;M是2×3的仿射矩阵;dsize是输出图像的大小(二元元组);flags是插值法,有INTE_NEARESTINTE_LINEAR(默认)等;borderMode是填充模式,有BORDER_CONSTANT等;borderValue是当borderMode=BORDER_CONSTANT的时候的值。

为了使用方便,OpenCV还提供了另一个函数resize(InputArray src, OutputArray dst,Size dsize, double fx=0, double fy=0, int interpolation=INTER_LINEAR)实现缩放,其中dsize是输出图像的大小,二元元组;fx、fy分别是水平垂直方向上的缩放比例,默认为0;interpolation是插值法。

Mat I = imread("test.png");
Mat dst;
resize(I, dst, Size(I.cols/2, I.rows/2), 0.5, 0.5);

//resize也可写成下面这种形式
//resize(I, dst, Size(), 0.5, 0.5);

投影变换

如果物体在三维空间中发生旋转,这种变换通常成为投影变换。由于可能出现阴影或者遮挡,所以投影变换很难修正。但是如果物体是平面的,那么很容易通过二维投影变换对此物体三维变换进行模型化,这就是专用的二维投影变换,可以通过以下公式表述:
\left(\begin{matrix}\bar{x}\\\bar{y}\\\bar{z}\end{matrix}\right)=\left(\begin{matrix}a_{11}&a_{12}&a_{13}\\a_{21}&a_{22}&a_{23}\\a_{31}&a_{32}&a_{33}\end{matrix}\right)\left(\begin{matrix}x\\y\\z\end{matrix}\right)
OpenCV提供函数getPerspectiveTransform(src, dst)实现求解投影矩阵,需要输入四组对应的坐标。该函数的Python的API,src和dst分别是4×2的二维ndarray,数据类型必须是float32,否则会报错;返回的矩阵是float64类型的。

OpenCV提供函数warpPerspective(src, M, dsize[, dst[, flags[, borderMode[, borderValue ]]]])实现投影变换,参数说明和仿射变化类似。

极坐标变换

通常通过极坐标变化校正图像中的圆形物体或包含在圆环中的物体。

笛卡尔坐标转化为极坐标

对于xoy平面上的任意一点(x,y),以(x_0,y_0)为中心的极坐标转换公式为
r= \sqrt{(x-x_0)^2 + (y-y_0)^2} \\ \theta = 2 \pi + arctan2(y-y_0, x-x_0 ), y-y_0<=0\\ arctan2(y-y_0, x - x_0 ),y-y_0>0
以变换中心为圆心的同一个圆上的点,在极坐标系\theta or显示为一条直线

可以用以下代码实现

import math
r = math.sqrt(math.pow(x-x0)+math.pow(y-y0))
theta = math.atan2(y-y0, x-x0)/math.pi*180 # 转化为角度

OpenCV提供函数cartToPolar(x, y[, magnitude[, angle[, angleInDegrees ]]])实现将原点移动到变换中心后的笛卡尔坐标向极坐标转换,返回值magnitude和angle是和参数x,y相同尺寸和类型的ndarray,angleInDegrees为True时返回的angle是角度,否则为弧度;x是数组且数据类型必须为浮点型、float32或float64,y尺寸和类型和x一致,x、y分别代表x坐标和y坐标。

极坐标转化为笛卡尔坐标

转换公式为
x=x_0+rcos\theta\\y=y_0+rsin\theta
OpenCV提供函数polarToCart(magnitude, angle[, x[, y[, angleInDegrees ]]]),返回的x,y是以原点为中心的笛卡尔坐标,即已知(\theta,r)(x_0,y_0),计算出的是(x-x_0,y-y_0);magnitude对应r,angle对应\theta;参数说明和cartToPolar类似。

举例已知\theta or坐标系中的(30, 10),(31, 10), (30, 11), (31, 11),\theta以角度表示,问笛卡尔坐标系中的哪四个坐标以(-12, 15)为中心经过极坐标变换后得到这四个坐标,实现代码为:

import cv2 as cv
import numpy as np
# 也可以用np.array([30, 31, 30, 31]),只影响输出下x,y的格式
angle = np.array([[30, 31], [30, 31]], np.float32) 
r = np.array([[10, 10], [11, 11]], np.float32)
x, y = cv.polarToCart(r, angle, angleInDegrees, True)
# 计算出的是(x-x0, y-y0)的坐标
x += -12
y += 15

如果用C++实现,可以将角度和距离放在Mat中。

极坐标变换处理图像

假设要将与(x_0,y_0)的距离范围为 [r_{min},r_{max}] ,角度范围在 [\theta_{min},\theta_{max}] 内的点进行极坐标向笛卡尔坐标的转换,输出图像\textbf{O}(i,j)的值用以下公式计算:
\textbf{O}(i,j)=f_{\textbf{I}}(x_0+(r_{min}+r_{step}i) * cos(\theta_{min}+\theta_{step}j),y_0+(r_{min}+r_{step}i) * sin(\theta_{min}+\theta_{step}j))
其中,0 < r_{step} \leq{1}代表步长,\theta_{step}一般取值\frac{360}{180 * N}N\geq{2},输出图像矩阵宽w\approx\frac{r_{max}-r_{min}}{r_{step}}+1,高h\approx\frac{\theta_{max}-\theta_{max}}{\theta_{step}}+1

实现

首先了解以下Numpy的tile(a, (m, n))函数,此函数返回一个m×n个a平铺成的矩阵,垂直方向m个,水平方向n个,如:

a = np.array([[1, 2], [3, 4]])
print(np.tile(a, (2, 3)))

输出

array([[1, 2, 1, 2, 1, 2],
       [3, 4, 3, 4, 3, 4],
       [1, 2, 1, 2, 1, 2],
       [1, 2, 1, 2, 1, 2]])

对C++来说,OpenCV提供函数repeat(const Mat& src, int ny, int nx)实现类似的功能。

下面的代码是用Python实现极坐标变换,使用的是最近邻插值法,也可以使用其他插值方法:

def polar(I, center, r, theta=(0, 360), rstep=1.0, thetastep=360.0/(180*8)):
    """
    :param I: 输入的图像矩阵
    :param center: 极坐标变换中心
    :param r: 二元元组,代表最小和最大距离
    :param theta: 二元元组,角度范围
    :param rstep: r的变换步长
    :param thetastep: theta的变换步长
    :return 输出图像矩阵
    """
    # 得到距离的最小最大范围
    r_min, r_max = r
    # 角度最小最大范围
    theta_min, theta_max = theta
    # 输出图像的高、宽
    h = int((r_max - r_min)/rstep) + 1
    w = int((theta_max - theta_min)/thetastep) + 1
    O = 125 * np.ones((h, w), I.dtype)
    # 极坐标变换
    # linspace(start, stop, num=50)产生从start到stop,数量为num的等差数列
    r = np.linspace(r_min, r_max, h)
    r = np.tile(r, (w, 1))
    # 转置
    r = np.transpose(r)
    theta = np.linspace(theta_min, theta_max, w)
    theta = np.tile(theta, (h, 1))
    x, y = cv2.polarToCart(r, theta, angleInDegrees=True)
    # 最近邻插值法
    for i in range(h):
        for j in range(w):
            # round()按照指定精度四舍五入
            px = int(round(x[i][j])+cx)
            py = int(round(y[i][j])+cy)
            if (px>=0 and py <= w-1) and (py >=0 and py <= h-1):
                O[i][j] = I[py][px]
    return O

OpenCV3.X以上提供线性极坐标函数linearPolar(src, dst, Point2f center, double maxRadius, int flags);实现了我们上面的函数效果,其中center是变换中心,maxRadius是最大距离,flags是插值算法,值和函数resizewarpAffine一样。需要注意的是,linearPolar函数生成的极坐标,\theta在垂直方向上,r在水平方向上,和之前讨论的相反,旋转90°后得到的结果类似。

对数极坐标函数

OpenCV3.X中提供函数logPolar(src, dst, center, M, flags),其中M是系数,值大一些效果好;flags等于WARP_FILL_OUTLIERS代表笛卡尔坐标向对数坐标变换,等于WARP_INVERSE_MAP代表对数坐标向笛卡尔坐标变换。

将笛卡尔坐标转换为对数坐标的公式为:
\bar{r} = M * log \sqrt{(x-x_0)^2 + (y-y_0)^2}\\ \theta = 2 \pi + arctan2(y- y_0 , x-x_0 ),y - y_0 <= 0\\arctan2( y-y_0, x-x_0), y-y_0>0
反过来:
x= x_0 - exp(\frac{\bar{r}}{M}) cos \theta , y=y_0 - exp(\frac{\bar{r}}{M})) sin \theta
通过对M值的修改可以发现M值越大,在水平方向得到的信息越多。

参考

《OpenCV算法精解——基于Python和C++》(张平)第三章


一沙一世界,一花一天堂。君掌盛无边,刹那成永恒。