大津の手法を用いて3値化を行う

Question

実現したいこと

• 大津の2値化を拡張して，3値化を行いたい

前提

Pythonを用いて，大津の手法を実装しています。閾値の決定方法のみ，ライブラリを使わずにクラス間分散を最大にするやり方で行っています。参考論文はこちらにあります。https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4310076

発生している問題

エラーは出ませんが，二つ目の閾値がうまく計算できていないようです。
一つ目の閾値はおそらく正しいです(二値化だけで実装した時と同じになったため)。結果を見たところ，画像を変えても二つ目の閾値がほぼ250付近のかなり高い画素値に設定されてしまっています。

該当のソースコード

import numpy as np
import matplotlib.pyplot as plt
import cv2

def otsu_threshold(image):
    # グレースケール変換
    gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

    # 画像の総画素数
    total_pixels = gray.shape[0] * gray.shape[1]

    # ヒストグラムの計算
    hist = np.histogram(gray, bins=256, range=[0, 256])[0] / total_pixels
    #plt.figure()
    #plt.plot(hist)

    all_sum = np.sum(hist * np.arange(256))
    all_mean = all_sum / total_pixels

    # クラス間分散と閾値の初期化
    max_variance = 0
    threshold1 = 0
    threshold2 = 0

    # 0から255までの閾値を試す
    for t1 in range(256):
        for t2 in range(t1+1, 255):
            # クラス1とクラス2の画素数を計算(wに相当)
            class1_pixels = np.sum(hist[:t1])
            class2_pixels = np.sum(hist[t1:t2])
            class3_pixels = total_pixels - (class1_pixels + class2_pixels)

            # クラス1とクラス2の画素値の総和を計算
            class1_sum = np.sum(hist[:t1] * np.arange(t1))
            class2_sum = np.sum(hist[t1:t2] * np.arange(t1, t2))
            class3_sum = np.sum(hist[t2:] * np.arange(t2, 256))

            # クラス1とクラス2の平均を計算
            class1_mean = class1_sum / class1_pixels if class1_pixels > 0 else 0
            class2_mean = class2_sum / class2_pixels if class2_pixels > 0 else 0
            class3_mean = class3_sum / class3_pixels if class3_pixels > 0 else 0

            # クラス間分散を計算
            #variance = class1_pixels * class2_pixels * ((class1_mean - class2_mean) ** 2)
            variance = class1_pixels * (class1_mean - all_mean)**2 + class2_pixels * (class2_mean - all_mean)**2 + class3_pixels * (class3_mean - all_mean)**2

            # クラス間分散が最大となる閾値を更新
            if variance > max_variance:
                max_variance = variance
                threshold1 = t1
                threshold2 = t2

    print("threshold = {0},{1}".format(threshold1, threshold2))
    # 3値化処理
    #binary_image = np.where(gray > threshold, 255, 0).astype(np.uint8)
    binary_image = np.zeros_like(gray)
    binary_image[gray <= threshold1] = 0
    binary_image[(gray > threshold1) & (gray <= threshold2)] = 127
    binary_image[gray > threshold2] = 255

    return binary_image

# 画像の読み込み
image = cv2.imread("const.jpg")

# 大津の二値化
binary_image = otsu_threshold(image)

# 結果の表示
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.imshow(cv2.cvtColor(image, cv2.COLOR_BGR2RGB))
plt.title('Original Image')
plt.axis('off')

plt.subplot(1, 2, 2)
plt.imshow(binary_image, cmap='gray', vmin=0, vmax=255)
plt.title('Otsu Thresholding')
plt.axis('off')

plt.tight_layout()
plt.show()

試したこと

いろんな画像を入れてみましたが，どれも同じで，二つ目の閾値が非常に高い値になりました。
式の確認などは何度も行ってるので合っているとは思います。

お力添えいただけますと幸いです。よろしくお願い致します。

Comments

[winterboum]

プログラムはきちんと追いかけて居ませんが、
> 一つ目の閾値はおそらく正しいです(二値化だけで実装した時と同じになったため)。
が？？？
。。。。。。。。。。。｜。。。。。。。。。。。。。 ⇐ 2値化
。。。。。。。｜。。。。。。。。。。。｜。。。。。。 ⇐ ３値化

だと思うので、「二値化だけで実装した時よりちいさい」のでは？

[ikedas]

> for t1 in range(256):
> for t2 in range(t1+1, 255):
ですが、
| for t1 in range(255):
| for t2 in range(t1+1, 256):
のほうがよいのではないでしょうか。

それと、0→255の順で調べていますが、255→0でも同じ結果になるのだろうか、というのが気になりました。2値ならなりそうですが。

[winterboum]

３値でも同じにならないとおかしいですよね。
class間の分散を最大にする、のですから。

[ikedas]

そうなのかもしれません。
> However, they are very simple for M = 2 and 3, (...), so that a special method to reduce the search processes is hardly needed. (p.64)
で、"never" ではなく "hardly" と言っているのが気になったもので。

[winterboum]

？？
引用されてるところだけ、、ですと
プロセスを単純化するのは難しい
ですよね？
逆から計算しても同じか というのと関係ある？

[ikedas]

(hardではなくhardlyですがそれはおいといて)
ヒストグラムに複数の山があったらどうなる? とか考えていました。classが4以上ならどこから調べ始めるかで結果が変わることがあるのは想像がつきますが、3ではどうかなー、と。

[winterboum]

え？
山が複数有るときclass が 4 以上だと方向で変わりますか？
境界値を動かすと classの点の数も分散も単純に増減しますよね。傾きは山か谷かで変わりますが。
4classのとき、3つの境界値 x,y,z に対してクラス間分散 をplotしたとき、同じ分散になるのは山ほど出ると思いますが、同じ高さのpeakが複数出てくる、それも 最大高さのpeakが、というのはイメージできない。。。 頭硬くなったな。

あ？
分散は単純な増減にならない？
一番左の境界が右に行くと左のclass １は 点の数は増え、平均値も増える。平均値の差は減る。減り方と増え方のバランスで、分散は増えたり減ったりしうるか。。。。
感覚では捉えきれなくなりました。降参。

Answer 1

histがtotal_pixelsで割られていて正規化されているので、
下記の最後の行は、

class1_pixels = np.sum(hist[:t1])
class2_pixels = np.sum(hist[t1:t2])
class3_pixels = total_pixels - (class1_pixels + class2_pixels)

total_pixelsではなく、1から引かないといけないと思います。

class3_pixels = 1.0 - (class1_pixels + class2_pixels)