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

# ICV=PA∗(MA−M)2+PB∗(MB−M)2
# 每一个阈值将整个直方图分割成两部分
# 两部分各自的平均值成为 MA 和 MB
# A 部分里的像素数占总像素数的比例记作 PA，B部分里的像素数占总像素数的比例记作 PB。
# 整体图片的灰度值的均值为 M。
img = cv2.imread('ostu.jpg', 0)
# flatten() 将数组变成一维
hist, bins = np.histogram(img.flatten(), 256, [0, 256])
# 计算累积分布图
cdf = hist.cumsum()

cdf_normalized = cdf * hist.max() / cdf.max()
# 可以看下直方图
# plt.plot(cdf_normalized, color='b')
# plt.hist(img.flatten(), 256, [0, 256], color='r')
# plt.xlim([0, 256])
# plt.legend(('cdf', 'histogram'), loc='upper left')
# plt.show()

M = img.flatten().mean()
print(M)
n = img.flatten()
print(n.shape)

M = n.mean()
print(n.min(), n.max(), M)

# 手工实现OSTU
ICV = 0
Threshold = 0
for t in np.linspace(n.min(), n.max() - 1):
    # print(t)
    filter_arr = n > t
    newarr = n[filter_arr]
    PA = len(newarr) / len(n)
    MA = newarr.mean()
    filter_arr = n <= t
    newarr = n[filter_arr]
    PB = len(newarr) / len(n)
    MB = newarr.mean()
    # print(PA, MA, PB, MB)

    I = PA * (MA - M) ** 2 + PB * (MB - M) ** 2
    if I > ICV:
        ICV = I
        Threshold = t

# cv2的实现
ret, otsu = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
print(Threshold, ret)

ret, img2 = cv2.threshold(img, Threshold, 255, cv2.THRESH_BINARY)
plt.subplot(1, 2, 1)
plt.imshow(img, 'gray')
plt.xticks([]), plt.yticks([])

plt.subplot(1, 2, 2)
plt.imshow(img2, 'gray')
plt.xticks([]), plt.yticks([])
plt.show()