Otsu‘s Thresholding的工作原理

在双峰图像（像素的直方图只包含两个峰值），一个好的阈值将在这两个值的中间。Otsu算法就是从图像的像素直方图中确定了一个最佳的全局阈值 t t t，可以最小化weighted within-class variance加权类内方差。

σ w 2 ( t ) = q 1 ( t ) σ 1 2 ( t ) + q 2 ( t ) σ 2 2 ( t ) sigma_w^2(t) = q_1(t)sigma_1^2(t)+q_2(t)sigma_2^2(t) σw2​(t)=q1​(t)σ12​(t)+q2​(t)σ22​(t)
权重：
q 1 ( t ) = ∑ i = 1 t P ( i ) & q 2 ( t ) = ∑ i = t + 1 I P ( i ) q_1(t) = sum_{i=1}^{t} P(i) quad & quad q_2(t) = sum_{i=t+1}^{I} P(i) q1​(t)=i=1∑t​P(i)&q2​(t)=i=t+1∑I​P(i)
均值：
μ 1 ( t ) = ∑ i = 1 t i P ( i ) q 1 ( t ) & μ 2 ( t ) = ∑ i = t + 1 I i P ( i ) q 2 ( t ) mu_1(t) = sum_{i=1}^{t} frac{iP(i)}{q_1(t)} quad & quad mu_2(t) = sum_{i= t+1}^{I} frac{iP(i)}{q_2(t)} μ1​(t)=i=1∑t​q1​(t)iP(i)​&μ2​(t)=i=t+1∑I​q2​(t)iP(i)​
方差：
σ 1 2 ( t ) = ∑ i = 1 t [ i − μ 1 ( t ) ] 2 P ( i ) q 1 ( t ) & σ 2 2 ( t ) = ∑ i = t + 1 I [ i − μ 2 ( t ) ] 2 P ( i ) q 2 ( t ) sigma_1^2(t) = sum_{i=1}^{t} [i-mu_1(t)]^2 frac{P(i)}{q_1(t)} quad & quad sigma_2^2(t) = sum_{i=t+1}^{I} [i-mu_2(t)]^2 frac{P(i)}{q_2(t)} σ12​(t)=i=1∑t​[i−μ1​(t)]2q1​(t)P(i)​&σ22​(t)=i=t+1∑I​[i−μ2​(t)]2q2​(t)P(i)​

import cv2 as cv
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

img = cv.imread('./images/kakaka.jpg')
img = cv.cvtColor(img, cv.COLOR_BGR2RGB)
img_gray = cv.cvtColor(img, cv.COLOR_RGB2GRAY)
#img_blur = cv.GaussianBlur(img_gray, (5, 5), 0)

hist_cv = cv.calcHist([img_gray],[0],None,[256],[0,256]) # （256，1）
hist_np = np.bincount(img_gray.ravel(), minlength=256)   # （256，）
hist_norm = hist_np/hist_np.sum()
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])
    q1, q2 = Q[i], 1.0 - Q[i]    # Q[255] = 1.0
    if q1 < 1.e-6 or q2 < 1.e-6:
        continue
    b1, b2 = np.hsplit(bins, [i])
    # 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
    # weighted within-class variances
    fn = v1*q1 + v2*q2
    if fn < fn_min:
        fn_min = fn
        thresh = i
print("threshold t : %d"%thresh)

threshold t : 115

mask = np.zeros_like(img_gray)
mask[np.where(img_gray >= 115)] = 1
plt.imshow(img*mask.reshape(245, 250, 1))

ret, otsu = cv.threshold(img_gray, 0, 255, cv.THRESH_BINARY+cv.THRESH_OTSU)
print("threshold t : %d"%ret)

threshold t : 114

img_roi = cv.bitwise_and(img, img ,mask = otsu)