Learn OpenCV ---- 大津法（Otsu‘s）阈值

阅读原文

什么是图像阈值

    图像阈值化是一种基于像素强度大图像二值化方法。这种方法的输入通常是一个灰度图和一个阈值，输出是一个二值图像。
    如果图像中某个像素的强度大于阈值，则该像素被标记为白色（前景）；如果图像中的某个像素的强度小于或者等于阈值，则输出像素位置被标记为黑色（背景）。
    图像阈值通常作为预处理步骤，用于使得更为关注的信息突出化。
    简单的图像阈值可以通过手动来设置阈值，进而得到一个好的图像阈值表示。但是这对于大量的图像来说，这是及其乏味和无效的。
    为了克服这一缺点，大津阈值算法被提出，它是一种根据图像来自动确定阈值的算法。在深入了解大津算法之前，我们先了解图像阈值设置与图像分割之间的关系。

图像阈值与图像分割

    图像分割是指将图像分割成不同区域或像素组的一类算法。因此，在某种意义上，图像阈值是图像分割的一种简单形式，即将图像分成两组，一组为前景，一组为背景。下图给出来图像分割的一个小家族。

    图像阈值可以分为局部阈值算法和全局阈值算法。对于全局阈值算法而言，是使用某个固定的阈值来对整个图像进行处理；而局部阈值算法可以利用图像的局部区域信息来选择不同的阈值。值得注意的是，大津算法属于全局阈值算法。

大津（Otsu）算法

    全自动的全局阈值算法通常由下述几步组成：

1. 预处理输入图像。
2. 获取图像直方图（像素分布）。
3. 计算阈值$T$。
4. 将图像中像素大于$T$的区域替换成白色，其余部分替换成黑色。

    可以看到，最主要的步骤在第三步，因此，不同的算法在第三步中有很大的不同。大津法在处理图像直方图后，通过最小化每个类上的方差来分割目标。通常，这种技术为双峰图像产生适当的结果。该类图像的直方图中有两个表达清晰的峰值，代表不同强度值的范围。其核心思想是将图像直方图划分为两个类，并定义一个阈值，使这些类的加权方差${\sigma }_w^2(t)$最小。该加权方差的计算公式为：${\sigma }_w^2(t)=w_1(t){\sigma }_1^2(t)+w_2(t){\sigma }_2^2(t)$，其中，两个权值是两个类别的概率除以阈值$t$，该阈值在[0, 255]范围内取值，计算公式如下：


其中，$P(i)$的计算公式如下：‘

其中，$n$为像素的强度。接着，我们需要求公式中的两个方差，这两个方差的计算公式如下：


可以看到，如果我们想要求得两个方差，我们需要先获得$\mu$的值：


自此，我们就可以得到了所需的加权方差。其实，还有另一种计算方法，这种方法也是更为常用，如下：

具体代码如下：

import cv2
import numpy as np


def otsu_implementation(img_title="boat.jpg", is_normalized=False, is_reduce_noise=False):
    # 使用灰度模式读取图像
    image = cv2.imread(img_title, 0)

    # 是否使用高斯模糊
    # 这里不使用
    if is_reduce_noise:
        image = cv2.GaussianBlur(image, (5, 5), 0)

    # 设置直方图bin的数量
    bins_num = 256

    # 获取图像的直方图
    hist, bin_edges = np.histogram(image, bins=bins_num)

    # 是否对图像进行正则化
    # 这里不进行
    if is_normalized:
        hist = np.divide(hist.ravel(), hist.max())

    # 计算bin的中心
    bin_mids = (bin_edges[:-1] + bin_edges[1:]) / 2.

    # 迭代累加得到权值
    weight1 = np.cumsum(hist)
    weight2 = np.cumsum(hist[::-1])[::-1]

    # 计算得到均值
    mean1 = np.cumsum(hist * bin_mids) / weight1
    mean2 = (np.cumsum((hist * bin_mids)[::-1]) / weight2[::-1])[::-1]

    # 计算最终方差
    inter_class_variance = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:]) ** 2

    # 取最大值对应的下标
    index_of_max_val = np.argmax(inter_class_variance)
    # 从bin中取得阈值
    threshold = bin_mids[:-1][index_of_max_val]
    print("Otsu's algorithm implementation thresholding result: ", threshold)
    return threshold


def main():
    otsu_implementation()


if __name__ == "__main__":
    main()

当然，Opencv也提供了大津算法来进行图像阈值的计算，因此，我们在此也给出使用了Opencv接口的代码：

import cv2
from matplotlib.ticker import FuncFormatter
from matplotlib import pyplot as plt
from otsu_implementation import otsu_implementation


def call_otsu_threshold(img_title="boat.jpg", is_reduce_noise=False):
    # 读取二值图像
    image = cv2.imread(img_title, 0)

    if is_reduce_noise:
        image = cv2.GaussianBlur(image, (5, 5), 0)

    # 现实图像的直方图
    plt.hist(image.ravel(), 256)
    plt.xlabel('Colour intensity')
    plt.ylabel('Number of pixels')
    plt.savefig("image_hist.png")
    plt.show()
    plt.close()

    # 图像二值化，其中THRESH_OTSU就为大津算法
    otsu_threshold, image_result = cv2.threshold(
        image, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU,
    )
    print("Obtained threshold: ", otsu_threshold)

    # 显示两类大概率的直方图
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.hist(image_result.ravel(), 256)
    ax.set_xlabel('Colour intensity')
    ax.set_ylabel('Number of pixels')
    ax.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: ('%1.1fM') % (x*1e-6)))
    plt.savefig("image_hist_result.png")
    plt.show()
    plt.close()

    # 将二值化后的图像显示
    cv2.imshow("Otsu's thresholding result", image_result)
    cv2.waitKey(0)
    cv2.destroyAllWindows()


def main():
    call_otsu_threshold()
    otsu_implementation()


if __name__ == "__main__":
    main()