[PYTHON] Zusammenfassung der Methoden zur automatischen Ermittlung von Schwellenwerten

In letzter Zeit wurde ich oft nach dem Wunsch gefragt, den Schwellenwert automatisch zu bestimmen, deshalb habe ich ihn zusammengefasst. Das Schwellenwertdesign entspricht einem eindimensionalen unbeaufsichtigten Clustering mit zwei Gruppen. Ich werde. Dieses Mal werde ich fünf klassische und nützliche Methoden am Beispiel der Binarisierung von Bildern zusammenfassen.

** Hinweis: ** Wenn das Verteilungssystem unbekannt ist, gibt es keine korrekte Antwort für die Schwellenwertbestimmung

Vorbereitung

Dieses Mal werde ich dieses Foto entsprechend verwenden.

Das Helligkeitshistogramm dieses Bildes ist wie folgt.

Das Problem besteht darin, ein Binärbild zu erstellen, indem aus dieser Helligkeitsverteilung der optimale Schwellenwert festgelegt wird.

K-means [1] Die bekannteste Methode für unbeaufsichtigtes Lernen ist K-means ([wiki](https://ja.wikipedia.org/wiki/K Durchschnittsmethode)). K-means zielt darauf ab, den Mittelpunkt {m1, ..., mM} jeder Klasse zu finden, der die folgenden Optimierungsprobleme löst.

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

Hier sind {x1, ..., xN} die angegebenen Daten.

Der Mittelpunkt der Klasse {m1, ..., mM)} kann durch Wiederholen der folgenden zwei Schritte ermittelt werden.

Berechnen Sie den Durchschnittswert in der Klasse.
Ordnen Sie die Daten der nächstgelegenen zentralen Klasse neu zu.

K-Mittel können auch bei mehreren Dimensionen angewendet werden, beim Problem der Schwellenwertgestaltung ist jedoch M = 2.

Das Ergebnis dieser Klassifizierung ist in der folgenden Abbildung dargestellt.

Die Merkmale von K-Mitteln sind wie folgt.

Für jede Aktualisierung muss der Abstand der N × M-Reihenfolge berechnet werden.
Es besteht eine Anfangswertabhängigkeit.

import cv2
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
#Daten lesen
img = cv2.imread('img/main.png',0)
data = img.reshape(-1)

#Etiketteninitialisierung
labels = np.random.randint(0,2,data.shape[0])

#Ausgangsbedingungen
OPTIMIZE_EPSILON = 1

m_0_old = -np.inf
m_1_old = np.inf

for i in range(1000):
#Berechnung jedes Durchschnitts
    m_0 = data[labels==0].mean()
    m_1 = data[labels==1].mean()
#Neuberechnung des Etiketts
    labels[np.abs(data-m_0) < np.abs(data-m_1)] = 0
    labels[np.abs(data-m_0) >= np.abs(data-m_1)] = 1
#Ausgangsbedingungen
    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
#Da sich die Klasse abhängig vom Anfangswert ändert, wird die Klasse mit der kleineren Obergrenze übernommen.
thresh_kmeans = np.minimum(data[labels==0].max(),data[labels==1].max())

Otsu-Methode [2]

Die Otsu-Methode ist eine Methode zur Maximierung der Differenz zwischen den gewichteten Durchschnittswerten jeder Klasse (wiki). Insbesondere entspricht es dem Finden des Schwellenwerts T, der das folgende Optimierungsproblem löst.

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

Hier, $ 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} $ Es wird.

Das Problem der Maximierung von σb oben ist die gewichtete Summe der Varianzen jeder Klasse.

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

In Übereinstimmung mit dem Problem der Minimierung. Daher können wir das Problem der Minimierung von σa und das Problem der Minimierung von σa / σb lösen.

Das Ergebnis dieser Klassifizierung ist in der folgenden Abbildung dargestellt.

Der Code in Python ist unten geschrieben.

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

Methode von Sezan et al. [3]

Dies ist eine Schwellenwertentwurfsmethode, die die Tatsache ausnutzt, dass es zwei Spitzen in der Verteilung der Bildhelligkeit gibt. Entwerfen Sie den Schwellenwert, indem Sie den Wert des inneren Randes der beiden Spitzen ermitteln.

Zunächst wird das geglättete Histogramm für die Polarberechnung verwendet, um den Peak ganz links (kleiner Wert) und den Peak ganz rechts (großer Wert) zu erfassen und als [m0, m1] festzulegen. In ähnlicher Weise werden der linke und der rechte Saum jedes Peaks durch die Polarberechnung berechnet, und der linke Saum ist [s0, s1] und der rechte Saum ist [e0, e1].

Zu diesem Zeitpunkt wird angenommen, dass der Schwellenwert zwischen dem rechten Saum e0 des linken Peaks und dem linken Saum s1 des rechten Peaks liegt. Daher wird der Schwellenwert unter Verwendung des Entwurfsparameters γ von 0 oder mehr und 1 oder weniger entworfen.

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

Das Ergebnis dieser Klassifizierung ist in der folgenden Abbildung dargestellt.

Die Merkmale der binären Methode von Sezan et al. Sind wie folgt.

Es gibt viele Hyperparameter, die angepasst werden müssen, aber Sie können sie je nach Situation richtig verwenden.
Schwierig, es mehrdimensional zu machen

Der Python-Code sieht folgendermaßen aus:

#Glättungsparameter
sigma = 5.0

#Prozentsatz der Wichtigkeit
# gamma = 1 :Nehmen Sie einen weiten schwarzen Bereich
# gamma = 0 :Nehmen Sie eine große weiße Fläche
gamma = 0.5
#Gaußscher Kern zum Glätten
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))
#Stellen Sie so ein, dass die Summe 1 ist
    G_kernel = G_kernel/G_kernel.sum()

    return G_kernel
#Erstellen eines Gaußschen Kernels
kernel = getGaus(55,sigma)
#Histogramm
num_hist, range_hist = np.histogram(data, bins= 256)
mean_hist = (range_hist[1:] + range_hist[:-1]) / 2

#Glätten
hist_bar = np.convolve(num_hist,kernel,'same')
#Berechnung der Differenz
d_hist = hist_bar[:-1] - hist_bar[1:]
#Kantenbearbeitung
d_hist = np.r_[[0],d_hist,[0]]

#Peakerkennung
m = np.where((d_hist[1:] >=0) & (d_hist[:-1] <=0))[0]
#Lokale Erkennung
es =np.where((d_hist[1:] <=0) & (d_hist[:-1] >=0))[0]

#Maximale und minimale Spitzen
m0 = m.min()
m1 = m.max()

#Lokalität vor und nach dem Gipfel
s0 = es[es<m0].max()
e0 = es[es>m0].min()
s1 = es[es<m1].max()
e1 = es[es>m1].min()

#Schwellenwertbestimmung
thresh_Sezan = (1 - gamma) * mean_hist[e0] + gamma * mean_hist[s1]

Minimieren Sie die Menge an Informationen über Calvac Libra [4]

Dies ist eine Methode, die die Menge an Informationen zur Calvac Libra (KL) verwendet, die Forschern der Informationstheorie vertraut ist. Es gibt viele Artikel, die auf der Grundlage anderer Informationsmengenstandards als der KL-Informationsmenge untersucht wurden.

Zunächst wird die Wahrscheinlichkeitsverteilung, die durch Normalisieren des Helligkeitshistogramms erhalten wird, als p definiert. Als nächstes wird die Wahrscheinlichkeitsrelationsverteilung q (T), die der durch den Schwellenwert T bestimmten Binärisierung entspricht, wie folgt definiert. $ q(i\leq T) = 0,\quad q(i>T) = 1/M $

Wobei M die Anzahl der Elemente ist, die größer als der Schwellenwert sind. q(T)KL Informationsmenge D von und p(q(T)||p)Der minimierte Wert ist der Schwellenwert. Mit anderen Worten, es löst das folgende Optimierungsproblem. $ {\rm minimize}_T\quad D(q(T)||p) $ Das Ergebnis ist in der folgenden Abbildung dargestellt.

Die Funktionen sind wie folgt.

Praktisch zum Schätzen von Hintergrundgeräuschen
Es kann nur angewendet werden, wenn die Verteilung voreingenommen ist.

Unten ist der Python-Code. Beachten Sie, dass im Code die KL-Informationsmenge wie folgt transformiert wird, um "0log0" zu vermeiden.

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)

Gemischte Gaußsche Modellanpassung

Dies ist eine Methode zum Anpassen einer gemischten Gaußschen Verteilung, vorausgesetzt, die Verteilungen der beiden Klassen sind Gaußsche Verteilungen. Es ist bekannt, dass es eine Gaußsche Verteilung hat, und es wird empfohlen, es nur zu verwenden, wenn genügend Daten vorhanden sind.

Der Schnittpunkt der beiden Gaußschen Verteilungen ist die Schwelle.

Als Problem:

Wenn es keine Gaußsche Form + langer Schwanz ist, kann es nicht gut geschätzt werden.
Es besteht eine Anfangswertabhängigkeit

Es gibt ein Problem.

Unten finden Sie das Ergebnis und den Python-Code. Der Pyhton-Code verwendet keine Pakete zum Lernen. Wenn Sie jedoch keine besondere Tendenz haben, gehen Sie zu hier. Bitte benutzen Sie es.

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

Zusammenfassung

Die Schwellenwerte unter Verwendung der obigen vier Methoden sind in der folgenden Tabelle aufgeführt.

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

Das binärisierte Bild ist in der folgenden Abbildung dargestellt.

Ich kann nicht sagen, welches besser ist, also lasst es uns je nach Situation richtig verwenden! Vorerst empfehle ich es von oben in diesem Artikel.

Einige Artikel beziehen sich auf [5]. Es fühlt sich gut an, gut organisiert zu sein.

Codedetails

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]