[PYTHON] Résumé des méthodes pour déterminer automatiquement les seuils

Récemment, on m'a souvent posé des questions sur le désir de déterminer automatiquement le seuil, alors je l'ai résumé. La conception du seuil correspond à un regroupement à deux groupes unidimensionnel non supervisé. Je vais. Cette fois, je résumerai cinq méthodes classiques et utiles en utilisant la binarisation d'images à titre d'exemple.

** Remarque: ** Si le système de distribution est inconnu, il n'y a pas de réponse correcte pour la détermination du seuil

Préparation

Cette fois, j'utiliserai cette photo prise de manière appropriée.

L'histogramme de luminosité de cette image est le suivant.

Le problème est de créer une image binaire en définissant le seuil optimal à partir de cette distribution de luminosité.

K-means [1] La méthode la plus connue pour l'apprentissage non supervisé est K-means ([wiki](https://ja.wikipedia.org/wiki/K méthode moyenne)). K-means vise à trouver le point central {m1, ..., mM} de chaque classe qui résout les problèmes d'optimisation suivants.

{\rm minimize}\_{m_1,\ldots, m_M} \quad \sum_{i=1}^N \min_j |x_i - m_j|

Ici {x1, ..., xN} sont les données données.

Le point central de la classe {m1, ..., mM)} peut être trouvé en répétant les deux étapes suivantes.

Calculez la valeur moyenne de la classe.
Réaffectez les données à la classe centrale la plus proche.

K-means peut être appliqué même dans le cas de dimensions multiples, mais dans le problème de la conception de seuil, M = 2.

Le résultat de cette classification est illustré dans la figure suivante.

Les caractéristiques de K-means sont les suivantes. --Il est nécessaire de calculer la distance d'ordre N × M pour chaque mise à jour. --Il existe une dépendance de valeur initiale.

import cv2
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
#Lire les données
img = cv2.imread('img/main.png',0)
data = img.reshape(-1)

#Initialisation de l'étiquette
labels = np.random.randint(0,2,data.shape[0])

#Conditions de sortie
OPTIMIZE_EPSILON = 1

m_0_old = -np.inf
m_1_old = np.inf

for i in range(1000):
#Calcul de chaque moyenne
    m_0 = data[labels==0].mean()
    m_1 = data[labels==1].mean()
#Recalcul des étiquettes
    labels[np.abs(data-m_0) < np.abs(data-m_1)] = 0
    labels[np.abs(data-m_0) >= np.abs(data-m_1)] = 1
#Conditions de sortie
    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
#Puisque la classe change en fonction de la valeur initiale, celle avec la limite supérieure la plus petite est adoptée.
thresh_kmeans = np.minimum(data[labels==0].max(),data[labels==1].max())

Méthode Otsu [2]

La méthode Otsu est une méthode pour maximiser la différence entre les valeurs moyennes pondérées de chaque classe (wiki). Plus précisément, cela correspond à trouver le seuil T qui résout le problème d'optimisation suivant.

{\rm maximize}\_T\quad \sigma_b^2 (T):= w_1 (m_1 - m)^2 +w_2 (m_2 - m)^2

ici, $ 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} $ Il devient.

Le problème de la maximisation de σb ci-dessus est la somme pondérée des variances de chaque classe.

\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}

Conforme au problème de la minimisation. Par conséquent, nous pouvons résoudre le problème de la minimisation de σa et le problème de la minimisation de σa / σb.

Le résultat de cette classification est illustré dans la figure suivante.

Le code en python est écrit ci-dessous.

# 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)

Méthode de Sezan et al. [3]

Il s'agit d'une méthode de conception de seuil qui utilise le fait qu'il existe deux pics dans la distribution de la luminosité de l'image. Concevez le seuil en trouvant la valeur de la jupe intérieure des deux pics.

Tout d'abord, l'histogramme lissé est utilisé pour le calcul polaire pour détecter le pic le plus à gauche (petite valeur) et le pic le plus à droite (grande valeur), et les définir comme [m0, m1], respectivement. De même, l'ourlet gauche et droit de chaque pic sont calculés par le calcul polaire, et l'ourlet gauche est [s0, s1] et l'ourlet droit est [e0, e1].

A ce moment, la valeur de seuil est considérée comme étant entre l'ourlet droit eO du pic gauche et l'ourlet gauche s1 du pic droit. Par conséquent, le seuil est conçu en utilisant le paramètre de conception γ de 0 ou plus et de 1 ou moins.

{\rm Sezan~threshold} := (1 - \gamma )e_0 + \gamma s_1

Le résultat de cette classification est illustré dans la figure suivante.

Les caractéristiques de la méthode binaire de Sezan et al. Sont les suivantes.

Il existe de nombreux hyper paramètres à régler, mais vous pouvez les utiliser correctement en fonction de la situation.
Difficile de le rendre multidimensionnel

Le code python ressemble à ceci:

#Paramètres de lissage
sigma = 5.0

#Pourcentage d'importance
# gamma = 1 :Prenez une large zone de noir
# gamma = 0 :Prenez une large zone de blanc
gamma = 0.5
#Noyau gaussien pour le lissage
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))
#Ajustez pour que la somme soit 1
    G_kernel = G_kernel/G_kernel.sum()

    return G_kernel
#Créer un noyau gaussien
kernel = getGaus(55,sigma)
#Histogramme
num_hist, range_hist = np.histogram(data, bins= 256)
mean_hist = (range_hist[1:] + range_hist[:-1]) / 2

#Lissage
hist_bar = np.convolve(num_hist,kernel,'same')
#Calcul de la différence
d_hist = hist_bar[:-1] - hist_bar[1:]
#Traitement des bords
d_hist = np.r_[[0],d_hist,[0]]

#Détection de pic
m = np.where((d_hist[1:] >=0) & (d_hist[:-1] <=0))[0]
#Détection locale
es =np.where((d_hist[1:] <=0) & (d_hist[:-1] >=0))[0]

#Pics maximum et minimum
m0 = m.min()
m1 = m.max()

#Localité avant et après le pic
s0 = es[es<m0].max()
e0 = es[es>m0].min()
s1 = es[es<m1].max()
e1 = es[es>m1].min()

#Détermination du seuil
thresh_Sezan = (1 - gamma) * mean_hist[e0] + gamma * mean_hist[s1]

Minimisez la quantité d'informations sur Calvac Libra [4]

Il s'agit d'une méthode qui utilise la quantité d'informations Calvac Libra (KL) familière aux chercheurs en théorie de l'information. Il existe de nombreux articles qui ont été étudiés sur la base de normes de quantité d'information autres que la quantité d'information KL.

Premièrement, la distribution de probabilité obtenue en normalisant l'histogramme de luminosité est définie comme p. Ensuite, la distribution de relation de probabilité q (T) correspondant à la binarisation déterminée par la valeur de seuil T est définie comme suit. $ q(i\leq T) = 0,\quad q(i>T) = 1/M $

Où M est le nombre d'éléments supérieur au seuil. q(T)KL information montant D par et p(q(T)||p)La valeur qui minimise est la valeur seuil. En d'autres termes, il résout le problème d'optimisation suivant. $ {\rm minimize}_T\quad D(q(T)||p) $ Le résultat est illustré dans la figure ci-dessous.

Les caractéristiques sont les suivantes.

Pratique pour estimer le bruit de fond --Il ne peut être appliqué que si la distribution est biaisée.

Vous trouverez ci-dessous le code python. Notez que dans le code, le montant de l'information KL est transformé comme suit afin d'éviter "0log0".

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)

Raccord de modèle gaussien mixte

Il s'agit d'une méthode d'ajustement d'une distribution gaussienne mixte, en supposant que les distributions des deux classes sont des distributions gaussiennes. Il est connu pour avoir une distribution gaussienne et son utilisation n'est recommandée que lorsque les données sont suffisantes.

L'intersection des deux distributions gaussiennes est le seuil.

En tant que problème: --Si ce n'est pas une forme gaussienne + longue queue, il ne peut pas être bien estimé. --Il existe une dépendance de valeur initiale

Il ya un problème.

Vous trouverez ci-dessous le résultat et le code python. Le code pyhton n'utilise pas de package pour étudier, mais à moins que vous n'ayez un biais particulier, allez sur ici. Veuillez l'utiliser.

# 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()

Résumé

Les seuils utilisant les quatre méthodes ci-dessus sont indiqués dans le tableau ci-dessous.

	K-means	Otsu	Sezan	KL divergence	Gauss Fitting
Threshold	109.0	109.0	90.558594	145.0	147.0

L'image binarisée est illustrée dans la figure ci-dessous.

Je ne peux pas dire lequel est le meilleur, alors utilisons-le correctement en fonction de la situation! Pour le moment, je le recommande depuis le haut de cet article.

Certains articles font référence à [5]. Ça fait du bien d'être bien organisé.

Détails du code

https://github.com/yuji0001/Threshold_Technic

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 [email protected]