【youcans 的 OpenCV 例程200篇】164.使用 Laplace 边缘信息改进全局阈值处理

欢迎关注『youcans 的 OpenCV 例程 200 篇』系列，持续更新中
OpenCV 例程200篇总目录-202205更新

【youcans 的 OpenCV 例程200篇】164.使用 Laplace 边缘信息改进全局阈值处理

3.4 全局阈值处理改进方法

在实际的图像处理中，噪声严重影响阈值处理的结果，严重的噪声会把简单的阈值处理问题变为不能解决的问题。

例程 11.21：使用 Laplace 边缘信息改进全局阈值处理

对于酵母细胞图像，希望通过全局阈值处理等等图像中与亮点对应的区域。如果直接使用 OTSU 方法可以分割细胞区域，但不能检测亮点。

使用 Laplace 算子计算梯度，可以得到亮点的边缘像素，忽略背景区域像素对直方图的贡献，可以改善直方图的分布，从而便于通过阈值处理进行分割。

具体步骤如下：

（1）计算图像 f ( x , y ) f(x,y)f(x,y) 的 Laplace 算子，得到梯度图像；
（2）以灰度值的 99.5% 分位为阈值，对梯度图像进行二值处理，作为遮罩模板，以排除无效背景像素的影响；
（3）基于遮罩模板计算图像 f ( x , y ) f(x,y)f(x,y) 的直方图分布，即只对 g T ( x , y ) = 1 g_T(x,y)=1gT(x,y)=1 的像素进行统计计算；
（4）基于遮罩模板的直方图分布，采用 OTSU 算法计算最佳分割阈值；
（5）使用 OTSU 算法得到的最佳分割阈值，对图像 f ( x , y ) f(x,y)f(x,y) 进行全局阈值处理。

注意本例中用 OTSU 算法求非零像素的最佳分割阈值，不能通过调用 cv2.threshold() 获得。

    # 11.21 使用 Laplace 边缘信息改进全局阈值处理
    img = cv2.imread("../images/Fig1043a.tif", flags=0)

    # # 全局阈值处理，作为参照比较
    histCV1 = cv2.calcHist([img], [0], None, [256], [0, 256])  # 灰度直方图
    ret1, imgOtsu = cv2.threshold(img, 127, 255, cv2.THRESH_OTSU)  # 阈值分割, thresh=T

    # (1) 计算 Laplacian 梯度算子
    laplace = cv2.Laplacian(img, cv2.CV_32F, ksize=3)  # Laplace 卷积算子
    grad = cv2.convertScaleAbs(laplace)
    gradMax = np.int(np.max(grad))

    # (2) 以灰度值的 99.5% 分位为阈值, 对边缘图像进行二值处理, 作为遮罩模板
    per995 = np.percentile(grad, q=99.5)  # 99.5 分位的灰度值, [0, per995] 占比99.5%
    _, gradPer995 = cv2.threshold(np.uint8(grad), per995, 1, cv2.THRESH_BINARY)  # 对边缘图像二值处理

    # (3) 计算基于遮罩模板的直方图分布，以排除无效背景像素的影响
    fp = np.uint8(img * gradPer995)
    histCV2 = cv2.calcHist([fp], [0], None, [256], [0, 256])
    histCV2[0] = 0  # fp 非零像素直方图

    # (4) OTSU 算法计算 fp 非零像素的最佳分割阈值
    # nonzeroPixels = np.count_nonzero(gradPer995)  # 非零像素总数
    nonzeroPixels = sum(histCV2[1:])  # 非零像素总数
    totalGray = np.dot(histCV2[:,0], range(256))  # 内积, 总和灰度值
    mG = totalGray / nonzeroPixels  # 平均灰度
    icv = np.zeros(256)
    numFt, sumFt = 0, 0
    for t in range(0, 256):  # 遍历灰度值
        numFt += histCV2[t,0]   # F(t) 像素数量
        sumFt += histCV2[t,0] * t  # F(t) 灰度值总和
        pF = numFt / nonzeroPixels  # F(t) 像素数占比
        mF = (sumFt/numFt) if numFt>0 else 0  # F(t) 平均灰度
        numBt = nonzeroPixels-numFt  # B(t) 像素数量
        sumBt = totalGray - sumFt  # B(t) 灰度值总和
        pB = numBt / nonzeroPixels  # B(t) 像素数占比
        mB = (sumBt/numBt) if numBt>0 else 0  # B(t) 平均灰度
        icv[t] = pF * (mF-mG)**2 + pB * (mB-mG)**2  # OTSU 算法: 灰度 t 的类间方差
    maxIcv = max(icv)  # ICV 的最大值
    maxIndex = np.argmax(icv)  # 最大值的索引
    print(per995, nonzeroPixels, maxIcv, maxIndex)

    # 使用 fp 非零像素的最佳分割阈值，对原始图像进行固定阈值处理
    ret, imgBin = cv2.threshold(img, maxIndex, 255, cv2.THRESH_BINARY)  # 以 maxIndex 作为最优阈值

    plt.figure(figsize=(10, 7))
    plt.subplot(231), plt.axis('off'), plt.title("Origin"), plt.imshow(img, 'gray')
    plt.subplot(232,yticks=[]), plt.axis([0,255,0,np.max(histCV1)])
    plt.bar(range(256), histCV1[:,0]), plt.title("Gray Hist")
    plt.subplot(233), plt.title("OTSU binary(T={})".format(round(ret1))), plt.axis('off')
    plt.imshow(imgOtsu, 'gray')
    plt.subplot(234), plt.axis('off'), plt.title("Threshold of Laplacian")
    plt.imshow(gradPer995, cmap='gray')  # 遮罩模板，Laplacian 995 分位
    plt.subplot(235, yticks=[]), plt.title("Hist of boundries")  # 直方图
    plt.bar(range(256), histCV2[:,0])
    plt.subplot(236), plt.title("OTSU by Laplacian(T={})".format(maxIndex)), plt.axis('off')
    plt.imshow(imgBin, 'gray')
    plt.show()1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

在这里插入图片描述

（本节完）