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

2019-05-04 本文已影响0人 Aengus_Sun

此篇文章是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_x$ 、 $s_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]$ 表示不大于 $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_NEAREST、INTE_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=\left\{\begin{array}22\pi+arctan2(y-y_0,x-x_0),y-y_0<=0\\arctan2(y-y_0,x-x_0),y-y_0>0 \end{array}\right.$

以变换中心为圆心的同一个圆上的点，在极坐标系 $\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}<=1$ 代表步长， $\theta_{step}$ 一般取值 $\frac{360}{180*N}$ ， $N>=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是插值算法，值和函数resize、warpAffine一样。需要注意的是，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=\left\{\begin{array}22 \pi+ arctan2(y-y_0,x-x_0),y-y_0<=0\\arctan2(y-y_0,x-x_0),y-y_0>0 \end{array}\right.$

反过来：

$x=x_0-exp(\frac{\bar{r}}{M})cos\theta,y=y_0-exp(\frac{\bar{r}}{M}))sin\theta$

通过对M值的修改可以发现M值越大，在水平方向得到的信息越多。