OpenCV算法学习笔记之几何变换
in OpenCV with 1 comment

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

in OpenCV with 1 comment

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

仿射变换

二维空间的仿射变换公式为:

$$ \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,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=\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或float63,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是插值算法,值和函数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=\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值越大,在水平方向得到的信息越多。

参考

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

Responses
  1. Aengus

    test

    Reply