最近，しきい値を自動的に決定したいという話をよく相談されるので，まとめてみました．
しきい値設計とは，1次元の教師なし2グループのクラスタリングに対応し
ます．
今回は，画像の2値化を例として，5個の古典的かつ有用な手法をまとめます．
注意: 分布系が未知ならば，しきい値決定には正解はありません
準備
今回は適当にとったこの写真を使います．
この画像の輝度のヒストグラムは以下の様になります．
この輝度の分布から最適な閾値を設定し，バイナリー画像を作るという問題になります．
K-means [1]
教師なし学習で最も有名な手法は，K-meansです (wiki).
K-meansは以下の最適化問題の解となるそれぞれのクラスの中心点 {m1,...,mM}を探すことを目的としています．
$$
{\rm minimize}_{m_1,\ldots, m_M} \quad \sum_{i=1}^N \min_j |x_i - m_j|
$$
ここで{x1,...,xN}は与えられたデータになります．
クラスの中心点{m1,...,mM)}は以下の2つを繰り返すことによって求められます．
1. クラス内の平均値の計算する．
2. データをもっとも近い中心のクラスに割り当て直す．
K-meansは多次元の場合でも適用できますが，閾値設計の問題ではM=2となります．
この分類の結果は次の図のようになります．
K-meansの特徴としては，以下のような特徴があります．
- 更新ごとにN×Mオーダーの距離を計算する必要がある．
- 初期値依存性が存在する．
import cv2
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# データの読み込み
img = cv2.imread('img/main.png',0)
data = img.reshape(-1)
# ラベルの初期化
labels = np.random.randint(0,2,data.shape[0])
# 終了条件
OPTIMIZE_EPSILON = 1
m_0_old = -np.inf
m_1_old = np.inf
for i in range(1000):
# それぞれの平均の計算
m_0 = data[labels==0].mean()
m_1 = data[labels==1].mean()
# ラベルの再計算
labels[np.abs(data-m_0) < np.abs(data-m_1)] = 0
labels[np.abs(data-m_0) >= np.abs(data-m_1)] = 1
# 終了条件
if np.abs(m_0 - m_0_old) + np.abs(m_1 - m_1_old) < OPTIMIZE_EPSILON:
break
m_0_old = m_0
m_1_old = m_1
# 初期値によって，クラスが変化するため上界の小さい方を採用
thresh_kmeans = np.minimum(data[labels==0].max(),data[labels==1].max())
大津法 [2]
大津法は各クラスの重み付き平均値の差を最大化する方法です(wiki)．
具体的には以下の最適化問題を解くしきい値Tを見つけることに対応します．
$$
{\rm maximize}_T\quad \sigma_b^2 (T):= w_1 (m_1 - m)^2 +w_2 (m_2 - m)^2
$$
ここで，
$$
w_1 = \frac{ | {x_1,\ldots, x_N}\leq T|}{N},\quad w_2 = \frac{ | {x_1,\ldots, x_N} > T|}{N}
$$
$$
m_1 = \frac{\sum_{x_i\leq T} x_i}{\sum x_i},\quad m_2 = \frac{\sum_{x_i>T} x_i}{\sum x_i}
$$
となります．
上記のσbを最大化する問題は，各クラスの分散の重み付き和
$$
\sigma_a^2 (T) := w_1 \frac{\sum_{x_i\leq T} (x_i - m_1)^2}{N} +w_2 \frac{\sum_{x_i>T} (x_i - m_2)^2}{N}
$$
を最小化する問題と一致します．
そのため，σa最小化の問題やσa/σbを最小化する問題を解くこともあります．
この分類の結果は次の図のようになります．
以下にpythonでのコードを記述します．
# Define Otsu scoring function
def OtsuScore(data,thresh):
w_0 = np.sum(data<=thresh)/data.shape[0]
w_1 = np.sum(data>thresh)/data.shape[0]
# check ideal case
if (w_0 ==0) | (w_1 == 0):
return 0
mean_all = data.mean()
mean_0 = data[data<=thresh].mean()
mean_1 = data[data>thresh].mean()
sigma2_b = w_0 *((mean_0 - mean_all)**2) + w_1 *((mean_1 - mean_all)**2)
return sigma2_b
# Callculation of Otsu score and analyze the optimal
scores_otsu = np.zeros(256)
for i in range(scores_otsu.shape[0]):
scores_otsu[i] = OtsuScore(data,i)
thresh_otsu = np.argmax(scores_otsu)
Sezanらによる方法 [3]
画像の輝度の分布などで，ピークが2つあるということを利用したしきい値設計手法です．
2つのピークの内側のすその値を求めることにより，しきい値を設計します．
まず，平滑化したヒストグラムを極地計算により，一番左側（値が小さい）のピークと一番右側（値が大きい）のピークを検出し，それぞれ[m0,m1]とします．
同じく極地の計算によって，各ピークの左右の裾を計算し，左の裾を[s0,s1], 右の裾を[e0,e1]とします．
このとき，しきい値は左側のピークの右側の裾e0と右側のピークの左側の裾s1の間にあると考えられます．
そのため，0以上1以下の設計パラメーターγを用いて，しきい値を設計します．
$$
{\rm Sezan~threshold} := (1 - \gamma )e_0 + \gamma s_1
$$
この分類の結果は次の図のようになります．
Sezanらによる2値の手法の特徴として以下があります．
- 調整するハイパーパラメーターが多いが状況に合わせて使い分けることができる
- 多次元化が難しい
pythonコードは以下のようになります．
# 平滑化パラメータ
sigma = 5.0
# 重要視する割合
# gamma = 1 : 黒の領域を広めに取る
# gamma = 0 : 白の領域を広めに取る
gamma = 0.5
# 平滑化のためのガウスカーネル
def getGaus(G_size,G_sigma):
G_kernel = np.zeros(G_size)
G_i0 = G_size//2
for i in range(G_size):
G_kernel[i] = np.exp(-(i-G_i0)**2 / (2*G_sigma**2))
# 和が1になるように調整
G_kernel = G_kernel/G_kernel.sum()
return G_kernel
# ガウスカーネルの作成
kernel = getGaus(55,sigma)
# ヒストグラム化
num_hist, range_hist = np.histogram(data, bins= 256)
mean_hist = (range_hist[1:] + range_hist[:-1]) / 2
# 平滑化
hist_bar = np.convolve(num_hist,kernel,'same')
# 差分の計算
d_hist = hist_bar[:-1] - hist_bar[1:]
# 端の処理
d_hist = np.r_[[0],d_hist,[0]]
# ピーク検出
m = np.where((d_hist[1:] >=0) & (d_hist[:-1] <=0))[0]
# 局地検出
es =np.where((d_hist[1:] <=0) & (d_hist[:-1] >=0))[0]
# 最大ピークと最小ピーク
m0 = m.min()
m1 = m.max()
# ピーク前後の局地
s0 = es[es<m0].max()
e0 = es[es>m0].min()
s1 = es[es<m1].max()
e1 = es[es>m1].min()
# 閾値決定
thresh_Sezan = (1 - gamma) * mean_hist[e0] + gamma * mean_hist[s1]
カルバック・ライブラー情報量最小化 [4]
情報理論の研究者におなじみのカルバック・ライブラー(KL)情報量を使う手法です．
KL情報量以外の他の情報量基準で研究されている論文も多数あります．
まず，輝度のヒストグラムを正規化した確率分布をpと定義します．
次に，しきい値Tによって定まる2値化に対応する確率関分布q(T)を以下のように定義します．
$$
q(i\leq T) = 0,\quad q(i>T) = 1/M
$$
ここで，Mはしきい値より大きい要素の数です．
q(T)とpによるKL情報量D(q(T)||p)を最小化する値をしきい値とします．
つまり，以下の最適化問題を解きます．
$$
{\rm minimize}_T\quad D(q(T)||p)
$$
この結果を以下の図にしまします．
特徴としては，以下があります．
- 背景ノイズの推定に便利
- 分布が偏っていないと適用することができない．
以下に，pythonコードを示します．
コード内では，「0log0」を回避するために，KL情報量を以下のような変形を行っていことに注意してください．
$$
D(q(T)||p) = \sum_{i=1}^N q(T) \ln \frac{q(T)}{p}= -\ln M - \sum_{i=1}^N q(T) \ln p \
$$
def InfoScore(data,thresh,bins=100):
num_hist, range_hist = np.histogram(data, bins= bins)
mean_hist = (range_hist[1:] + range_hist[:-1]) / 2
p = num_hist/num_hist.sum()
N = p.shape[0]
M = np.sum(mean_hist>thresh)
if (M== 0):
return np.inf
q = np.zeros(N)
q[mean_hist>thresh] = 1 / M
Dqp = - np.log(M) - np.sum(q*np.log(p))
return Dqp
# Callculation of Otsu score and analyze the optimal
scores_info = np.zeros(256)
for i in range(scores_info.shape[0]):
scores_info[i] = InfoScore(data,i,bins = 190)
thresh_info = np.argmin(scores_info)
混合ガウスモデルのフィッティング
2つのクラスの分布がそれぞれガウス分布であることを仮定し，混合ガウス分布をフィッティンフする手法です．
ガウス分布であることが既知であり，データが十分にあるときのみの使用をオススメします．
2つのガウス分布の交点がしきい値となります．
問題点として：
- ガウスの形状でない＋long tailの場合はうまく推定できない．
- 初期値依存性がある
という問題があります．
以下が結果とpythonコードになります．
pyhtonコードは勉強のためパッケージを使っていませんが，特に偏った心情がなければこちらをおつかいください．
# def gaus func
def gaus(x,mu,sigma2):
return 1/np.sqrt(2 * np.pi * sigma2) *np.exp(-(x-mu)**2/(2*sigma2))
# Class Number
M=2
# Optimmization Condition
OPTIMIZE_EPS = 0.01
# Init
y = data.copy()
mu = 256*np.random.rand(M)
sigma2 = 100*np.random.rand(M)
w = np.random.rand(M)
w = w / w.sum()
for cycle in range(1000):
# E step
gamma_tmp = np.zeros((M,y.shape[0]))
for i in range(M):
gamma_tmp[i] = w[i] * gaus(y,mu[i],sigma2[i])
gamma = gamma_tmp/ gamma_tmp.sum(axis = 0).reshape(1,-1)
Nk = gamma.sum(axis=1)
# M step
mus = (gamma * y).sum(axis = 1)/Nk
sigma2s = (gamma *((y.reshape(1,-1)- mu.reshape(-1,1))**2)).sum(axis = 1)/Nk
ws = Nk / Nk.sum()
# check break condition
if (np.linalg.norm(mu-mus)<OPTIMIZE_EPS)& (np.linalg.norm(sigma2-sigma2s)<OPTIMIZE_EPS) & (np.linalg.norm(w-ws)<OPTIMIZE_EPS):
break
# updata
mu = mus
sigma2 = sigma2s
w = ws
print(cycle)
# make distribution
x = np.arange(0,256,1)
p = np.zeros((M,x.shape[0]))
for i in range(M):
p[i] = w[i] * gaus(x,mu[i],sigma2[i])
# find threshold
if m[0]<m[1]:
thresh_gaus = np.where(p[0]<p[1])[0].min()
else:
thresh_gaus = np.where(p[0]>p[1])[0].min()
まとめ
上記の4つの手法を用いたしきい値は下の表の用になりました．
| K-means | Otsu | Sezan | KL divergence | Gauss Fitting | |
|---|---|---|---|---|---|
| Threshold | 109.0 | 109.0 | 90.558594 | 145.0 | 147.0 |
また，2値化画像は下の図のようになります．
どれがいいとは一概に良い言えませんので，状況に合わせて使い分けましょう!
とりあえず，この記事の上からオススメします．
いくつかの論文は[5]を参考にしています．
よくまとまっていていい感じです．
Code詳細
Reference
[1] H, Steinhaus,“Sur la division des corps matériels en parties” (French). Bull. Acad. Polon. Sci. 4 (12): 801–804 (1957).
[2] Nobuyuki Otsu. "A threshold selection method from gray-level histograms". IEEE Trans. Sys. Man. Cyber. 9 (1): 62–66 (1979).
[3] M. I. Sezan, ‘‘A peak detection algorithm and its application to histogram-based image data reduction,’’ Graph. Models Image Process. 29, 47–59(1985).
[4] C. H. Li and C. K. Lee, ‘‘Minimum cross-entropy thresholding,’’ Pattern Recogn. 26, 617–625 (1993).
[5] Sezgin, M. Survey over image thresholding techniques and quantitative performance evaluation. Journal of Electronic Imaging, 13(1), 220(2004).
Author
Yuji Okamoto yuji.0001@gmail.com