165.多阈值 OTSU 处理方法

9. 阈值处理与边缘检测索引

一、多阈值处理方法

• 在实际的图像处理中，噪声严重影响阈值处理的结果，严重的噪声会把简单的阈值处理问题变为不能解决的问题。
  OTSU 方法使用最大化类间方差（intra-class variance）作为评价准则，基于对图像直方图的计算，可以给出类间最优分离的最优阈值。
  OTSU 方法可以扩展到任意数量的阈值。假设有 K 个分类 时，可以定义类间方差为：
  

  
  常用地，考虑由 3个灰度区间组成的 3个类，可以由 2个阈值分割，定义类间方差为：
  
  使类间方差 ICV 最大化的灰度值 就是最优阈值。
  阈值处理后的图像由下式给出：
  

• 需要指出的是，双阈值处理涉及两重循环，算法耗时较长，但仍是基于直方图进行运算的，因此用时与图像无关。对于双重循环可以进行优化，就不在本节介绍了。
  如果需要使用更多个阈值进行分割，通常不再用循环遍历的方法，而是采用聚类或启发式方法来获得分割阈值。

二、例程

• 11.22：阈值处理之多阈值 OTSU

    # 11.22 阈值处理之多阈值 OTSU
    def doubleThreshold(img):
        histCV = cv2.calcHist([img], [0], None, [256], [0, 256])  # 灰度直方图
        grayScale = np.arange(0, 256, 1)  # 灰度级 [0,255]
        totalPixels = img.shape[0] * img.shape[1]  # 像素总数
        totalGray = np.dot(histCV[:,0], grayScale)  # 内积, 总和灰度值
        mG = totalGray / totalPixels  # 平均灰度，meanGray
        varG = sum(((i-mG)**2 * histCV[i,0]/totalPixels) for i in range(256))

        T1, T2, varMax = 1, 2, 0.0
        # minGary, maxGray = np.min(img), np.max(img)  # 最小灰度，最大灰度
        for k1 in range(1, 254):  # k1: [1,253], 1<=k1<k2<=254
            n1 = sum(histCV[:k1, 0])  # C1 像素数量
            s1 = sum((i * histCV[i, 0]) for i in range(k1))
            P1 = n1 / totalPixels  # C1 像素数占比
            m1 = (s1 / n1) if n1 > 0 else 0  # C1 平均灰度

            for k2 in range(k1+1, 256):  # k2: [2,254], k2>k1
                # n2 = sum(histCV[k1+1:k2,0])  # C2 像素数量
                # s2 = sum( (i * histCV[i,0]) for i in range(k1+1,k2) )
                # P2 = n2 / totalPixels  # C2 像素数占比
                # m2 = (s2/n2) if n2>0 else 0  # C2 平均灰度
                n3 = sum(histCV[k2+1:,0])  # C3 像素数量
                s3 = sum((i*histCV[i,0]) for i in range(k2+1,256))
                P3 = n3 / totalPixels  # C3 像素数占比
                m3 = (s3/n3) if n3>0 else 0  # C3 平均灰度

                P2 = 1.0 - P1 - P3  # C2 像素数占比
                m2 = (mG - P1*m1 - P3*m3)/P2 if P2>1e-6 else 0  # C2 平均灰度

                var = P1*(m1-mG)**2 + P2*(m2-mG)**2 + P3*(m3-mG)**2
                if var>varMax:
                    T1, T2, varMax = k1, k2, var

        epsT = varMax / varG  # 可分离测度
        print(totalPixels, mG, varG, varMax, epsT, T1, T2)
        return T1, T2, epsT

    img = cv2.imread("../images/Fig1043a.tif", flags=0)
    # img = cv2.imread("../images/Fig1045a.tif", flags=0)
    histCV = cv2.calcHist([img], [0], None, [256], [0, 256])  # 灰度直方图

    T1, T2, epsT = doubleThreshold(img)
    print("T1={}, T2={}, esp={:.4f}".format(T1, T2, epsT))

    binary = img.copy()
    binary[binary<T1] = 0
    binary[binary>T2] = 255

    ret, imgOtsu = cv2.threshold(img, 127, 255, cv2.THRESH_OTSU)  # OTSU 阈值分割
    ret1, binary1 = cv2.threshold(img, T1, 255, cv2.THRESH_TOZERO)  # 小于阈值置 0，大于阈值不变
    ret2, binary2 = cv2.threshold(img, T2, 255, cv2.THRESH_TOZERO)

    plt.figure(figsize=(9, 6))
    plt.subplot(231), plt.axis('off'), plt.title("Origin"), plt.imshow(img, 'gray')
    plt.subplot(232,yticks=[]), plt.axis([0,255,0,np.max(histCV)])
    plt.bar(range(256), histCV[:,0]), plt.title("Gray Hist")
    plt.subplot(233), plt.title("OTSU binary(T={})".format(round(ret))), plt.axis('off')
    plt.imshow(imgOtsu, 'gray')
    plt.subplot(234), plt.axis('off'), plt.title("Threshold(T={})".format(T1))
    plt.imshow(binary1, 'gray')
    plt.subplot(235), plt.axis('off'), plt.title("Threshold(T={})".format(T2))
    plt.imshow(binary2, 'gray')
    plt.subplot(236), plt.axis('off'), plt.title("DoubleT({},{})".format(T1,T2))
    plt.imshow(binary, 'gray')
    plt.show()

运行结果：
Fig1043a.tif：T1=35, T2=101, esp=0.8733

Fig1045a.tif：T1=81, T2=177, esp=0.9540

三、资料

youcans_的博客：
https://blog.csdn.net/youcans/article/details/124400439