canny边缘检测算子

理论原理

Canny 算子是常用的边缘检测算法，其执行步骤如下:

应用高斯滤波器平滑图像以去除噪声
计算图像的强度梯度
应用非最大抑制（Non-maximum suppression）消除边缘检测的虚假响应
应用双阈值确定潜在边缘
通过滞后 (hysteresis ) 方法跟踪边缘：通过抑制所有其他弱边缘和未连接到强边缘的边缘，完成边缘检测

高斯滤波

通过高斯滤波去除图像噪声的影响，常用 $5 \times 5$ 大小的高斯核，设置 $\sigma=1.4$
$$
\mathbf{B}=\frac{1}{159}\left[\begin{array}{ccccc}
2 & 4 & 5 & 4 & 2 \\
4 & 9 & 12 & 9 & 4 \\
5 & 12 & 15 & 12 & 5 \\
4 & 9 & 12 & 9 & 4 \\
2 & 4 & 5 & 4 & 2
\end{array}\right] * \mathbf{A}
$$

高斯滤波核尺寸越大，越能够平滑噪声影响，但与此同时 Canny 算子的边缘检测的性能会降低

图像梯度计算

通过 Sobel 算子计算图像的水平和垂直反向导数，然后计算梯度大小和方向
$$
\begin{aligned}
& G=\sqrt{G_x^2+G_y^2} \
& \theta=\arctan \left(\frac{G_y}{G_x}\right)
\end{aligned}
$$

将梯度方向通过四舍五入方法归入到水平/垂直/对角 $\left(0^{\circ}, 45^{\circ}, 90^{\circ}\right.$ 和 $\left.135^{\circ}\right)$ ，比如 $\left[0^{\circ}, 22.5^{\circ}\right]$ 和 $\left[157.5^{\circ}, 180^{\circ}\right]$ 映射为 $0^{\circ}$

非最大抑制

非最大抑制是一种边缘细化技术。进行梯度计算后的图像边缘仍旧很模糊，边缘拥有多个候选位置，所以需要应用非最大抑制来寻找“最大”像素点，即具有最大强度值变化的位置，移除其他梯度值，保证边缘具有准确的响应。其原理如下

将当前像素的边缘强度与像素在正梯度方向和负梯度方向上的边缘强度进行比较
如果与相同方向的掩模中的其他像素相比，当前像素的边缘强度是最大的(例如，指向 $\mathrm{y}$ 轴方向的像素将与垂直轴上方和下方的像素相比较) ，则该值将被保留。否则，该值将被抑制(去除梯度值为 0 )

具体实现时，将连续梯度方向分类为一组小的离散方向，然后在上一步的输出(即边缘强度和梯度方向 )上移动 $3 \times 3$ 滤波器。在每个像素处，如果中心像素的大小不大于渐变方向上两个相邻像素的大小，它将抑制中心像素的边缘强度(通过将其值设置为 0 )

如果梯度角度为 $0^{\circ}$ （即边缘在南北方向），如果其梯度大小大于东西方向像素处的大小，则该点将被视为在边缘上
如果梯度角度为 $45^{\circ}$ （即边缘位于西北-东南方向），如果其梯度大小大于东北和西南方向像素处的大小，则该点将被视为位于边缘上
如果梯度角度为 $90^{\circ}$ （即边缘在东西方向），如果其梯度大小大于南北方向像素处的大小，则该点将被视为在边缘上
如果梯度角度为 $135^{\circ}$ （即边缘位于东北-西南方向），如果其梯度大小大于西北和东南方向像素处的大小，则该点将被视为位于边缘上

双边阈值

通过非最大抑制，可以有效确定边缘的最大像素点，剩余的边缘像素提供了图像中真实边缘的更精确表示。但是，由于噪声和颜色变化，一些边缘像素仍然存在。为了去除这些虚假响应，必须滤除具有弱梯度值的边缘像素，并保留具有高梯度值的边缘像素。这是通过选择高阈值和低阈值来实现的。如果边缘像素的渐变值高于高阈值，则将其标记为强边缘像素。如果边缘像素的渐变值小于高阈值且大于低阈值，则将其标记为弱边缘像素。如果边缘像素的值小于低阈值，它将被抑制

两个阈值是通过经验确定的，其定义将取决于给定输入图像的内容。通过其比率（ $upper: lower$ ）设置为 $2: 1$ 和 $3: 1$ 之间

通过滞后方法进行边缘追踪

经过上述步骤处理后，结果图像仅包含了强边缘像素和弱边缘像素。对于弱边缘像素而言，这些像素既可以从真实边缘提取，也可以从噪声/颜色变化中提取

通常，由真实边缘引起的弱边缘像素将连接到强边缘像素，而噪声响应则不连接。为了跟踪边缘连接，通过观察弱边缘像素及其 8 个邻域像素来进行blob分析。只要 blob 中包含一个强边缘像素，就可以将该弱边缘点识别为一个应该保留的点

具体内容

import cv2
import numpy as np
import sys
sys.setrecursionlimit(3000)
def trace(M_supp,edge,i,j,LT):
    if edge[i,j]==0:
        edge[i,j]=255
        for di,dj in[(-1,-1),(-1,0),(-1,1),(0,-1),(0,1),(1,-1),(1,0),(1,1)]:
            if M_supp[i+di,j+dj]>=LT:
                trace(M_supp,edge,i+di,j+dj,LT)
HT=200
LT=50
src=cv2.imread('E:\shabi.jpg?raw=true')
dst=cv2.cvtColor(src,cv2.COLOR_BGRA2GRAY)
blur_dst=cv2.GaussianBlur(dst,(5,5),0.5)
dx=cv2.Sobel(blur_dst,cv2.CV_16S,1,0)
dy=cv2.Sobel(blur_dst,cv2.CV_16S,0,1)
M=abs(dx)+abs(dy)
theta=np.arctan2(dx,dy)
M_supp=np.zeros(M.shape)
for i in range(1,M.shape[0]-1):
    for j in range(1,M.shape[1]-1):
        x=dx[i,j]
        y=dy[i,j]
        angle=theta[i,j]/np.pi
        mag=M[i,j]
        if abs(angle)<=1/8. or abs(angle)>=7/8.:
            if mag>=M[i-1,j] and mag>=M[i+1,j]:
                M_supp[i,j]=mag
        elif abs(angle-1/2.)<=1/8. or abs(angle+1/2.)<=1/8.:
            if mag>=M[i,j-1] and mag>=M[i,j+1]:
                M_supp[i,j]=mag
        elif abs(angle-3/4.)<=1/8. or abs(angle+3/4.)<=1/8.:
            if mag>=M[i+1,j-1] and mag>=M[i-1,j+1]:
                M_supp[i,j]=mag
        else:
            if mag>=M[i+1,j+1] and mag>=M[i-1,j-1]:
                M_supp[i,j]=mag
M_supp=np.pad(M_supp,((1,1)),'constant')
edge=np.zeros(M_supp.shape,dtype=np.uint8)
for i in range(1,M_supp.shape[0]-1):
    for j in range(1,M_supp.shape[1]-1):
        if M_supp[i,j]>=HT:
            trace(M_supp,edge,i,j,LT)
edge=edge[1:-1,1:-1]
cv2.imshow('edge',edge)
cv2.imwrite('edge.jpg',edge)
cv2.waitKey(0)

原始图像

边缘检测

基于阈值的二值化分割

人工选取阈值

def human(gray_img,threhold):
    dst=gray_img.copy()
    dst[dst<threhold]=0
    dst[dst>=threhold]=255
    return dst
dst2=human(src,150)

分割后图像

迭代法

def iteration(gray_img):
    dst=gray_img.copy()
    threhold=128
    threhold1=0
    while(1):
        if(abs(threhold-threhold1)<=1):break
        else: threhold=threhold1
        n0=dst[np.where(dst<threhold)]
        n1=dst[np.where(dst>=threhold)]
        miu0=n0.mean() if len(n0)>0 else 0
        miu1=n1.mean() if len(n1)>0 else 0
        threhold1=(miu0+miu1)/2
    dst[dst<threhold]=0
    dst[dst>=threhold]=255
    return dst

分割结果

OSTU算法

理论原理

大津算法（OTSU算法）是一种常用的图像二值化方法，用于将灰度图像转化为二值图像。该算法由日本学者大津展之于1979年提出，因此得名。

大津算法的核心思想是通过寻找一个阈值，将图像的像素分为两个类别：前景和背景。具体步骤如下：

统计图像的灰度直方图，得到每个灰度级的像素数目。
遍历所有可能的阈值（0到255），计算根据该阈值将图像分为前景和背景的类内方差。
根据类内方差的最小值确定最佳阈值。
在大津算法中，类内方差是衡量前景和背景之间差异的度量。通过选择使类内方差最小的阈值，可以实现最佳的图像分割效果。

大津算法的优点是简单易懂，计算效率高。它适用于灰度图像的二值化处理，特别是对于具有双峰直方图的图像效果更好。然而，该算法对于具有非双峰直方图的图像可能产生较差的分割结果。因此，在应用大津算法之前，需要对图像的直方图进行分析，确保适用性。

大津算法在图像处理中被广泛应用，例如在文档图像处理、目标检测、图像分割等领域。

下面推导类间方差函数：
设阈值为灰度 $\mathrm{k}(k \in[0, L-1], L=256)$ 。这个阈值把图像像素分割成两类，C1类像素小于等于 $\mathrm{k} ， \mathrm{C} 2$ 类像素大于 $\mathrm{k}$ 。设这两类像素各自的均值为 $m_1, m_2$ ，图像全局均值为 $m_G$ 。同时像素被分为 $\mathrm{C} 1$ 和 $\mathrm{C} 2$ 类的概率分别为 $p_1, p_2$ 。则有:
$$
\begin{gathered}
p_1 m_1+p_2 m_2=m_G \\
p_1+p_2=1
\end{gathered}
$$

根据方差的概念，类间方差表达式为:
$$
\sigma^2=p_1\left(m_1-m_G\right)^2+p_2\left(m_2-m_G\right)^2
$$

展开:
$$
\sigma^2=p_1 m_1^2+p_1 m_G^2-2 p_1 m_1 m_G+p_2 m_2^2+p_2 m_G^2-2 p_2 m_2 m_G
$$

合并2，5及3，6项可得:
$$
\sigma^2=p_1 m_1^2+p_2 m_2^2+m_G^2-2 m_G^2=p_1 m_1^2+p_2 m_2^2-m_G^2
$$

我们再把 $m_G=p_1 m_1+p_2 m_2$带回得
$$
\sigma^2=\left(p_1-p_1^2\right) m_1^2+\left(p_2-p_2^2\right) m_2^2-2 p_1 p_2 m_1 m_2
$$

再注意到 $p_1+p_2=1$ ，所以 $p_1-p_1^2=p_1\left(1-p_1\right)=p_1 p_2 ， p_2-p_2^2=p_2\left(1-p_2\right)=p_1 p_2$ ，从而得到:
$$
\sigma^2=p_1 p_2\left(m_1-m_2\right)^2
$$

对于给定的阈值 $k$ ，我们可以统计出灰度级的分布列:

灰度值	0	1	$\ldots$	255
$p_i$	$p_0$	$p_1$	$\ldots$	$p_{255}$

显然根据分布列性质有 $\sum_{i=0}^{L-1} p_i=1$ (请注意这里的 $p_1, p_2$ 是分布列中的，不是上面的定义)那么有:
$$
p_1=\sum_{i=0}^{k-1} p_i, p_2=\sum_{i=k}^{L-1} p_i, m_1=\sum_{i=0}^{k-1} i p_i, m_2=\sum_{i=k}^{L-1} i p_i
$$

将 $\mathrm{k}$ 从 $[0, L-1]$ 遍历，找出使得 $\sigma^2$ 最大的 $\mathrm{k}$ 值，这个 $\mathrm{k}$ 值就是阈值。
对于分割，这个分割就是二值化

具体内容

def ostu(gray_img):
    dst=gray_img.copy()
    h=dst.shape[0]
    w=dst.shape[1]
    threhold_t=0
    max_g=0
    for t in range(255):
        n0=dst[np.where(dst<t)]
        n1=dst[np.where(dst>=t)]
        w0=len(n0)/(h*w)
        w1=len(n1)/(h*w)
        u0=np.mean(n0) if len(n0)>0 else 0.
        u1=np.mean(n1) if len(n1)>0 else 0.
        g=w0*w1*(u0-u1)**2
        if g>max_g:
            max_g=g
            threhold_t=t
    dst[dst<threhold_t]=0
    dst[dst>=threhold_t]=255
    return dst

分割图像