计算视觉——图像梯度、图像边缘、几何特征、检测与提取在我很多期博客中，如基于内容的图像缩放、图像手绘/素描风格转换、图像

写在前面

在我很多期博客中，如基于内容的图像缩放、图像手绘/素描风格转换、图像的局部特征信息及全景图像拼接，都有涉及到图像梯度的概念，因为梯度确实是一种具有高信息熵且易于提取的信息，在图像处理、计算机视觉、机器视觉领域研究中应用广泛，我在很多完全不同的图像研究方向的文献中见到过图像梯度的身影。可以说，绝大多数图像相关的研究都可以利用图像的梯度去实现特定效果的优化，只是一些情况下该优化的作用很小。但是，这并不影响图像梯度的价值。而图像边缘、图像几何特征更是当前计算机视觉、模式识别研究的重点关注对象，图像边缘信息对行人检测、人脸识别等人工智能主流研究方向有着不可替代的作用。

在本期博客中，我们将探讨如下问题：

图像梯度与图像边缘
图像几何检测原理及其实现

Let's begin!!

1.图像梯度与图像边缘

我们来看一下维基百科对于Image Gradient(图像梯度)的解释。

"An image gradient is a directional change in the intensity or color in an image. The gradient of the image is one of the fundamental building blocks in image processing. For example, the Canny edge detector uses image gradient for edge detection. In graphics software for digital image editing, the term gradient or color gradient is also used for a gradual blend of color which can be considered as an even gradation from low to high values, as used from white to black in the images to the right. Mathematically, the gradient of a two-variable function (here the image intensity function) at each image point is a 2D vector with the components given by the derivatives in the horizontal and vertical directions. At each image point, the gradient vector points in the direction of largest possible intensity increase, and the length of the gradient vector corresponds to the rate of change in that direction."

"图像梯度是指图像强度或颜色的方向变化。图像梯度是图像处理的基础框架之一。例如，Canny边缘检测器使用图像梯度进行边缘检测。在用于数字图像编辑的图形软件中，梯度或颜色梯度一词也用于表示颜色的逐渐混合，可以认为是从低到高值的均匀渐变，如下图中从白色到黑色所使用的渐变。数学上，一个双变量函数(这里是图像强度函数)在每个图像点处的梯度是一个二维向量，其分量由水平和垂直方向上的导数给出。在每个图像点上，梯度向量指向可能最大强度增加的方向，梯度向量的长度对应于该方向的变化率。"

上述内容非常严谨地定义了图像和视觉中的图像梯度，但是仍然需要指出，蓝字部分的内容存在一些难以理解的部分，如颜色的逐渐混合，其实这里指的是Photoshop（或其他的图形编辑软件）中渐变的效果，当然，这些渐变也属于梯度。

维基百科用了较通俗的语言去描述图像梯度，相比之下，百度百科对图像梯度的解释就显得更专业，更容易理解，如下：
"图像梯度可以把图像看作二维离散函数，图像梯度其实就是这个二维离散函数的求导：
图像梯度：

其中， $I(i,j)$ 为灰度图像 $I$ 在第 $i$ 行第 $j$ 个像素上的像素值。
上面说的是简单的梯度定义，其实还有更多更复杂的梯度公式。"

我更喜欢百度百科更接近本质的解释，简单来说，图像梯度就是图像在不同方向上的变化趋势，可以将图像视为只存在水平和垂直两方向存在梯度变化，进行临近像素的梯度提取。当然，这么做会存在一定的问题，那就是这样求解得到的梯度信息对噪声敏感，即对噪声鲁棒性差。

接下来，将以实际的图像为例，更进一步的解释 1.1 图像梯度、图像边缘及梯度提取算子，并介绍目前最广泛应用的 1.2 canny边缘检测算法原理及实现。

1.1图像梯度、图像边缘及梯度提取算子

说起图像梯度，其实我们得好好谈一谈图像边缘，从许多曾经的视觉研究实验中，我们可以发现，边缘是很特殊的一种信息，因为人眼识别物体，靠的正是边缘，而非颜色、纹理。在绘画中我们可以通过线条去绘制任何一样物体，但若该手绘物体缺失了部分边缘，我们就难以识别。但相反的，我们对缺失了纹理、色彩的物体，却依然具备强大的辨别能力。

Walther等人在2011年进行了这样一个视觉研究，让视觉受试者观看海滩、城市街道、森林、高速公路、山脉和办公室的照片和线条图时收集功能磁共振成像数据。尽管在场景统计数据上有明显的差异，但研究者们还是能够通过PPA(the parahippocampal place area,海马旁区和RSC(the ret-rosplenial cortex,脾叶皮层)解码初级视觉皮层中线条图的fMRI(Functional magnetic resonance imaging,功能磁共振成像)数据和彩色照片的活动数据的场景类别。更值得注意的是，在PPA和RSC中，从线条图识别的错误模式与从彩色照片识别的错误模式非常相似。这些数据表明，在这些区域，用于区分场景类别的信息对于线条图和照片来说是相似的。

而我们在谈论图像梯度的时候，很多情况下关注的是图像边缘。因此，图像的边缘检测的目的，即为确定图像突变(不连续)的地方

显然，图像内高信息熵信息隐藏在边缘信息当中
边缘传递的信息多余像素点本身

边缘的由来主要源于：

表面的不连续

深度(距离)的不连续

表面颜色的不连续

光照的不连续

已经介绍完梯度与边缘，我们现在就先来看看梯度的计算方式。在先前度娘百科(百度百科)对于图像梯度已经有了一个数理上的定义，当然，对于图像梯度的表示和计算方式有很多，对于一幅给定的数字图像，我们将其视作二维离散函数 $f(x,y)$ ，那么常见的梯度表征和计算方式有：

梯度向量

\triangledown f(x,y)=\begin{bmatrix} \frac{\partial f(x,y)}{\partial x} \\ \frac{\partial f(x,y)}{\partial y} \end{bmatrix}=\begin{bmatrix} f_x \\ f_y \end{bmatrix}

梯度幅值

\begin{vmatrix} \triangledown f(x,y) \end{vmatrix}=\sqrt{f^2_x+f^2_y}

梯度方向

而在实际编程运算时，我们通常采用算子通过卷积获得一幅图像的梯度信息，如果不明白何为图像中的卷积运算，可以看看下图。

实际上，上图中的原图像大小为 $5\times 5$ （左侧），而卷积核为 $3\times 3$ （中间），因此得到的卷积结果是 $3\times 3$ 的（右侧）。原图像大小为 $m\times m$ (m为奇数)的图像，经过 $k\times k$ (k为奇数)的卷积核卷积后，得到的卷积结果图像大小为 $(m-k+1)\times (m-k+1)$ 。

以最简单的Prewitt算子为例，利用算子卷积的图像梯度结果如下：

不难发现，在梯度图像中，值越高的部分，越接近人眼视觉所观察到的边缘。因此，我们可以用梯度更直接的解释边缘：边缘是图像能量急剧变化的地方，是图像梯度高的地方。

但这也存在一个问题，图像的梯度是对邻域像素值变化敏感的，而图像内的噪声则会导致部分非边缘区域同样具有较高的梯度，这会导致边缘检测的噪声鲁棒性差。

为了提高噪声鲁棒性，Canny在1986年提出了一种边缘检测算法，该算法是边缘检测中最经典，最具影响力的方法，至今还有非常广泛的应用，是计算机视觉领域的基础框架组成之一，论文被引超20000余次。

1.2 Canny边缘检测算法原理

Canny边缘检测算法可以分为四步：

噪声抑制
计算图像梯度及方向
非极大值抑制
双阈值检测
边缘连接

对于噪声抑制和梯度及梯度方向的计算，Canny边缘检测算法并没有进行创新，只是做了单纯的高斯核卷积平滑，以此进行噪声抑制，并通过差分运算得到梯度及梯度方向。Canny边缘检测突破性的创新是集中在后三步上的。

非极大值抑制

非极大值抑制的思想在于：当图像的梯度到达局部(特定大小的邻域)最大值的时候，我们认为它是边缘；对于非极大值的点，即使它超过了我们所定义的“边缘梯度阈值”，仍不认为它是边缘。并且，预定义了八个可能的方向（左上、上、右上、右、右下、下、左下、左），对于该方向上的梯度值，将非极大值归零。

具体操作：1.筛选幅值大于阈值的像素点，设为待验证点；2.遍历待验证点，若为局部最大值，保留，反之去归零；3.保留下来的点就是非极大值抑制后的边缘点。

双阈值检测

说到双阈值检测，我们得先说采用单一阈值的缺陷，如下图所示。

显而易见，采用单一阈值会有两种问题随之产生：1.阈值过高，会导致部分梯度较低的边缘被当成噪声剔除，导致边缘不连续；2.阈值过低，会有过多不应该作为边缘的梯度信息被误判为边缘而包含进边缘图。

而双阈值解决了单阈值的上述问题。双阈值定义了两种阈值，低阈值与高阈值：若梯度小于低阈值，则不是边缘；若梯度大于高阈值，则为强边缘；若梯度在高低阈值间，则为弱边缘。

至此，就是双阈值检测的部分，乍一看仿佛与单一阈值没有区别，只要高于低阈值就能作为弱边缘，那这么一来好像跟单阈值就没什么区别了。诚然，如果缺少了“边缘连接”这下一步骤，双阈值和单阈值的本质确实是一样的。

边缘连接

在开始解释“边缘连接”这一步骤前，我们不妨来仔细观察一下“单阈值”与“双阈值”的图像，这对于“边缘连接”原理的理解非常有帮助。

在单阈值图像中，存在部分“离群”的线条，其中有一部分其实并非真实的边缘，而是被误判的高梯度位置；其中也有一部分不连通的真实边缘，这又导致了边缘的不连续。

而在双阈值图像中，较亮的线条为强边缘，稍暗一些的线条为弱边缘。不难发现，一条不与强边缘相连的弱边缘，基本都是被误判为边缘的区域。相反，与强边缘相连的弱边缘，则囊括了绝大多数单阈值所没有判定包含的“真实边缘”。

那么现在，我们可以明白边缘连接的原理了：利用高阈值来确定主要边缘轮廓（强边缘），并利用与强边缘连接的低阈值所得弱边缘来填补缺失边缘，并防止引入噪声（不与强边缘相连的弱边缘就归零）。

由此，完整地解释了Canny边缘检测算法的原理，我们进入代码部分。

1.3 代码实战部分

### edge.py
import numpy as np
import math
def conv(image, kernel):
    """卷积的一个实现.

    对于任意一个像素点，该本版采用了点积运算以及 np.sum()来进行快速的加权求和

    Args:
        图像: 尺寸为(Hi, Wi)的numpy数组.
        卷积核(kernel): 尺寸为(Hk, Wk)的numpy数组.

    Returns:
        out: 尺寸为(Hi, Wi)的numpy数组.
    """
    Hi, Wi = image.shape
    Hk, Wk = kernel.shape
    out = np.zeros((Hi, Wi))

    # 对于本次作业, 我们将使用边界值来对图像进行填充.
    # 这是因为，如果使用0进行填充，会使得在图像边界上的导数非常大,
    # 但通常而言，我们希望忽略图像边界上的边缘（图像边界永远不会被认为是边缘）.
    pad_width0 = Hk // 2
    pad_width1 = Wk // 2
    pad_width = ((pad_width0,pad_width0),(pad_width1,pad_width1))
    padded = np.pad(image, pad_width, mode='edge')

    ### YOUR CODE HERE
    for i in range(Hi):
        for j in range(Wi):
            out[i][j] = (padded[i:i + Hk,j:j + Wk] * kernel).sum()
    ### END YOUR CODE

    return out

def gaussian_kernel(size, sigma):
    """ 生成高斯卷积核.

    该函数按照高斯核的公式，生成了一个卷积核矩阵。

    提示:
    - 利用 np.pi 以及 np.exp 当计算 pi 以及 exp()函数.

    参数:
        size: 一个整数，表示输出的kernel的尺寸.
        sigma: 一个float，对应于高斯公式中的sigma，用来控制权重的分配.

    返回值:
        kernel: 尺寸为(size, size)的numpy数组.
    """

    kernel = np.zeros((size, size))

    ### YOUR CODE HERE
    k=(size-1)/2
    for i in range(size):
        for j in range(size):
            kernel[i][j]=1/(2*np.pi*sigma**2)*np.exp(((i-k)**2+(j-k)**2)/(-2*sigma**2))
            #kenel[i][j]=1/(2*np.pi*power(sigma,2))*np.exp(((power(i-k,2)+(power(j-k,2))/(-2*power(sigma,2)))
    ### END YOUR CODE

    return kernel

def partial_x(img):
    """ 计算输入图像水平方向的偏导.

    提示:
        - 你可以利用你在前面完成的conv函数.

    参数:
        img: 尺寸为(Hi, Wi)的numpy数组.
    输出:
        out: 水平方向的梯度图像
    """
    
    Hi, Wi = img.shape
    out = None
    padded = np.pad(img, 1, mode='edge')
    padded1 = np.zeros(padded.shape)
    ### YOUR CODE HERE
    for i in range(1,Hi+1):
        for j in range(1,Wi+1):
            padded1[i][j]=(padded[i][j+1]-padded[i][j-1])/2
    ### END YOUR CODE
    out=padded1[1:Hi+1,1:Wi+1]
    return out

def partial_y(img):
    """ 计算输入图像竖直方向的偏导.

    提示:
        - 你可以利用你在前面完成的conv函数/或者用前面刚开发的partial_x函数.

    参数:
        img: 尺寸为(Hi, Wi)的numpy数组.
    输出:
        out: 竖直方向的梯度图像
    """
    Hi, Wi = img.shape
    out = None
    padded = np.pad(img, 1, mode='edge')
    padded1 = np.zeros(padded.shape)
    ### YOUR CODE HERE
    for i in range(1,Hi+1):
        for j in range(1,Wi+1):
            padded1[i][j]=(padded[i+1][j]-padded[i-1][j])/2
    ### END YOUR CODE
    out=padded1[1:Hi+1,1:Wi+1]
    return out

def gradient(img):
    """ 计算输入图像的梯度大小和方向.

    参数:
        img: 灰度图. 尺寸为 (H, W) 的Numpy数组.

    返回值:
        G: 输入图像的梯度值图像。它的每个像素点值都是该像素点的梯度值.
            尺寸为 (H, W) 的Numpy数组.
        theta: 输入图像的梯度方向图像(角度, 0 <= theta < 360)，每个像素点都代表该像素点的梯度方向.
            尺寸为 (H, W) 的Numpy数组.

    提示:
        - 可以使用 np.sqrt 以及 np.arctan2 来计算平方根以及反正切值
    """
    Hi, Wi = img.shape
    G = np.zeros(img.shape)
    theta = np.zeros(img.shape)
    outx=partial_x(img)
    outy=partial_y(img)
    ### YOUR CODE HERE
    for i in range(Hi):
        for j in range(Wi):
            G[i][j]=np.sqrt(outx[i][j]**2+outy[i][j]**2)
    for i in range(Hi):
        for j in range(Wi):
            theta[i][j]=np.arctan2(outy[i][j],outx[i][j])* 180 / np.pi
            if theta[i][j]<0:
                theta[i][j]=theta[i][j]+360
            if theta[i][j]>360:
                theta[i][j]=theta[i][j]-360
    ### END YOUR CODE
    
    return G, theta


def non_maximum_suppression(G, theta):
    """ 实现非极大值抑制

    对于任意一个像素点，该函数都完成了在梯度方向上(theta)的非极大值抑制

    参数:
        G: 梯度幅值图像，尺寸为 (H, W)的Numpy数组.
        theta: 梯度方向图像，尺寸为 (H, W)的Numpy数组.

    Returns:
        out: 非极大值抑制以后的图像。如果像素点不是局部极大值点，则置0，否则保留原梯度幅值.
    """
    padded = np.pad(G, 1, mode='edge')
    H, W = G.shape
    out = np.zeros((H, W))

    # 将梯度方向round到最近的45°上来
    theta = np.floor((theta + 22.5) / 45) * 45
    ### BEGIN YOUR CODE
    for i in range(1,H+1):
        for j in range(1,W+1):
            thetaij=theta[i-1][j-1]
            if thetaij==45:
                if padded[i][j]>padded[i+1][j+1] and padded[i][j]>padded[i-1][j-1]:
                    out[i-1][j-1]=padded[i][j]
                else:
                    out[i-1][j-1]=0
            if thetaij==225:
                if padded[i][j]>padded[i+1][j+1] and padded[i][j]>padded[i-1][j-1]:
                    out[i-1][j-1]=padded[i][j]
                else:
                    out[i-1][j-1]=0
                                        
            if thetaij==90:
                if padded[i][j]>padded[i+1][j] and padded[i][j]>padded[i-1][j]:
                    out[i-1][j-1]=padded[i][j]
                else:
                    out[i-1][j-1]=0
            if thetaij==270:
                if padded[i][j]>padded[i+1][j] and padded[i][j]>padded[i-1][j]:
                    out[i-1][j-1]=padded[i][j]
                else:
                    out[i-1][j-1]=0
                                 
            if thetaij==135:
                if padded[i][j]>padded[i-1][j+1] and padded[i][j]>padded[i+1][j-1]:
                    out[i-1][j-1]=padded[i][j]
                else:
                    out[i-1][j-1]=0
            if thetaij==315:
                if padded[i][j]>padded[i-1][j+1] and padded[i][j]>padded[i+1][j-1]:
                    out[i-1][j-1]=padded[i][j]
                else:
                    out[i-1][j-1]=0
                    
            if thetaij==0:
                if padded[i][j]>padded[i][j-1] and padded[i][j]>padded[i][j+1]:
                    out[i-1][j-1]=padded[i][j]
                else:
                    out[i-1][j-1]=0
            if thetaij==180:
                if padded[i][j]>padded[i][j-1] and padded[i][j]>padded[i][j+1]:
                    out[i-1][j-1]=padded[i][j]
                else:
                    out[i-1][j-1]=0
            if thetaij==360:
                if padded[i][j]>padded[i][j-1] and padded[i][j]>padded[i][j+1]:
                    out[i-1][j-1]=padded[i][j]
                else:
                    out[i-1][j-1]=0
    ### END YOUR CODE

    return out


def double_thresholding(img, high, low):
    """
    参数:
        img: 非极大值抑制以后的梯度幅值图像.
        high: 对于强边缘(strong edge)而言的高阈值.
        low: 对于弱边缘(strong edge)而言的高阈值.

    返回值:
        strong_edges: 布尔类型的数组，代表强边缘。
            强边缘指的是梯度的幅值大于高阈值的像素点
        弱边缘: 布尔类型的数组，代表弱边缘。
            强边缘指的是梯度的幅值小于或者等于高阈值，并且大于低阈值的像素点。
    """

    strong_edges = np.zeros(img.shape, dtype=np.bool)
    weak_edges = np.zeros(img.shape, dtype=np.bool)
    H, W = img.shape
    ### YOUR CODE HERE
    for i in range(H):
        for j in range(W):
            if img[i][j]>high:
                strong_edges[i][j]=img[i][j]
            if img[i][j]<=high and img[i][j]>low:
                weak_edges[i][j]=img[i][j]
                
    ### END YOUR CODE

    return strong_edges, weak_edges


def get_neighbors(y, x, H, W):
    """ 返回坐标为 (y, x) 的像素点点的邻居(neighbor).

    对于一幅尺寸为(H, W)的图像，返回像素点 (y, x) 的邻居的所有有效索引.
    一个有效的索引 (i, j) 应该满足:
        1. i >= 0 and i < H
        2. j >= 0 and j < W
        3. (i, j) != (y, x)

    参数:
        y, x: 像素点的位置
        H, W: 图像的尺寸
    返回值:
        neighbors: 该像素点的邻居的索引值[(i, j)]所组成的list .
    """
    neighbors = []

    for i in (y-1, y, y+1):
        for j in (x-1, x, x+1):
            if i >= 0 and i < H and j >= 0 and j < W:
                if (i == y and j == x):
                    continue
                neighbors.append((i, j))

    return neighbors

def link_edges(strong_edges, weak_edges):
    """ 找出与真实的边缘相连接的弱边缘，并将它们连接到一起.
    
    对于每个强边缘点，它们都是真实的边缘。
    我们需要遍历每个真实的边缘点，然后在弱边缘中找出与之相邻的像素点，并把他们连接起来。
    在这里，我们认为如果像素点(a, b)与像素点(c, d)相连，只要(a, b)位于(c, d)的八邻域内。

    Args:
        strong_edges: 尺寸为 (H, W)的二值图像.
        weak_edges: 尺寸为 (H, W)的二值图像.
    
    Returns:
        edges: 尺寸为 (H, W)的二值图像.
        
    提示：
        弱边缘一旦与强边缘相连，那它就是真实的边缘，与强边缘的地位一样。
        这句话的意思是，一旦一个弱边缘像素点被检测成为了一个真实的边缘点，所有与它相连的其它弱边缘也应该被归为真实的边缘。
        所以在编程的时候，你只遍历一遍强边缘是不够的，因为此时可能有新的弱边缘点被标记为真实的边缘点。
    """

    H, W = strong_edges.shape
    indices = np.stack(np.nonzero(strong_edges)).T
    edges = np.zeros((H, W), dtype=np.bool)

    # 生成新的拷贝
    weak_edges = np.copy(weak_edges)
    edges = np.copy(strong_edges)
    g=0;
    padded_weak = np.pad(weak_edges, 1, mode='edge')
    edges1=np.pad(strong_edges, 1, mode='edge')
    ### YOUR CODE HERE
    while(g!=-1):
        g=0;
        k=0;
        for i in range(H):              
            for j in range(W):
                #if strong_edges[i][j]!=0:
                if edges1[i+1][j+1]!=0:
                    if padded_weak[i][j]!=0 and edges1[i][j]!=padded_weak[i][j]:
                        edges1[i][j]=padded_weak[i][j]
                        k=k+1
                    if padded_weak[i][j+1]!=0 and edges1[i][j+1]!=padded_weak[i][j+1]:
                        edges1[i][j+1]=padded_weak[i][j+1]
                        k=k+1
                    if padded_weak[i][j+2]!=0 and edges1[i][j+2]!=padded_weak[i][j+2]:
                        edges1[i][j+2]=padded_weak[i][j+2]
                        k=k+1
                    
                    if padded_weak[i+1][j]!=0 and edges1[i+1][j]!=padded_weak[i+1][j]:
                        edges1[i+1][j]=padded_weak[i+1][j]
                        k=k+1
                    if padded_weak[i+1][j+2]!=0 and edges1[i+1][j+2]!=padded_weak[i+1][j+2]:
                        edges1[i+1][j+2]=padded_weak[i+1][j+2]
                        k=k+1
                        
                    if padded_weak[i+2][j]!=0 and edges1[i+2][j]!=padded_weak[i+2][j]:
                        edges1[i+2][j]=padded_weak[i+2][j]
                        k=k+1
                    if padded_weak[i+2][j+1]!=0 and edges1[i+2][j+1]!=padded_weak[i+2][j+1]:
                        edges1[i+2][j+1]=padded_weak[i+2][j+1]
                        k=k+1
                    if padded_weak[i+2][j+2]!=0 and edges1[i+2][j+2]!=padded_weak[i+2][j+2]:
                        edges1[i+2][j+2]=padded_weak[i+2][j+2]
                        k=k+1
        if k==0:
            g=-1;
        else:
            strong_egdes=np.copy(edges1[1:H+1,1:W+1])
    
    edges=edges1[1:H+1,1:W+1]           
                        
                    
            
    ### END YOUR CODE

    return edges

def canny(img, kernel_size=5, sigma=1.4, high=20, low=15):
    """ 将上面所完成的函数组合到一起，完成canny边缘检测器.

    Args:
        img: 输入图像
        kernel_size: int， 表示kernel的大小
        sigma: float, 用来计算kernel.
        high: 为强边缘而设的高阈值.
        low: 为弱边缘而设的低阈值.
    Returns:
        edge: 输出的边缘图像
    """
    ### YOUR CODE HERE
    
    kernel = gaussian_kernel(kernel_size, sigma)
    smoothed = conv(img, kernel)
    G,theta=gradient(smoothed)
    nms = non_maximum_suppression(G, theta)
    strong_edges, weak_edges = double_thresholding(nms, high, low)
    edge = link_edges(strong_edges, weak_edges)
    
    ### END YOUR CODE

    return edge

1.4 实施例

初始化

# 初始化配置
import numpy as np
import matplotlib.pyplot as plt
from time import time
from skimage import io


%matplotlib inline
plt.rcParams['figure.figsize'] = (15.0, 12.0) # 设置默认尺寸
plt.rcParams['image.interpolation'] = 'nearest'
plt.rcParams['image.cmap'] = 'gray'

# 自动重装外部模块
%load_ext autoreload
%autoreload 2

首先利用高斯滤波器来对图片进行平滑处理。高斯核大小为 $(2k+1)\times (2k+1)$ 的高斯滤波器，该滤波器可以由以下方程给出

h_{ij}=\frac{1}{2\pi \sigma^2}\exp\bigg(-\frac{(i-k)^2+(ji-k)^2}{2\sigma^2}\bigg),0\le i,j\le 2k+1

from edge import conv, gaussian_kernel

# 定义一个 3x3 的高斯 kernel，并将其sigma值设为 1
kernel = gaussian_kernel(3, 1)
kernel_test = np.array(
    [[ 0.05854983, 0.09653235, 0.05854983],
     [ 0.09653235, 0.15915494, 0.09653235],
     [ 0.05854983, 0.09653235, 0.05854983]]
)
print(kernel)
# 检测生成的高斯kernel是否正确
if not np.allclose(kernel, kernel_test):
    print('Incorrect values! Please check your implementation.')

输出结果：

[[0.05854983 0.09653235 0.05854983]
 [0.09653235 0.15915494 0.09653235]
 [0.05854983 0.09653235 0.05854983]]

找一幅图像进行测试。

# 用不同的尺寸以及sigma值来进行测试
kernel_size = 5
sigma = 1.4

# 载入图片
img = io.imread('iguana.png', as_grey=True)

# 生成高斯kernel
kernel = gaussian_kernel(kernel_size, sigma)

# 利用kernel来对图片进行平滑
smoothed = conv(img, kernel)

plt.subplot(1,2,1)
plt.imshow(img)
plt.title('Original image')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(smoothed)
plt.title('Smoothed image')
plt.axis('off')

plt.show()

这里有个小问答：改变kernel_size 以及 sigma分别具有什么效果？

改变核的大小：

会导致中心像素点受临域像素值的影响更大，从到导致高斯滤波的模糊效果更佳明显。

改变sigma的大小：

若减小，会使高斯滤波核中心权重占比变大，中心位置临域的权重占比都减小，若sigma趋于0，则高斯滤波基本没有效果；若增大，会使高斯滤波核的中心位置权重占比减小，中心位置领域的权重占比增大，若sigma趋于无穷，则高斯滤波等于均值滤波。

我们回到实施例部分，找出求取图像差分所对应的kernel $D_x$ 和 $D_y$

from edge import partial_x, partial_y

# 测试案例
I = np.array(
    [[0, 0, 0],
     [0, 1, 0],
     [0, 0, 0]]
)

# 希望的输出结果
I_x_test = np.array(
    [[ 0, 0, 0],
     [ 0.5, 0, -0.5],
     [ 0, 0, 0]]
)

I_y_test = np.array(
    [[ 0, 0.5, 0],
     [ 0, 0, 0],
     [ 0, -0.5, 0]]
)

# 计算梯度
I_x = partial_x(I)
I_y = partial_y(I)
print(I_x)
print(I_y)
# 确定partial_x and partial_y是否编写正确
if not np.all(I_x == I_x_test):
    print('partial_x incorrect')
    
if not np.all(I_y == I_y_test):
    print('partial_y incorrect')

输出结果：

[[ 0.   0.   0. ]
 [ 0.5  0.  -0.5]
 [ 0.   0.   0. ]]
[[ 0.   0.5  0. ]
 [ 0.   0.   0. ]
 [ 0.  -0.5  0. ]]

计算x、y方向上的梯度

# 计算平滑后的图像的差分
Gx = partial_x(smoothed)
Gy = partial_y(smoothed)

plt.subplot(1,2,1)
plt.imshow(Gx)
plt.title('Derivative in x direction')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(Gy)
plt.title('Derivative in y direction')
plt.axis('off')

plt.show()

结果如下

现在，让我们用两个方向的差分来计算图片的梯度大小以及方向，完成edge.py所定义的梯度和梯度方向的求解。

from edge import gradient

G, theta = gradient(smoothed)

if not np.all(G >= 0):
    print('Magnitude of gradients should be non-negative.')
    
if not np.all((theta >= 0) * (theta < 360)):
    print('Direction of gradients should be in range 0 <= theta < 360')

plt.imshow(G)
plt.title('Gradient magnitude')
plt.axis('off')
plt.show()

可以发现，找出来的梯度图像比较粗大并且模糊，并不适合直接用来表示的边缘。此时我们可以利用非极大值抑制的方法来找出“尖锐”的边缘。简单来说，这可以通过保留（梯度）图像中的局部极大值点并且抛弃其他点来进行。

from edge import non_maximum_suppression

# 测试例
g = np.array(
    [[0.4, 0.5, 0.6],
     [0.3, 0.5, 0.7],
     [0.4, 0.5, 0.6]]
)

# 输出非极大值抑制的结果
# 改变梯度方向：通过四个方向来测试，即0,45,90,135。你可以将输出结果与自己笔算的结果做对比。
for angle in range(0, 180, 45):
    print('Thetas:', angle)
    t = np.ones((3, 3)) * angle # Initialize theta
    #print(t)
    print(non_maximum_suppression(g, t))

输出结果：

Thetas: 0
[[0. 0. 0.]
 [0. 0. 0.]
 [0. 0. 0.]]
Thetas: 45
[[0.  0.  0. ]
 [0.  0.  0.7]
 [0.  0.  0. ]]
Thetas: 90
[[0.  0.  0. ]
 [0.  0.  0.7]
 [0.  0.  0. ]]
Thetas: 135
[[0.  0.  0. ]
 [0.  0.  0.7]
 [0.  0.  0. ]]

nms = non_maximum_suppression(G, theta)

plt.imshow(nms)
plt.title('Non-maximum suppressed')
plt.axis('off')
plt.show()

在经过非极大值抑制后，仍然存在许多像素点。这些像素点有些是边缘，但有些是由噪声或者颜色变化（例如粗糙的表面）所导致的。要消除这部分影响，最简单的方式是增加一个阈值判定，只允许相应强度大于阈值的点被标记成边缘。Canny边缘检测算法采用了双阈值算法。大于高阈值的像素点被标记为强边缘，小于低阈值的像素点可以认为是非边缘并且被移除，而在高低阈值之间的点被标记为弱边缘。

from edge import double_thresholding

low_threshold = 0.02
high_threshold = 0.03

strong_edges, weak_edges = double_thresholding(nms, high_threshold, low_threshold)
assert(np.sum(strong_edges & weak_edges) == 0)

edges=strong_edges * 1.0 + weak_edges * 0.5

plt.subplot(1,2,1)
plt.imshow(strong_edges)
plt.title('Strong Edges')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(edges)
plt.title('Strong+Weak Edges')
plt.axis('off')

plt.show()

通常而言，强边缘被认为是“肯定的边缘”, 因此可以直接认为该像素点就是边缘. 而弱边缘则需要进一步判断。如果弱边缘点与边缘相连接，则认为该像素点是边缘。这背后的逻辑是这样的：噪声或者颜色变化点不太可能产生强边缘（如果设定好恰当对的阈值范围），因此强边缘只能由源图像中的边缘产生。而弱边缘可由边缘产生，也可以由噪声和颜色变化产生。而由噪声和颜色变化产生的弱边缘点，一般来说都均有分布在整幅图片上，只有一小部分是与强边缘相连接的，大部分与强边缘相连接的点都是真实的边缘像素点。

先用一个简单测试例来展示强弱边缘连接及判别的过程。

from edge import get_neighbors, link_edges

test_strong = np.array(
    [[1, 0, 0, 0],
     [0, 0, 0, 0],
     [0, 0, 0, 0],
     [0, 0, 0, 1]]
)

test_weak = np.array(
    [[0, 0, 0, 1],
     [0, 1, 0, 0],
     [1, 0, 0, 0],
     [0, 0, 1, 0]]
)

test_linked = link_edges(test_strong, test_weak)

plt.subplot(1, 3, 1)
plt.imshow(test_strong)
plt.title('Strong edges')

plt.subplot(1, 3, 2)
plt.imshow(test_weak)
plt.title('Weak edges')

plt.subplot(1, 3, 3)
plt.imshow(test_linked)
plt.title('Linked edges')
plt.show()

接下来，我们在图像中进行边缘连接

edges = link_edges(strong_edges, weak_edges)

plt.imshow(edges)
plt.axis('off')
plt.show()

测试一组不同超参数选取下的Canny边缘检测结果

from edge import canny

# 载入图像
img = io.imread('iguana.png', as_grey=True)

# 运行canny边缘检测器
edges = canny(img, kernel_size=5, sigma=1.4, high=0.03, low=0.02)
print (edges.shape)

plt.subplot(1, 3, 1)
plt.imshow(edges)
plt.axis('off')
plt.title('Your result')

2 图像几何特征检测

图像的特征主要有图像的颜色特征、纹理特征、几何特征和空间关系特征。其中，几何特征也称形状特征。在物体识别、检测领域有非常重要的作用。Hough变换是其中非常经典的方法。

2.1 Hough变换简介

Hough变换是在1962年由Hough所发明，随后在1972年Duda第一次用它检测图像中的直线。该变换的目的是为了寻找图像中的直线结构，但需要注意的是，Hough变换同样可以检测圆等其他几何结构，只要其参数方程是已知的。

优点

概念简单
操作简单
被遮挡的图形仍具备几何检测的交过
对于有参数的几何特征均有效

缺点

对于参数较多的图形，计算复杂度高
每次只能检测一种几何特征
以直线检测为例，不能确定线段的长度与位置，不能区分两段共线的线段

假设

假设我们现在已经通过边缘检测算法得到了图像完整的边缘结构，可以发现，某些像素点构成了图像的几何结构。

接下来，我们将以直线检测为例，说明Hough变换的原理。

2.2 Hough变换原理——以直线检测为例

对于任意一个像素点 $(x_i,y_i)$ ，可能存在许多直线都通过该点，且所有通过该点的直线均具有如下方程

一般而言，都具有两个位置参数 $(a,b)$ ，则该方程可以改为

由此，我们将 $x_i,y_i$ 视作参数， $a,b$ 视作变量，那么在xy空间并且经过点 $(x_i,y_i)$ 的直线，在ab空间都经过点 $(a_i,b_i)$

那么这是一条在 $ab$ 空间，由 $x$ 和 $y$ 作为参数的直线，因此，一个在 $xy$ 空间的点 $(x_1,y_1)$ 给定了一条在 $ab$ 空间的直线；而另一个点 $(x_2,y_2)$ 将会给出 $ab$ 空间的另一条直线。若 $(x_1,y_1)$ 和 $(x_2,y_21)$ 是共线的，那么它们在 $ab$ 空间的直线一定相交，交点坐标就是直线的参数 $(a_i,b_i)$

如果我们将 $xy$ 空间中直线上的点全部映射到 $ab$ 空间，会发现几个非常有趣的现象。

在 $xy$ 空间的两点 $(x_1,y_1)$ 和 $(x_2,y_2)$ 确定了一条直线，该两点对应与 $ab$ 空间的两条直线。
在 $ab$ 空间，这些直线会交于 $(a',b')$ ，该交点就是 $xy$ 空间上， $(x_1,y_1)$ 和 $(x_2,y_2)$ 所连成的直线的参数。
所有在 $xy$ 空间中 $(x_1,y_1)$ 与 $(x_2,y_2)$ 连线上的点，都对应于一条在 $ab$ 空间内的直线，并且这些直线都相交于同一点 $(a',b')$

基于上述特性，我们可以实现一种基于Hough变换的直线检测，步骤如下：

将参数空间 $ab$ 离散取值，分为一个个单元，该单元被称为累积单元(accumulator cells)
计算直线与单元相交次数
- 对于任意一组被检测为边缘的点对 $(x_1,y_1)$ 和 $(x_2,y_2)$ 检测它们在 $(a,b)$ 空间内的交点 $(a',b')$
- 将交点所在的单元的计数+1
- 若单元的计数超过某一设定的计数(votes,此处的votes也可以理解为阈值，但是解释成计数更合适)，那么该单元对应于 $xy$ 空间的一条直线

下图是一个Hough变换的结果例，右图显示了计数值前20的直线

由此，完整地解释了Hough变换在直线检测中的原理，在其中若对参数方程进行修改，就可以实现对其他几何形状的检测。

2.3 代码实战

### hough.py
import numpy as np
import math
def hough_transform(img):
    """ 

    利用下面的参数化方程:
        rho = x * cos(theta) + y * sin(theta)
    将点 (x,y) 转到Hough空间中类似于正弦函数( sine-like function)的一条曲线.

    参数:
        img: 尺寸为 (H, W)的二值图像.
        
    返回值:
        accumulator: 形状为 (m, n)的Numpy数组.
        rhos: 形状为 (m, )的Numpy数组.
        thetas: 形状为 (n, )的Numpy数组.
    """
    # 设置 rho 和 theta 的范围
    W, H = img.shape
    diag_len = int(np.ceil(np.sqrt(W * W + H * H)))
    rhos = np.linspace(-diag_len, diag_len, diag_len * 2.0 + 1)
    thetas = np.deg2rad(np.arange(-90.0, 90.0))
    # 计算一些会重复使用的值
    cos_t = np.cos(thetas)
    sin_t = np.sin(thetas)
    num_thetas = len(thetas)

    # 在Hough 空间中初始化 accumulator 
    accumulator = np.zeros((2 * diag_len + 1, num_thetas), dtype=np.uint64)
    ys, xs = np.nonzero(img)

    # 将图像中的每一个像素点(x,y)转到对应的Hough空间中
    # 找出每个thetas所对应的rhos
    # 并在accumulator的相应位置上加1
    ### YOUR CODE HERE
    for i in range(len(xs)):
        x=xs[i]
        y=ys[i]
        for j in range(num_thetas):
            rho=np.floor(x*cos_t[j]+y*sin_t[j])+diag_len
            
            accumulator[int(rho),j]+=1
    ### END YOUR CODE

    return accumulator, rhos, thetas

2.4 道路检测实施例

通常而言，车道都是细长的直线，并且具有较高的亮度。这种条件保证了它能被我们的算法较好地检测出来。运行下面代码，载入图片并检测图片中的边缘。

from edge import canny

# 载入图像
img = io.imread('road.jpg', as_grey=True)

# 利用Canny检测器进行边缘检测
edges = canny(img, kernel_size=5, sigma=1.4, high=0.03, low=0.02)

plt.subplot(211)
plt.imshow(img)
plt.axis('off')
plt.title('Input Image')

plt.subplot(212)
plt.imshow(edges)
plt.axis('off')
plt.title('Edges')
plt.show()

可以发现，Canny算法能够找出图片中的道路。但是，我们也可以看到这里面出现了很多我们不需要的物体的边缘。考虑到所抓取的图片的空间位置，我们清楚道路总是在图片的下半部分，因此可以利用这一点消除一部分不希望得到的边缘。下面的代码定义了一个二值化的模板(mask)，用来抓取ROI(Region of interest, 感兴趣区域)的边缘。

H, W = img.shape

# Generate mask for ROI (Region of Interest)
mask = np.zeros((H, W))
for i in range(H):
    for j in range(W):
        if i > (H / W) * j and i > -(H / W) * j + H:
            mask[i, j] = 1

# Extract edges in ROI
roi = edges * mask

plt.subplot(1,2,1)
plt.imshow(mask)
plt.title('Mask')
plt.axis('off')

plt.subplot(1,2,2)
plt.imshow(roi)
plt.title('Edges in ROI')
plt.axis('off')
plt.show()

尽管我们利用了mask把边缘限制在了ROI区域，但是所找出的边缘仍旧是一团像素点，而不是连在一起的直线。一般的，我们习惯于用参数方程 $y=ax+b$ 来表示直线。它的斜率是 $a$ 而截距是 $b$ 。在这里，我们利用Hough变换来寻找这条直线。

通常，直线 $y=ax+b$ 可以表示成参数空间 $(a,b)$ 的一个点。但它不能表示竖直方向的直线，所以我们采用了类似于极坐标的参数表达形式 $θ∈[−π,π]$ 和 $ρ∈ℝ$ 表示直线:

利用极坐标的表示方法，我们就将 $xy$ 空间转到了 $θρ$ 空间(或者说Hough空间)。带入不同边缘像素点的坐标 $[x,y]$ ，我们就可以得到θρθρ-空间的一系列曲线，曲线的交点就是我们所求的参数 $θρ$ 。这样，我们就求得了参数坐标下直线的方程，也可以进一步地将其转成 $xy$ 空间的直线方程。

在求取曲线的交点的过程中，我们可以避免直接对交点进行求取，而是某个点被多少条曲线经过。例如，曲线1经过了点 $p_1,p_2,p_3$ ,曲线2经过了点 $p_3,p_4,p_5$ ,那么点 $p_3$ 的计数值就是2，而点 $p_1,p_2,p_4,p_5$ 的计数值就是1，这样通过比较计数值，我们可以确定点 $p_3$ 就是曲线的交点。

from hough import hough_transform

# 对ROI区域的边缘像素点作Hough变换
acc, rhos, thetas = hough_transform(roi)

# 右侧的车道
xs_right = []
ys_right = []

# 左侧的车道
xs_left = []
ys_left = []

#我们认为两条道路在前20条检测到的直线里面
for i in range(20):
    idx = np.argmax(acc)
    r_idx = idx // acc.shape[1]
    t_idx = idx % acc.shape[1]
    acc[r_idx, t_idx] = 0 # 将最大值置0，防止一直检测到同一条直线

    rho = rhos[r_idx]
    theta = thetas[t_idx]
    
    # 将Hough空间（ρ-θ空间）的一个点变换到x-y空间的一条线.
    a = - (np.cos(theta)/np.sin(theta)) # 斜率
    b = (rho/np.sin(theta)) # 截距

    # 检查是否检测到了左右两条车道
    if xs_right and xs_left:
        break
    
    if a < 0: # 左侧车道
        if xs_left:
            continue
        xs = xs_left
        ys = ys_left
    else: # 右侧车道
        if xs_right:
            continue
        xs = xs_right
        ys = ys_right

    for x in range(img.shape[1]):
        y = a * x + b
        if y > img.shape[0] * 0.6 and y < img.shape[0]:
            xs.append(x)
            ys.append(int(round(y)))

plt.imshow(img)
plt.plot(xs_left, ys_left, linewidth=5.0)
plt.plot(xs_right, ys_right, linewidth=5.0)
plt.axis('off')
plt.show()

由此，我们完成了本期对于图像梯度、图像边缘、几何特征、检测与提取的原理与实现。

你可能还会感兴趣

科研基本功——高效文献检索与文献阅读保姆级教程

计算视觉——图像、质量、评价

颜色视觉——杂谈篇

计算视觉——基于暗通道先验的图像去雾算法

写在最后

创作不易，卑微求赞，您的点赞能让更多人看到这篇文章~

谢谢！