OpenCV算法学习笔记之边缘检测(二)
in OpenCV with 0 comment

OpenCV算法学习笔记之边缘检测(二)

in OpenCV with 0 comment

Canny边缘检测

原理

基于卷积运算的边缘检测算法(即上一篇OpenCV笔记中提到的),有以下两个缺点:

(1)没有充分利用边缘的梯度方向;

(2)最后输出的边缘二值图,只是简单地利用阈值进行处理,边缘信息损失程度和阈值相关;

Canny算法基于这两点做出了改进,主要为:

(1)基于边缘梯度方向的非极大值抑制;

(2)双阈值的滞后阈值处理;

Canny算法的步骤近似如下:

第一步:图像矩阵$I$分别与水平方向上的卷积核$sobel_x$和垂直方向上的卷积核$sobel_y$卷积得到$dx$和$dy$,然后利用平方和的开方$magnitude = \sqrt{dx^2 + dy^2}$得到边缘强度,这一步和Sobel边缘检测一样,也可以用Prewitt核代替。边界的处理方式是补零

假设有以下矩阵:

第一步

与sobel算子进行卷积后,得到$dx$与$dy$:

第二步

然后计算平方和并开方,得到边缘强度(梯度强度),这里只写整数部分:

第三步

Sobel边缘检测是对magnitude的结果大于255的值截断为255,然后转换为8位图就得到了边缘强度图的灰度显示。而Canny算法的处理方式与此不同。

第二步:利用第一步得到的$dx$与$dy$,计算出梯度方向$angle=arctan2(dy,dx)$,即对每一个位置$(r,c)$,$angle(r,c) = arctan2(dy(r,c), dx(r,c))$代表该位置的梯度方向,一般用角度表示,即$angle(r,c) \in [0, 180] \cup [-180, 0] $。可以利用Python的函数包math中的atan2进行求取,也可以利用Numpy的函数arctan2,代码如下:

y = np.array([[1, 1], [-1, -1]])
x = np.array([[1, -1], [-1, 1]])
angle = np.arctan2(y, x)/np.pi*180

结果为:

array([[45., 135.],
       [-135., -45.]])

用第一步得到的$dx$和$dy$,计算得到每一个位置的梯度角:

第四步

第三步:对每一个位置进行非极大值抑制处理,非极大值抑制处理操作操作返回的仍然是一个矩阵,假设为nonMaxSup。在上面对示例中,边界的处理是补零,导致所得到对magnitude产生额外的边缘响应;如果采用对是以边界为对称的边界扩充方式,那么卷积结果的边界全是0。在非极大值抑制这一步中,对边界不进行任何处理。

所谓对非极大值抑制,就是如果$magnitude(r,c)$在沿着梯度方向$angle(r,c)$上的邻域内是最大的则将$nonMaxSup(r,c)$设置为最大值,否则将其设置为0。

首先将nonMaxSup初始化为以下矩阵:

初始化

接下来以填充$nonMaxSup(1,1)$为例:首先在$magnitude(1,1)$上放置倒置对坐标轴,然后对应到$angle$,发现$angle(1,1)=133$,再按照该梯度方向画出直线,最后以$magnitude(1,1)$为中心对$3*3$邻域内,大体定位出梯度方向上对领域,即右上方和左下方。如下图所示:

梯度邻域

由于$912>292$且$912>276$,故令$nonManSup(1,1)=magnitude(1,1)$。

一般可以将梯度离散为以下几种情况:

对应的领域如下所示:

对应情况

第四步:双阈值对滞后阈值处理。经非极大值抑制后得到的边缘强度图一般需要阈值化处理,常用对是全局阈值分割和局部自适应阈值分割。还有另外一种方法:滞后阈值处理,它使用高阈值和低阈值,按照下列规则进行阈值化处理:

  1. 边缘强度大于高阈值对那些点作为确定边缘点;
  2. 边缘强度比低阈值小的那些点立即被剔除;
  3. 处于高低阈值之间对那些点,只有这些点能按某一路径与确认边缘点相连时,才可以作为边缘点被接受。组成这一路径的所有点都要比低阈值大。实际中,可以先选定边缘强度大于高阈值的所有确定边缘点,然后在边缘强度大于低阈值的情况下尽可能的延长边缘;

Python实现

最好利用dx与dy的平方和开方对形式来衡量边缘强度。

def non_maximum_suppression(dx, dy):
    """
    非极大值抑制
    :param dx: 与水平卷积核卷积运算后的结果
    :param dy: 与垂直卷积核卷积运算后对结果
    :return 非极大值抑制后的矩阵
    """
    # 边缘强度
    edge_mag = np.sqrt(np.power(dx, 2.0) + np.power(dy, 2.0))
    # 高、宽
    rows, cols = dx.shape
    # 梯度方向
    grad_direction = np.zeros(dx.shape)
    # 边缘强度的非极大值抑制
    edge_mag_non_max_sup = np.zeros(dx.shape)
    for r in range(1, rows-1):
        for c in range(1, cols-1):
            # angle的范围[0, 180],[-180, 0]
            angle = math.atan2(dy[r][c], dx[r][c]) / math.pi*180
            grad_direction[r][c] = angle
            # 左右方向
            if abs(angle) < 22.5 or abs(angle) > 157.5:
                if edge_mag[r][c] > edge_mag[r][c-1] and edge_mag[r][c] > edge_mag[r][c+1]:
                    edge_mag_non_max_sup[r][c] = edge_mag[r][c]
            # 左上右下方向
            if (abs(angle) >= 22.5 or abs(angle) < 67.5) or (-angle > 112.5 and -angle <= 157.5):
                if edge_mag[r][c] > edge_mag[r-1][c-1] and edge_mag[r][c] > edge_mag[r+1][c+1]:
                    edge_mag_non_max_sup[r][c] = edge_mag[r][c]
            # 上下方向
            if abs(angle) >= 67.5 and abs(angle) <= 112.5:
                if edge_mag[r][c] > edge_mag[r-1][c] and edge_mag[r][c] > edge_mag[r+1][c]:
                    edge_mag_non_max_sup[r][c] = edge_mag[r][c]
            # 右上左下方向
            if (abs(angle) >= 112.5 or abs(angle) <= 157.5) or (-angle >= 22.5 and -angle < 67.5):
                if edge_mag[r][c] > edge_mag[r-1][c+1] and edge_mag[r][c] > edge_mag[r+1][c-1]:
                    edge_mag_non_max_sup[r][c] = edge_mag[r][c]
    return edge_mag_non_max_sup

接下来是滞后阈值处理:

def check_in_range(r, c, rows, cols):
    """
    判断一个点的坐标是否在图像范围内
    """
    if (r >= 0 and r < rows) and (c >= 0 and c < cols):
        return True
    else:
        return False


def trace(edge_mag_non_max_sup, edge, lower_thresh, r, c, rows, cols):
    """
    滞后阈值处理对第四步:延长边缘
    """
    # 大于高阈值的点为确定边缘点
    if edge[r][c] == 0:
        edge[r][c] = 255
        for i in range(-1, 2):
            for j in range(-1, 2):
                if check_in_range(r+i, c+j, rows, cols) and edge_mag_non_max_sup[r+i][c+j] >= lower_thresh:
                    tarce(edge_mag_non_max_sup, edge, lower_thresh, r+i, c+j, rows, cols)
                    
                    
def hysteresis_threshold(edge_non_max_sup, lower_thresh, upper_thresh):
    """
    滞后阈值处理
    :param edge_non_max_sup: 待处理边缘强度图
    :param lower_thresh: 低阈值
    :param upper_thresh: 高阈值
    :return 处理后对边缘强度图
    """
    # 高宽
    rows, cols = edge_non_max_sup.shape
    edge = np.zeros(edge_non_max_sup.shape, np.uint8)
    for r in range(1, rows-1):
        for c in range(1, cols-1):
            # 大于高阈值的点为确定边缘点,而且以该点为起始点延长边缘
            if edge_non_max_sup[r][c] >= upper_thresh:
                trace(edge_non_max_sup, edge, lower_thresh, r, c, rows, cols)
            # 小于低阈值的被剔除
            if edge_non_max_sup[r][c] < lower_thresh:
                edge[r][c] = 0
    return edge

OpenCV提供函数void Canny(InputArray image, OutputArray edges, double threshold1, double threshold2, int apertureSize=3, bool L2gradient=false),其中threshold1代表低阈值,threshold2代表高阈值,apertureSize代表Sobel核的窗口尺寸,L2gradient代表计算边缘强度时使用的方式,值等于true时用的是平方和形式,否则使用的是绝对值和的形式。

Laplacian边缘检测

原理

二维的拉普拉斯变换由以下公式定义:

$$ \begin{align}\nabla^2f(x,y) &= \frac{\partial^2f(x,y)}{\partial x^2} + \frac{\partial^2f(x,y)}{\partial y^2}\\ &\approx \frac{\partial (f(x+1,y) - f(x,y))}{\partial x} + \frac{\partial (f(x,y+1) - f(x,y))}{\partial y} \\ &\approx f(x+1,y) + f(x-1,y)+f(x,y-1)+f(x,y+1)-4f(x,y) \end{align} $$

推广到二维矩阵,即矩阵进行拉普拉斯变换,也就是与以下卷积核进行卷积(这两种都可以):

$$ l_0 = \left( \begin{matrix} 0&-1&0 \\ -1&4&-1 \\ 0 &-1&0 \end{matrix}\right) , l_{0^-} = \left( \begin{matrix} 0&1&0 \\ 1&-4&1 \\ 0 &1&0 \end{matrix}\right) $$

锚点在中心位置。图像与拉普拉斯核进行卷积运算的本质是计算任意位置的值与其在水平方向和垂直方向上四个相邻点平均值之间的差值(只是相差一个4的倍数)。

拉普拉斯边缘检测没有进行平滑处理,所以对噪点比较敏感;优点是只有一个卷积核,运算量比较低。

拉普拉斯核还有以下几种形式:

$$ l_1 = \left( \begin{matrix} -1&-1&-1 \\ -1&8&-1 \\ -1&-1&-1 \end{matrix}\right) , l_2 = \left( \begin{matrix} 2&-1&2 \\ -1&-4&-1 \\ 2&-1&2 \end{matrix}\right) , l_3 = \left( \begin{matrix} 0&2&0 \\ 2&8&2 \\ 0&2&0 \end{matrix}\right), l_4 = \left( \begin{matrix} 2&0&2 \\ 0&-8&0 \\ 2&0&2 \end{matrix}\right) $$

拉普拉斯核内所有元素的和必须为0,这样就使得在恒等灰度值区域不会产生错误的边缘。上述的拉普拉斯核均是不可分离的。

Python实现

def laplacian(image, _boundary='fill', _fillvalue=0):
    # 第一步:图像与拉普拉斯核卷积
    laplacian_kernel = np.array([[0, -1, 0], [-1, 4, -1], [0, -1, 0]], np.float32)
    # 卷积
    i_conv_lap = signal.convolve2d(image, laplacian_kernel, mode='same', boundary=_boundary, fillvalue=_fillvalue)
    # 第二步:进行阈值化处理
    i_conv_lap[i_conv_lap > 0] = 255
    i_conv_lap[i_conv_lap <= 0] = 0
    i_conv_lap.astype(np.uint8)
    return i_conv_lap

在阈值化处理时,也可以用以下公式进行处理以得到水墨画效果

$$ \begin{equation} abstraction(r,c) = \left\{ \begin{array}{lcl} 255,\quad i\_conv\_lap(r,c) > 0 \\ 0, \quad \quad i\_conv\_lap(r,c) \leq 0 \end{array} \right. \end{equation} $$

OpenCV提供函数void Laplacian(InputArray src, OutputArray dst, int ddepth, int ksize=-1, double scale=1, double delta=0, int borderType=BORDER_DEFAULT)实现拉普拉斯变换,参数解释如下所示:

参数解释
src输入矩阵
dst输出矩阵
ddepth输出矩阵的数据类型(位深)
ksize拉普拉斯核的类型
scale比例系数
delta平移系数
borderType边界扩充类型

参数ksize等于1时采用的是$l_{0^-}$形式的拉普拉斯核,等于3时采用的是$l_4$形式的核。

高斯拉普拉斯(LoG)边缘检测

由于拉普拉斯变换对于噪点比较敏感,在进行变换前先进行高斯平滑再进行卷积,这就是LoG边缘检测。虽然平滑和变换需要两次卷积运算,但是这里我们可以利用一次卷积运算即可完成。

原理

二维高斯函数为

$$ gauss(x,y,\sigma)=\frac{1}{2\pi \sigma^2}\exp(-\frac{(x^2+y^2)}{2\sigma^2}) $$

其拉普拉斯变化如下(省略化简过程):

$$ \begin{align}\nabla^2(gauss(x,y,\sigma)) &= \frac{\partial^2(gauss(x,y,\sigma))}{\partial x^2} + \frac{\partial^2(gauss(x,y,\sigma))}{\partial y^2}\\ &= \frac{1}{2\pi \sigma^4}(\frac{x^2+y^2}{\sigma^2} - 2)\exp(-\frac{x^2+y^2}{2\sigma^2}) \end{align} $$

其步骤如下:

第一步:构建$H×W$、标准差为$\sigma$的LoG卷积核:

$$ LoG_{H×W}=[\nabla^2gauss(w-\frac{W-1}{2}, h-\frac{H-1}{2}, \sigma)]_{0 \leq h <H,0\leq w<W} $$

其中$H,W$均为奇数,且一般$H=W$,锚点位置在$(\frac{H-1}{2}, \frac{W-1}{2})$。

第二步:图像与$LoG_{H×W}$进行卷积,结果记为i_cov_log

第三步:对i_cov_log进行二值化处理,阈值为0;

由于

$$ \nabla^2(gauss(x,y,\sigma))=\frac{1}{\sigma^2}[(\frac{x^2}{\sigma^2} - 1)gauss(x,\sigma)]gauss(y,\sigma)+\frac{1}{\sigma^2}[(\frac{y^2}{\sigma^2} - 1)gauss(y,\sigma)]gauss(x,\sigma) $$

故可以将高斯拉普拉斯核分离为两个卷积核,运算时,图像先与水平方向核卷积,再与垂直卷积核卷积得到结果1;然后图像先与垂直卷积核卷积,再与水平卷积核卷积得到结果2,结果1与结果2相加再进行阈值化处理得到最终结果。

Python实现

这里给出非分离的边缘检测,可以把$\nabla^2(gauss(x,y,\sigma))$前面的系数$\frac{1}{2\pi \sigma^4}$去掉。

def create_LoG_kernel(sigma, size):
    h, w = size
    r, c = np.mgrid[0:h:1, 0:w:1]
    r -= (h-1) / 2
    c -= (w-1) / 2
    # 方差
    sigma2 = pow(sigma, 2.0)
    # 高斯拉普拉斯核
    norm2 = np.power(r, 2.0) + np.power(c, 2.0)
    log_kernel = (norm2/sigma2 - 2) * np.exp(-norm2/(2*sigma2))
    return log_kernel


def LoG(image, sigma, size, _boundary='symm'):
    """
    :param _boundary: symm代表反射扩充
    """
    log_kernel = create_LoG_kernel(sigma, size)
    # 卷积
    img_conv_log = signal.convolve2d(image, log_kernel, 'same', boundary=_boundary)
    return img_conv_log

对于高斯拉普拉斯核的尺寸,一般取$(6*\sigma+1)×(6*\sigma+1)$,即大于$6\sigma$的最小奇数,这样得到的边缘效果会比较好。随着尺度(标准差)的增大,得到的边缘尺度也越来越大,越来越失去边缘的细节。

高斯差分(DoG)边缘检测

可以用差分代替卷积来近似高斯拉普拉斯变换以减小计算量。

原理

由于

$$ \sigma \nabla^2(gauss(x,y,\sigma)) = \frac{\partial gauss(x,y,\sigma)}{\partial \sigma} $$

以及一阶导数的定义:

$$ \begin{align}\frac{\partial gauss(x,y,\sigma)}{\partial \sigma} &= \lim_{k \to 1}\frac{gauss(x,y,k*\sigma) - gauss(x,y,\sigma)}{k*\sigma - \sigma} \\ &\approx \frac{gauss(x,y,k*\sigma) - gauss(x,y,\sigma)}{k*\sigma - \sigma} \end{align} $$

可以得到高斯拉普拉斯的近似值:

$$ \nabla^2 gauss(x,y,\sigma) = \frac{gauss(x,y,k*\sigma) - gauss(x,y,\sigma)}{\sigma^2 (k-1)} $$

该近似值常称为高斯差分。

步骤如下:

第一步:构建$H×W$、标准差为$\sigma$的DoG卷积核:

$$ DoG_{H×W}=[DoG(w-\frac{W-1}{2}, h-\frac{H-1}{2}, \sigma)]_{0 \leq h <H,0\leq w<W} $$

其中$H,W$均为奇数,且一般$H=W$,锚点位置在$(\frac{H-1}{2}, \frac{W-1}{2})$。

第二步:图像与$DoG_{H×W}$进行卷积,结果记为i_cov_dog

第三步:对i_cov_dog进行二值化处理,阈值为0;

参考

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

Responses