Numpyで画像の二値化を実装してみます。
まず、二値化する画像を読み込み、グレースケール画像に変換します。
original_image = plt.imread(image_name)
if np.issubdtype(original_image.dtype, np.floating):
original_image = (original_image * 255).astype(np.uint8)
gray_image = (0.2116 * original_image[:,:,0] + 0.7152 * original_image[:,:,1] + 0.0722 * original_image[:,:,2]).astype(np.uint8)
plt.imshow(gray_image, cmap='gray')
任意の閾値で二値化をする関数を実装してみます。これはnumpy.whereを用いると簡単に実装できます。
def thresholding(image, threshold):
return np.where(image <= threshold, 0, 255)
threshold_image = thresholding(gray_image, 128)
plt.imshow(threshold_image, cmap='gray')
次に大津の二値化(判別分析法)を実装してみます。
大津の二値化ではクラス間分散$\sigma_b^2$とクラス内分散$\sigma_w^2$の比である分離度$\frac{\sigma_b^2}{\sigma_w^2}$が最大になる閾値で二値化を行います。
黒画素クラスの平均と分散を$m_1$と$\sigma_1^2$、白画素クラスの平均と分散を$m_2$と$\sigma_2^2$とし、黒画素クラスと白画素クラスの所属する画素数をそれぞれ$\omega_1$、$\omega_2$とすると、クラス間分散$\sigma_b^2$とクラス内分散$\sigma_w^2$はそれぞれ以下のようになります。
\sigma_w^2=\frac{\omega_1\sigma_1^2+\omega_2\sigma_2^2}{\omega_1\omega_2} \\
\sigma_b^2=\frac{\omega_1\omega_2(m_1-m_2)^2}{(\omega_1+\omega_2)^2}
全分散$\sigma_t^2$は閾値によらず一定で、クラス間分散$\sigma_b^2$とクラス内分散$\sigma_w^2$と以下のように関係にあります。
\sigma_t^2=\sigma_w^2 + \sigma_b^2
このことから分離度$\frac{\sigma_b^2}{\sigma_w^2}$は次のようになり、クラス間分散$\sigma_b^2$が最大のときに分離度$\frac{\sigma_b^2}{\sigma_w^2}$も最大になることがわかります。
\frac{\sigma_b^2}{\sigma_w^2}=\frac{\sigma_b^2}{\sigma_t^2-\sigma_b^2}
実装としては、すべての閾値に対してクラス間分散$\sigma_b^2$を計算して、クラス間分散$\sigma_b^2$が最大となるとき閾値で二値化を行います。クラス間分散$\sigma_b^2$の分母は一定なので、実装では計算を省略します。
def otsu_thresholding(image):
between_class_variance = np.zeros(256)
for i in range(np.min(image), np.max(image)):
class1 = image <= i
class2 = ~class1
average1 = np.average(image[class1])
average2 = np.average(image[class2])
num1 = np.count_nonzero(class1)
num2 = np.count_nonzero(class2)
between_class_variance[i] = (num1 * num2 * (average1 - average2) ** 2)
threshold = np.argmax(between_class_variance)
return np.where(image <= threshold, 0, 255)
plt.imshow(otsu_thresholding(gray_image), cmap='gray')
実装したコードはGoogle Colaboratoryに置いてあります。