从头开始python实现Otsu二值化

Question

看来我的实现不正确，不确定我到底做错了什么：

这是我的图像的直方图：

所以阈值应该在 170 左右？我将阈值设为 130。

这是我的代码：

#Otsu in Python

import numpy as np
from PIL import Image
import matplotlib.pyplot as plt  


def load_image(file_name):
    img = Image.open(file_name)
    img.load()
    bw = img.convert('L')
    bw_data = np.array(bw).astype('int32')
    BINS = np.array(range(0,257))
    counts, pixels =np.histogram(bw_data, BINS)
    pixels = pixels[:-1]
    plt.bar(pixels, counts, align='center')
    plt.savefig('histogram.png')
    plt.xlim(-1, 256)
    plt.show()

    total_counts = np.sum(counts)
    assert total_counts == bw_data.shape[0]*bw_data.shape[1]

    return BINS, counts, pixels, bw_data, total_counts

def within_class_variance():
    ''' Here we will implement the algorithm and find the lowest Within-  Class Variance:

        Refer to this page for more details http://www.labbookpages.co.uk
/software/imgProc/otsuThreshold.html'''

    for i in range(1,len(BINS), 1):         #from one to 257 = 256 iterations
       prob_1 =    np.sum(counts[:i])/total_counts
       prob_2 = np.sum(counts[i:])/total_counts
       assert (np.sum(prob_1 + prob_2)) == 1.0



       mean_1 = np.sum(counts[:i] * pixels[:i])/np.sum(counts[:i])
       mean_2 = np.sum(counts[i:] * pixels[i:] )/np.sum(counts[i:])
       var_1 = np.sum(((pixels[:i] - mean_1)**2 ) * counts[:i])/np.sum(counts[:i])
       var_2 = np.sum(((pixels[i:] - mean_2)**2 ) * counts[i:])/np.sum(counts[i:])


       if i == 1:
         cost = (prob_1 * var_1) + (prob_2 * var_2)
         keys = {'cost': cost, 'mean_1': mean_1, 'mean_2': mean_2, 'var_1': var_1, 'var_2': var_2, 'pixel': i-1}
         print('first_cost',cost)


       if (prob_1 * var_1) +(prob_2 * var_2) < cost:
         cost =(prob_1 * var_1) +(prob_2 * var_2)
         keys = {'cost': cost, 'mean_1': mean_1, 'mean_2': mean_2, 'var_1': var_1, 'var_2': var_2, 'pixel': i-1}  #pixels is i-1 because BINS is starting from one

    return keys







if __name__ == "__main__":

    file_name = 'fish.jpg'
    BINS, counts, pixels, bw_data, total_counts =load_image(file_name)
    keys =within_class_variance()
    print(keys['pixel'])
    otsu_img = np.copy(bw_data).astype('uint8')
    otsu_img[otsu_img > keys['pixel']]=1
    otsu_img[otsu_img < keys['pixel']]=0
    #print(otsu_img.dtype)
    plt.imshow(otsu_img)
    plt.savefig('otsu.png')
    plt.show()

生成的 otsu 图像如下所示：

这是鱼的图像（它有一个赤膊男子抱着一条鱼，所以工作可能不安全）：

链接：https://i.stack.imgur.com/EDTem.jpg

编辑：

原来通过将阈值改为255（差异更明显）

Answer 1

我在发布的答案中使用了实现@Jose A，它试图最大化类间方差。看起来 jose 忘记将强度级别乘以它们各自的强度像素计数（为了计算平均值），所以我更正了背景平均 mub 和前景平均 muf 的计算。我将此作为答案发布，并尝试编辑已接受的答案。

def otsu(gray):
    pixel_number = gray.shape[0] * gray.shape[1]
    mean_weight = 1.0/pixel_number
    his, bins = np.histogram(gray, np.arange(0,257))
    final_thresh = -1
    final_value = -1
    intensity_arr = np.arange(256)
    for t in bins[1:-1]: # This goes from 1 to 254 uint8 range (Pretty sure wont be those values)
        pcb = np.sum(his[:t])
        pcf = np.sum(his[t:])
        Wb = pcb * mean_weight
        Wf = pcf * mean_weight

        mub = np.sum(intensity_arr[:t]*his[:t]) / float(pcb)
        muf = np.sum(intensity_arr[t:]*his[t:]) / float(pcf)
        #print mub, muf
        value = Wb * Wf * (mub - muf) ** 2

        if value > final_value:
            final_thresh = t
            final_value = value
    final_img = gray.copy()
    print(final_thresh)
    final_img[gray > final_thresh] = 255
    final_img[gray < final_thresh] = 0
    return final_img

Answer 2

我不知道我的实现是否正常。但这就是我得到的：

def otsu(gray):
    pixel_number = gray.shape[0] * gray.shape[1]
    mean_weigth = 1.0/pixel_number
    his, bins = np.histogram(gray, np.array(range(0, 256)))
    final_thresh = -1
    final_value = -1
    for t in bins[1:-1]: # This goes from 1 to 254 uint8 range (Pretty sure wont be those values)
        Wb = np.sum(his[:t]) * mean_weigth
        Wf = np.sum(his[t:]) * mean_weigth

        mub = np.mean(his[:t])
        muf = np.mean(his[t:])

        value = Wb * Wf * (mub - muf) ** 2

        print("Wb", Wb, "Wf", Wf)
        print("t", t, "value", value)

        if value > final_value:
            final_thresh = t
            final_value = value
    final_img = gray.copy()
    print(final_thresh)
    final_img[gray > final_thresh] = 255
    final_img[gray < final_thresh] = 0
    return final_img

Answer 3

这是我刚刚从scikit image source code 修改的另一个实现。它是为一维数组设计的，因此您必须编写一个包装器才能使其与图像一起使用。

def threshold_otsu(x: Iterable, *args, **kwargs) -> float:
    """Find the threshold value for a bimodal histogram using the Otsu method.

    If you have a distribution that is bimodal (AKA with two peaks, with a valley
    between them), then you can use this to find the location of that valley, that
    splits the distribution into two.

    From the SciKit Image threshold_otsu implementation:
    https://github.com/scikit-image/scikit-image/blob/70fa904eee9ef370c824427798302551df57afa1/skimage/filters/thresholding.py#L312
    """
    counts, bin_edges = np.histogram(x, *args, **kwargs)
    bin_centers = (bin_edges[1:] + bin_edges[:-1]) / 2

    # class probabilities for all possible thresholds
    weight1 = np.cumsum(counts)
    weight2 = np.cumsum(counts[::-1])[::-1]
    # class means for all possible thresholds
    mean1 = np.cumsum(counts * bin_centers) / weight1
    mean2 = (np.cumsum((counts * bin_centers)[::-1]) / weight2[::-1])[::-1]

    # Clip ends to align class 1 and class 2 variables:
    # The last value of ``weight1``/``mean1`` should pair with zero values in
    # ``weight2``/``mean2``, which do not exist.
    variance12 = weight1[:-1] * weight2[1:] * (mean1[:-1] - mean2[1:]) ** 2

    idx = np.argmax(variance12)
    threshold = bin_centers[idx]
    return threshold