Otsu算法

简介

在图像处理中我们经常会遇到（灰度图）二值化的需求，一般情况下可以指定全局阈值进行二值化，那这种时候我们如何知道这个阈值的好坏呢？答案就是不停的尝试。如果是一副双峰图像（双峰图像是指图像直方图中存在两个峰），直观上应该是在两峰之间的峰谷附近取一个值。Otsu就是处理这个问题的，即根据图像直方图自动算出一个阈值。本文会先介绍和简单推导一下公式，然后再脑洞清奇地用代码探索一下Otsu的二值化和kmean的二类聚类方法的效果是否有差异。

原理与公式

原理

Otsu算法的原理非常简单，算法假定该图像根据双模直方图（前景像素和背景像素）把包含两类像素，于是它要计算能将两类分开的最佳阈值，使得它们的类内方差最小；由于两两平方距离恒定，所以即它们的类间方差最大。

公式推导

在大津算法中，我们穷举搜索能使类内方差最小的阈值\(t\)。类内方差定义为两个类的方差的加权和： \[\sigma_w^2(t) = \omega_1(t)\sigma_1^2(t)+\omega_2(t)\sigma_2^2(t) \tag{1}\] 权重 \(\omega_i\)是被阈值\(t\)分开的两个类的概率，而 \(\sigma^2_ i\)是这两个类的类内方差。有： \[\omega_1(t)=\sum_{i=0}^t p(i) \tag{2}\] \[\omega_2(t)=\sum_{i=t+1}^l p(i) \tag{3}\] 类均值\(\mu_i(t)\)为： \[\mu_1(t)=\left[\sum_{i=0}^t p(i)\,x(i)\right]/\omega_1 \tag{4}\] \[\mu_2(t)=\left[\sum_{i=t+1}^l p(i)\,x(i)\right]/\omega_2 \tag{5}\] 其中\(x(i)\)为第\(i\)个直方图面元中心的值。 而类间方差： \[\sigma_b^2(t) = \sigma^2-\sigma_w^2(t)=\omega_1(t)\omega_2(t)[\mu_1(t)-\mu_2(t)]^2 \tag{6}\] 可知类间方差最大与类内方差最小等价。

代码

我们就采用使得类内方差(1)式最小的方案，Python代码如下

import cv2
import numpy as np
img = cv2.imread('test.jpg', 0)
blur = cv2.GaussianBlur(img, (5, 5), 0)

# find normalized_histogram, and its cumulative distribution function
hist = cv2.calcHist([blur], [0], None, [256], [0, 256])
hist_norm = hist.ravel()/hist.max()
Q = hist_norm.cumsum()
bins = np.arange(256)
fn_min = np.inf
thresh = -1
for i in range(1, 256):
    p1, p2 = np.hsplit(hist_norm, [i])  # probabilities
    q1, q2 = Q[i], Q[255]-Q[i]  # cum sum of classes
    b1, b2 = np.hsplit(bins, [i])  # weights
    # finding means and variances
    m1, m2 = np.sum(p1*b1)/q1, np.sum(p2*b2)/q2
    v1, v2 = np.sum(((b1-m1)**2)*p1)/q1, np.sum(((b2-m2)**2)*p2)/q2
    # calculates the minimization function
    fn = v1*q1 + v2*q2
    if fn < fn_min:
        fn_min = fn
        thresh = i
print("loop:", i)
print("threshold:", thresh)

kmean二值化思考

这是笔者的一个想法，直接在直方图上做kmean二类聚类，看看是否能够达到和Otsu算法一样的效果。我们先按kmean的思路模拟一下算法流程：

1. 在直方图上任取两个整数点\(A_0\),\(B_0\)，\(A_0,B_0\in[0,255]\)；
2. 在直方图上的其余点按距离归属为二类，易知\(A_i\),\(B_i\)的中点\(M_i\)即为分界点，\(M_i=(A_i+B_i)/2\)；
3. 分别计算被\(M_i\)分割的两类的加权均值，并用两类的加权均值替换\(A_{i+1},B_{i+1}\)；
4. 重复第2步和第3步，直到\(M_i\)变化幅度小于\(eps\)。

kmean代码

import cv2
import numpy as np
img = cv2.imread('test.jpg', 0)
blur = cv2.GaussianBlur(img, (5, 5), 0)

# find normalized_histogram, and its cumulative distribution function
hist = cv2.calcHist([blur], [0], None, [256], [0, 256])
hist_norm = hist.ravel()/hist.max()
Q = hist_norm.cumsum()
bins = np.arange(256)
fn_min = np.inf
thresh = -1
A = 64
B = 192
eps = 1
count = 0
for t in range(1, 256):
    i = int(round((A+B)/2))
    p1, p2 = np.hsplit(hist_norm, [i])  # probabilities
    q1, q2 = Q[i], Q[255]-Q[i]  # cum sum of classes
    b1, b2 = np.hsplit(bins, [i])  # weights
    # finding means and variances
    A, B = np.sum(p1*b1)/q1, np.sum(p2*b2)/q2
    if abs(((A+B)/2)-i) < eps:
        count = t
        break
print("loop:", count)
print("threshold:", i)

结果比较

运行上面两段代码如下：

Otsu

loop: 255
threshold: 120

kmean

loop: 2
threshold: 120

从代码上分析也可以看到Otsu法事实上就是遍历直方图的横轴使得类内方差最小即可，迭代次数固定。而kmean方法一般情况下4次左右迭代即可找到与Otsu法一样的结果。事实上第二种算法也是用到了EM算法（可见最大期望算法原理与推导）的相关思想。

参考

大津算法