十三种基于直方图的图像全局二值化算法实现

  • 1. 什么是基于直方图的图像全局二值化算法
  • 2. 灰度平均值
  • 3. 百分比阈值（P-Tile法）
  • 3. 基于双峰的阈值
    • 3.1 基于平均值的阈值
    • 3.2 基于最小值的阈值方法
  • 4. 迭代最佳阈值
  • 5.大津(OTSU)法
  • 6.一维最大熵
  • 7.力矩保持法
  • 8.模糊集
  • 9. Kittler最小错误分类法
  • 10. ISODATA
  • 11. Shanbhag 法
  • 12. Yen法

1. 什么是基于直方图的图像全局二值化算法

本文的内容来自与十三种基于直方图的图像全局二值化算法原理、实现、代码及效果。本文的中的算法是基于图像的灰度直方图进行阈值分割。
所谓图像的灰度直方图是在将图像进行灰度处理后，计算灰度值的分布。也就是每个灰度值在灰度图像中有多少个点。一般情况下，图像灰度值的取值范围为0~255，所以灰度直方图是一个有256个元素的数组。
阈值分割就是找到一个合适阈值，对灰度图像进行处理。大于阈值的设置为255，小于阈值的设置为0。也可以根据需要设为其他数值。这样处理灰度图像，可以将前景与背景做一个区分。
此方法的关键就是计算阈值。在前面的链接中，介绍了13种方法并提供了C语言实现。本文将这些算法用Python进行了重写。

2. 灰度平均值

灰度平均值是图像总像素值/总像素数：
总像素值=∑g=0255g×h(g)总像素值=\displaystyle \sum_{g=0}^{255}g\times h(g)总像素值=g=0∑255​g×h(g)
总像素数=∑g=0255h(g)总像素数=\displaystyle \sum_{g=0}^{255}h(g)总像素数=g=0∑255​h(g)
总像素数其实就是图像的宽度X图像的长度，也是图像的大小。

# coding:utf8import numpy as np
import cv2
from matplotlib import pyplot as pltdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef GetMeanThreshold(H):Amount = np.sum(H)Gray_Sum = np.sum(H*np.arange(256))return int(float(Gray_Sum/Amount))def GrayThreshold(image, maxval=255):g = GrayHist(image)thresh = GetMeanThreshold(g)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th, img_new = GrayThreshold(img_gray)th1, img_new_1 = cv2.threshold(img_gray, 146, 255, cv2.THRESH_BINARY)print(th, th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()


对于小鸟的图像其灰度平均阈值为146。

3. 百分比阈值（P-Tile法）

此方法先设定一个阈值p(先验概率)，比如50%。然后从0开始对灰度值求和。计算到灰度值Y时，[0~Y]的综合>=总像素值*p，则设定灰度值Y为阈值。
该方法需要根据先验概率确定一个百分比。人的经验很重要。

# coding:utf8import numpy as np
import cv2
from matplotlib import pyplot as pltdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef GetPTileThreshold(H,ptile=0.5):Amount = np.sum(H)A = Amount * ptilepsum = 0for Y in range(256):psum += H[Y]if psum >= A:return Yreturn -1  # 没有符合条件的阈值def GrayThreshold(image, maxval=255):g = GrayHist(image)thresh = GetPTileThreshold(g,0.25)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th, img_new = GrayThreshold(img_gray)th1, img_new_1 = cv2.threshold(img_gray, 127, 255, cv2.THRESH_BINARY)print(th, th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()

当百分比设为25%时，阈值为123。这种方法需要认为调节阈值，否则很难得到理想的结果。

3. 基于双峰的阈值

此方法适合与有明显双峰值的图像，对于单峰值或者平坦的灰度直方图形并不适合。
其原理是通过迭代对直方图数据进行判断，判端是否是一个双峰的直方图。如果不是，则对直方图数据进行窗口为3的平滑。如果迭代了一定的数量比如1000次后仍未获得一个双峰的直方图，则函数执行失败。如成功获得，则最终阈值取两个双峰之间的谷底值作为阈值。
所谓窗口为3的平滑，就是将灰度p的值，用p-1，p和p+1的值进行平均。
直方图通过迭代形成双峰后，有2种方法获得阈值：
（1）双峰的平均值
（2）双峰之间的谷底值
第一种方法是找到双峰后，取平均值的整数。第二种方法找到谷底的最小值。

3.1 基于平均值的阈值

# coding:utf8
import numpy as np
import cv2
from matplotlib import pyplot as pltdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef GetIntermodesThreshold(H,method=0):Iter = 0HistGramC = np.array(H, dtype=np.float64)  # 基于精度问题，一定要用浮点数来处理，否则得不到正确的结果HistGramCC = np.array(H, dtype=np.float64)  # 求均值的过程会破坏前面的数据，因此需要两份数据# 通过三点求均值来平滑直方图while IsDimodal(HistGramCC) == False:  # 判断是否已经是双峰的图像了HistGramCC[0] = (HistGramC[0] + HistGramC[0] + HistGramC[1]) / 3  # 第一点for Y in range(1, 255):HistGramCC[Y] = (HistGramC[Y - 1] + HistGramC[Y] + HistGramC[Y + 1]) / 3  # 中间的点HistGramCC[255] = (HistGramC[254] + HistGramC[255] + HistGramC[255]) / 3  # 最后一点HistGramC = np.array(HistGramCC, dtype=np.float64)  # 备份数据Iter += 1if Iter >= 1000:return -1  # 直方图无法平滑为双峰的，返回错误代码if method > 0:Peakfound = Falsefor Y in range(1,255):if HistGramCC[Y - 1] < HistGramCC[Y] and HistGramCC[Y + 1] < HistGramCC[Y]:Peakfound = Trueif Peakfound and HistGramCC[Y - 1] >= HistGramCC[Y] and HistGramCC[Y + 1] >= HistGramCC[Y]:return Y - 1return -1else:# 阈值为两峰值的平均值Peak = np.zeros(2,dtype=np.uint16)Index = 0for Y in range(1, 255):if HistGramCC[Y - 1] < HistGramCC[Y] and HistGramCC[Y + 1] < HistGramCC[Y]:Peak[Index] = Y - 1Index += 1return int((Peak[0] + Peak[1]) / 2)def IsDimodal(H):  # 检测直方图是否为双峰的# 对直方图的峰进行计数，只有峰数位2才为双峰Count = 0for Y in range(1, 255):if H[Y - 1] < H[Y] and H[Y + 1] < H[Y]:Count += 1if Count > 2: return Falseif Count == 2:return Trueelse:return Falsedef GrayThreshold(image, maxval=255):g = GrayHist(img_gray)thresh = GetIntermodesThreshold(g,1)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th1, img_new_1 = cv2.threshold(img_gray, 0, 255, cv2.THRESH_TRIANGLE)th, img_new = GrayThreshold(img_gray)print(th, th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()

此种方法得到的小鸟图像阈值为159。原始直方图和平滑后的直方图（黄色）如图：

3.2 基于最小值的阈值方法

该方法是查找双峰之间的最小值。使用此方法需要调用GetIntermodesThreshold(g,1)。最后一个参数设置为1即可。

4. 迭代最佳阈值

此方法是先从左边找到HistGram中第一个数值不为0的点min，再从右边找到第一个数值不为0的点max。然后取阈值threshold为此2点的中间值。通过计算[min,threshold]和(threshold,max]平均值灰度值对应的点，进行迭代直到平均灰度值点与阈值点相等。
平均灰度值计算为：
Gmin=∑g=minTkg×h(g)∑g=minTkh(g)G_{min}=\frac{\sum_{g=min}^{T_k}g\times h(g)}{\sum_{g=min}^{T_k}h(g)}Gmin​=∑g=minTk​​h(g)∑g=minTk​​g×h(g)​
Gmax=∑g=Tkmaxg×h(g)∑g=Tkmaxh(g)G_{max}=\frac{\sum_{g=T_k}^{max}g\times h(g)}{\sum_{g=T_k}^{max}h(g)}Gmax​=∑g=Tk​max​h(g)∑g=Tk​max​g×h(g)​
Tk+1=Gmin+Gmax2T_{k+1}=\frac{G_{min}+G_{max}}{2}Tk+1​=2Gmin​+Gmax​​
重复进行运算，直到Tk=Tk+1。
需要注意的是可能存在不收敛的情况。

# coding:utf8import numpy as np
import cv2
from matplotlib import pyplot as pltdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef GetIterativeBestThreshold(HistGram):Iter = 0value = np.where(HistGram > 0)[0]  # 非0下标index = len(value)if index == 0:return 0elif index == 1 or index == 2:return value[0]else:MinValue = value[0]MaxValue = value[-1]Threshold = MinValueNewThreshold = int(MaxValue + MinValue) >> 1while Threshold != NewThreshold:  # 当前后两次迭代的获得阈值相同时，结束迭代Threshold = NewThreshold#[MinValue,Threshold]灰度平均值对应的点SUmInteralOne = np.sum(HistGram[MinValue:Threshold + 1] * np.arange(MinValue, Threshold + 1))SumOne = np.sum(HistGram[MinValue:Threshold + 1])MeanValueOne = SUmInteralOne / SumOne# [Threshold+1,MaxValue]灰度平均值对应的点SUmInteralTwo = np.sum(HistGram[Threshold + 1:MaxValue + 1] * np.arange(Threshold + 1, MaxValue + 1))SumTwo = np.sum(HistGram[Threshold + 1:MaxValue + 1])MeanValueTwo = SUmInteralTwo / SumTwoNewThreshold = int(MeanValueOne + MeanValueTwo) >> 1  # 求出新的阈值Iter += 1if Iter >= 1000:return -1return Thresholddef GrayThreshold(image, maxval=255):g = GrayHist(image)thresh = GetIterativeBestThreshold(g)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th, img_new = GrayThreshold(img_gray)th1, img_new_1 = cv2.threshold(img_gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_TRIANGLE)print(th, th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()

此方法求出的阈值为126。

5.大津(OTSU)法

大津法是日本学者大津在1979年提出的方法，也叫最大类间法。其对单峰值或者比较平坦的灰度值有很好的效果，对双峰或者多峰的情况则较差。其方法如下：
记t为前景与背景的分割阈值，前景点数占图像比例为w0，平均灰度为u0；背景点数占图像比例为w1，平均灰度为u1。
则图像的总平均灰度为：u=w0 *u0+w1*u1。
前景和背景图象的方差：g=w0*(u0-u)*(u0-u)+w1*(u1-u)*(u1-u)=w0*w1*(u0-u1)2,此公式为方差公式。
不断的迭代求解出最大的方差值，可以认为此时前景和背景的相差最大，从而得到分割图像的阈值。

# coding:utf8import numpy as np
import cv2
from matplotlib import pyplot as pltdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef GetOSTUThreshold(H):# 判断只有1种灰度、2种灰度以及多种灰度V = np.where(H > 0)[0]  # 非0下标I = len(V)if I == 0: return 0if I == 1: return V[0]if I == 2: return V[0] if H[V[0]] < H[V[1]] else V[1]MinValue = V[0]MaxValue = V[-1]PixelBack = 0PixelIntegralBack = 0Threshold = 0Amount = np.sum(H)PixelIntegral = np.sum(H * np.arange(256))SigmaB = -1for Y in range(MinValue, MaxValue):PixelBack = PixelBack + H[Y]PixelFore = Amount - PixelBackOmegaBack = PixelBack / AmountOmegaFore = PixelFore / AmountPixelIntegralBack += H[Y] * YPixelIntegralFore = PixelIntegral - PixelIntegralBackMicroBack = PixelIntegralBack / PixelBackMicroFore = PixelIntegralFore / PixelForeSigma = OmegaBack * OmegaFore * (MicroBack - MicroFore) * (MicroBack - MicroFore)if Sigma > SigmaB:SigmaB = SigmaThreshold = Yreturn Thresholddef GrayThreshold(image, maxval=255):g = GrayHist(image)thresh = GetOSTUThreshold(g)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th, img_new = GrayThreshold(img_gray)th1, img_new_1 = cv2.threshold(img_gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)print(th, th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()

自定义的OTSU方法与opencv中的OTSU方法得到的阈值是一样的，都是125。

6.一维最大熵

最大熵的思想是将图像分为两部分(可以认为是前景和背景)，循环计算两部分的灰度值的熵，使其和最大。取得最大值的灰度即为阈值。
其中熵的概念为：
hist=h(g)∑x=0255h(x)+最小正浮点值hist=\frac{h(g)}{\sum_{x=0}^{255}h(x)}+最小正浮点值hist=∑x=0255​h(x)h(g)​+最小正浮点值

S1=∑g=0Thist(g)S_1=\sum_{g=0}^{T}hist(g)S1​=∑g=0T​hist(g)

E1=−∑g=0Thist(g)S1×log(∑g=0Thist(g)S1)E_1=-\frac{\sum_{g=0}^{T}hist(g)}{S_1}\times log(\frac{\sum_{g=0}^{T}hist(g)}{S_1})E1​=−S1​∑g=0T​hist(g)​×log(S1​∑g=0T​hist(g)​)

E2=−∑g=T+1255hist(g)1−S1×log(∑g=T255hist(g)1−S1)E_2=-\frac{\sum_{g=T+1}^{255}hist(g)}{1-S_1}\times log(\frac{\sum_{g=T}^{255}hist(g)}{1-S_1})E2​=−1−S1​∑g=T+1255​hist(g)​×log(1−S1​∑g=T255​hist(g)​)
计算出最大的E1+E2的灰度值即为阈值。

# coding:utf8import numpy as np
import cv2
from matplotlib import pyplot as plt
import mathdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef Get1DMaxEntropyThreshold(H):V = np.where(H > 0)[0]  # 非0下标I = len(V)if I == 0: return 0if I == 1: return V[0]if I == 2: return V[0] if H[V[0]] < H[V[1]] else V[1]MinValue = V[0]MaxValue = V[-1]Threshold=MinValueHistGramD = np.zeros(256,dtype=np.float64)Amount = np.sum(H[MinValue:MaxValue + 1])eps = np.finfo(np.float64).epsHistGramD[MinValue:MaxValue+1] = H[MinValue:MaxValue+1]/Amount + epsMaxEntropy = 0.0for Y in range(MinValue+1,MaxValue):SumIntegral = 0.for X in range(MinValue,Y+1):SumIntegral += HistGramD[X]EntropyBack = 0.for X in range(MinValue,Y+1):EntropyBack += -HistGramD[X] / SumIntegral * np.log(HistGramD[X] / SumIntegral)EntropyFore = 0.for X in range(Y+1,MaxValue+1):EntropyFore += -HistGramD[X] / (1 - SumIntegral) * np.log(HistGramD[X] / (1 - SumIntegral))if MaxEntropy < EntropyBack + EntropyFore:Threshold = YMaxEntropy = EntropyBack + EntropyForereturn Thresholddef GrayThreshold(image,maxval=255):g = GrayHist(image)thresh = Get1DMaxEntropyThreshold(g)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th, img_new = GrayThreshold(img_gray)th1,img_new_1 = cv2.threshold(img_gray, 0, 255,  cv2.THRESH_BINARY + cv2.THRESH_TRIANGLE)print(th,th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()

计算出的阈值为104。

7.力矩保持法

力矩法的论文可以从此链接中获得力矩法论文。上述论文是英文版的。此方法比较复杂。
实现方法如下：

# coding:utf8import numpy as np
import cv2
from matplotlib import pyplot as plt
import mathdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef AH(H, index):return np.sum(H[:index+1])def BH(H, index):K = np.arange(index+1)return np.sum(H[:index+1]*K)def CH(H, index):K = np.float_power(np.arange(index+1),2)return np.sum(H[:index + 1] * K)def DH(H, index):K = np.float_power(np.arange(index+1),3)return np.sum(H[:index + 1] * K)def GetMomentPreservingThreshold(H):V = np.where(H > 0)[0]  # 非0下标I = len(V)if I == 0: return 0if I == 1: return V[0]if I == 2: return V[0] if H[V[0]] < H[V[1]] else V[1]MaxValue = V[-1]Avec = np.zeros(256, dtype=np.float64)Amount = np.sum(H)for Y in range(256):Avec[Y] = AH(H, Y) / AmountIndex = 0B = BH(H,255)A = AH(H,255)C = CH(H,255)D = DH(H,255)X2 = (B * C - A * D) / (A * C - B * B)X1 = (B * D - C * C) / (A * C - B * B)X0 = 0.5 - (B / A + X2 / 2) / math.sqrt(X2 * X2 - 4 * X1)Min = float(MaxValue)for Y in range(256):if math.fabs(Avec[Y] - X0) < Min:Min = math.fabs(Avec[Y] - X0)Index = Yreturn int(Index)def GrayThreshold(image, maxval=255):g = GrayHist(image)thresh = GetMomentPreservingThreshold(g)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th, img_new = GrayThreshold(img_gray)th1, img_new_1 = cv2.threshold(img_gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_TRIANGLE)print(th, th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()

8.模糊集

具体方法参考基于模糊集理论的一种图像二值化算法的原理、实现效果及代码

# coding:utf8import numpy as np
import cv2
from matplotlib import pyplot as plt
import mathdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef GetHuangFuzzyThreshold(H):Threshold = -1BestEntropy = np.finfo(np.float32).max# 找到第一个和最后一个非0的色阶值V = np.where(H > 0)[0]  # 非0下标I = len(V)if I == 0: return 0if I == 1: return V[0]if I == 2: return V[0] if H[V[0]] < H[V[1]] else V[1]First = V[0]Last = V[-1]# 计算累计直方图以及对应的带权重的累计直方图S = np.zeros(Last+1,dtype=np.int64)W = np.zeros(Last+1,dtype=np.int64)S[0] = H[0]start = First if First>1 else 1for Y in range(start,Last+1):S[Y] = S[Y - 1] + H[Y]W[Y] = W[Y - 1] + Y * H[Y]# 建立公式（4）及（6）所用的查找表Smu = np.zeros(Last+1-First,dtype=np.float32)for Y in range(1,Last+1-First):mu = 1. / (1. + float(Y) / (Last - First))               # 公式（4）Smu[Y] = -mu * math.log(mu) - (1 - mu) * math.log(1 - mu)      # 公式（6）# 迭代计算最佳阈值for Y in range(First,Last+1):Entropy = 0mu = int(np.round(float(W[Y] / S[Y])))          # 公式17for X in range(First,Y+1):Entropy += Smu[abs(X - mu)] * H[X]mu = int(np.round(float((W[Last] - W[Y]) / (S[Last] - S[Y])))) if Y < Last else 0for X in range(Y + 1,Last+1):Entropy += Smu[abs(X - mu)] * H[X]       # 公式8if BestEntropy > Entropy:BestEntropy = Entropy       #取最小熵处为最佳阈值Threshold = Yreturn Thresholddef GrayThreshold(image, maxval=255):g = GrayHist(image)thresh = GetHuangFuzzyThreshold(g)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th, img_new = GrayThreshold(img_gray)th1, img_new_1 = cv2.threshold(img_gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_TRIANGLE)print(th, th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()

9. Kittler最小错误分类法

具体方法见Kittler, J & Illingworth, J (1986), “Minimum error thresholding”, Pattern Recognition 19: 41-47

# coding:utf8import numpy as np
import cv2
from matplotlib import pyplot as plt
import mathdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef GetKittlerMinError(HistGram):value = np.where(HistGram>0)[0] #非0下标index = len(value)if index == 0:return 0elif index == 1:return value[0]elif index == 2:return value[0]else:MinValue = value[0]MaxValue = value[-1]Threshold = -1MinSigma = np.finfo(np.float64).epsfor Y in range(MinValue,MaxValue):PixelBack = 0PixelFore = 0OmegaBack = 0OmegaFore = 0for X in range(MinValue,Y+1):PixelBack += HistGram[X]OmegaBack = OmegaBack + X * HistGram[X]for X in range(Y+1,MaxValue+1):PixelFore += HistGram[X]OmegaFore = OmegaFore + X * HistGram[X]OmegaBack = OmegaBack / PixelBackOmegaFore = OmegaFore / PixelForeSigmaBack = 0SigmaFore = 0for X in range(MinValue,Y+1):SigmaBack = SigmaBack + (X - OmegaBack) * (X - OmegaBack) * HistGram[X]for X in range(Y+1,MaxValue):SigmaFore = SigmaFore + (X - OmegaFore) * (X - OmegaFore) * HistGram[X]if SigmaBack == 0 or SigmaFore == 0:if Threshold == -1:Threshold = Yelse:SigmaBack = math.sqrt(SigmaBack / PixelBack)SigmaFore = math.sqrt(SigmaFore / PixelFore)Sigma = 1 + 2 * (PixelBack * math.log(SigmaBack / PixelBack) + PixelFore * math.log(SigmaFore / PixelFore))if Sigma < MinSigma:MinSigma = SigmaThreshold = Yreturn Thresholddef GrayThreshold(image,maxval=255):g = GrayHist(image)thresh = GetKittlerMinError(g)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th, img_new = GrayThreshold(img_gray)th1,img_new_1 = cv2.threshold(img_gray, 0, 255,  cv2.THRESH_BINARY + cv2.THRESH_TRIANGLE)print(th,th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()

10. ISODATA

参考论文：

Ridler, TW & Calvard, S (1978), “Picture thresholding using an iterative selection method”, IEEE Transactions on Systems, Man and Cybernetics 8: 630-632, http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4310039

# coding:utf8import numpy as np
import cv2
from matplotlib import pyplot as plt
import mathdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef GetIsoDataThreshold(H):V = np.where(H > 0)[0]  # 非0下标I = len(V)if I == 0: return 0if I == 1: return V[0]if I == 2: return V[0] if H[V[0]] < H[V[1]] else V[1]MinValue = V[0]MaxValue = V[-1]threshold = MinValuewhile (True):SumOne = 0SumInteralOne = 0SumTwo = 0SumInteralTwo = 0for i in range(threshold):SumOne += H[i]SumInteralOne += H[i] * ifor i in range(threshold+1,MaxValue):SumTwo += H[i]SumInteralTwo += (H[i] * i)if SumOne > 0 and SumTwo > 0:SumInteralOne /= SumOneSumInteralTwo /= SumTwoif threshold == int(np.round((SumInteralOne + SumInteralTwo) / 2.0)):breakthreshold += 1if threshold > 253:return 0return thresholddef GrayThreshold(image,maxval=255):g = GrayHist(image)thresh = GetIsoDataThreshold(g)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th, img_new = GrayThreshold(img_gray)th1,img_new_1 = cv2.threshold(img_gray, 0, 255,  cv2.THRESH_BINARY + cv2.THRESH_TRIANGLE)print(th,th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()

11. Shanbhag 法

参考论文：

  Shanbhag, Abhijit G. (1994), "Utilization of information measure as a means of image thresholding", Graph. Models Image Process. (Academic Press, Inc.) 56 (5): 414--419, ISSN 1049-9652, DOI 10.1006/cgip.1994.1037

# coding:utf8import numpy as np
import cv2
from matplotlib import pyplot as plt
import mathdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef GetShanbhagThreshold(H):V = np.where(H > 0)[0]  # 非0下标I = len(V)if I == 0: return 0if I == 1: return V[0]if I == 2: return V[0] if H[V[0]] < H[V[1]] else V[1]P1 = np.zeros(256, dtype=np.float64)P2 = np.zeros(256, dtype=np.float64)Amount = np.sum(H)norm_histo = H / AmountP1[0] = norm_histo[0]P2[0] = 1.0 - P1[0]for ih in range(1, 256):P1[ih] = P1[ih - 1] + norm_histo[ih]P2[ih] = 1.0 - P1[ih]first_bin = 0for ih in range(1, 256):if not (math.fabs(P1[ih]) < np.finfo(np.float64).eps):first_bin = ihbreaklast_bin = 255for ih in range(255, first_bin - 1, -1):if not (math.fabs(P2[ih]) < np.finfo(np.float64).eps):last_bin = ihbreakthreshold = -1min_ent = np.finfo(np.float64).maxfor it in range(first_bin, last_bin + 1):ent_back = 0.0term = 0.5 / P1[it]for ih in range(1, it + 1):ent_back -= norm_histo[ih] * math.log(1.0 - term * P1[ih - 1]);ent_back *= terment_obj = 0.0term = 0.5 / P2[it]for ih in range(it + 1, 256):ent_obj -= norm_histo[ih] * math.log(1.0 - term * P2[ih])ent_obj *= termtot_ent = math.fabs(ent_back - ent_obj)if tot_ent < min_ent:min_ent = tot_entthreshold = itreturn thresholddef GrayThreshold(image, maxval=255):g = GrayHist(image)thresh = GetShanbhagThreshold(g)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th, img_new = GrayThreshold(img_gray)th1, img_new_1 = cv2.threshold(img_gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_TRIANGLE)print(th, th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()

12. Yen法

参考论文：

 1) Yen J.C., Chang F.J., and Chang S. (1995) "A New Criterion  for Automatic Multilevel Thresholding" IEEE Trans. on Image  Processing, 4(3): 370-3782) Sezgin M. and Sankur B. (2004) "Survey over Image Thresholding Techniques and Quantitative Performance Evaluation" Journal of  Electronic Imaging, 13(1): 146-165

# coding:utf8import numpy as np
import cv2
from matplotlib import pyplot as plt
import mathdef GrayHist(img):grayHist = np.zeros(256, dtype=np.uint64)for v in range(256):grayHist[v] = np.sum(img == v)return grayHistdef GetYenThreshold(H):V = np.where(H > 0)[0]  # 非0下标I = len(V)if I == 0: return 0if I == 1: return V[0]if I == 2: return V[0] if H[V[0]] < H[V[1]] else V[1]P1 = np.zeros(256, dtype=np.float64)P1_sq = np.zeros(256, dtype=np.float64)P2_sq = np.zeros(256, dtype=np.float64)norm_histo = H/ np.sum(H)P1[0] = norm_histo[0]for ih in range(1, 256):P1[ih] = P1[ih - 1] + norm_histo[ih]P1_sq[0] = norm_histo[0] * norm_histo[0]for ih in range(1, 256):P1_sq[ih] = P1_sq[ih - 1] + norm_histo[ih] * norm_histo[ih]P2_sq[255] = 0.0for ih in range(254, -1, -1):P2_sq[ih] = P2_sq[ih + 1] + norm_histo[ih + 1] * norm_histo[ih + 1]threshold = -1max_crit = np.finfo(np.float64).epsfor it in range(256):l1 = math.log(P1_sq[it] * P2_sq[it]) if P1_sq[it] * P2_sq[it] > 0.0 else 0.0l2 = math.log(P1[it] * (1.0 - P1[it])) if P1[it] * (1.0 - P1[it]) > 0.0 else 0.0crit = -1.0 * l1 + 2 * l2if crit > max_crit:max_crit = critthreshold = itreturn thresholddef GrayThreshold(image, maxval=255):g = GrayHist(image)thresh = GetYenThreshold(g)threshImage_out = image.copy()# 大于阈值的都设置为maxvalthreshImage_out[threshImage_out > thresh] = maxval# 小于阈值的都设置为0threshImage_out[threshImage_out <= thresh] = 0return thresh, threshImage_outif __name__ == "__main__":img = cv2.imread('bird.png')img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)th, img_new = GrayThreshold(img_gray)th1, img_new_1 = cv2.threshold(img_gray, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_TRIANGLE)print(th, th1)plt.subplot(131), plt.imshow(img_gray, cmap='gray')plt.title('Original Image'), plt.xticks([]), plt.yticks([])plt.subplot(132), plt.imshow(img_new, cmap='gray')plt.title('Image'), plt.xticks([]), plt.yticks([])plt.subplot(133), plt.imshow(img_new_1, cmap='gray')plt.title('CV2 Image1'), plt.xticks([]), plt.yticks([])plt.show()

十三种基于直方图的图像全局二值化算法原理、实现、代码及效果相关推荐

  1. 十三种基于直方图的图像全局二值化算法原理、实现、代码及效果（转）

    源:十三种基于直方图的图像全局二值化算法原理.实现.代码及效果.

  2. 十三种基于直方图的图像全局二值化算法原理、实现、代码及效果。

    图像二值化的目的是最大限度的将图象中感兴趣的部分保留下来,在很多情况下,也是进行图像分析.特征提取与模式识别之前的必要的图像预处理过程.这个看似简单的问题,在过去的四十年里受到国内外学者的广泛关注,产 ...

  3. 基于直方图的图像全局二值化算法原理、实现--基于谷底最小值的阈值

    1.描述: 此方法实用于具有明显双峰直方图的图像,其寻找双峰的谷底作为阈值,但是该方法不一定能获得阈值,对于那些具有平坦的直方图或单峰图像,该方法不合适. 2.实现过程: 该函数的实现是一个迭代的过程 ...

  4. 场景文本检测算法 可微分二值化DBNet原理与代码解析

    目录 原理介绍 Label Generation Loss函数 后处理 论文 https://arxiv.org/abs/1911.08947 代码 https://github.com/WenmuZ ...

  5. matlab最小分类错误全局二值化算法

    转自:http://download.csdn.net/detail/hupeng810/1511870 function imagBW = kittlerMet(imag) % KITTLERMET ...

  6. matlab确定灰度阈值T,matlab灰度图像二值化【灰度图像二值化算法研究】

    摘要: 在很多图像处理的过程中,经常需要对灰度图像进行二值化.本文对几种常用的图像二值化算法进行了阐述,并通过仿真,进行比较研究.根据实验结果,阐明了各种算法的优缺点. Abstract: The b ...

  7. 二值化算法OTSU源码解析

    点击上方"小白学视觉",选择加"星标"或"置顶" 重磅干货,第一时间送达 概述: 本文中小编将会跟大家分享一下OpenCV3.1.0中图像二 ...

  8. 【图像处理】——图像的二值化操作及阈值化操作（固定阈值法（全局阈值法——大津法OTSU和三角法TRIANGLE）和自适应阈值法（局部阈值法——均值和高斯法））

    目录 一.二值化的概念(实际上就是一个阈值化操作) 1.概念: 2.实现方法 3.常用方法 二.阈值类型 1.常见阈值类型(主要有五种类型) (1)公式描述 (2)图表描述 2.两种特殊的阈值算法(O ...

  9. python-opencv 图像阈值二值化

    本文讲解基于OpenCV-python的图像二值化API及浅显原理讲解 文章目录 一. 阈值 1. 简单阈值 2. 自适应阈值 二. 图像二值化 1. 全局图像二值化 2. 局部图像二值化 3. Ot ...

最新文章

  1. EEMD算法的基本原理
  2. Centos中Redis的下载编译与安装(超详细)
  3. python类百度百科_Python抓取百度百科数据
  4. C语言学习之插入排序
  5. linux x64下安装oracle 11g
  6. [原创]完美开启Win8中管理员Administrator帐户
  7. （41）zabbix监控api接口性能及可用性 天气预报api为例
  8. html文件form根目录,HTML ,form 和 link 使用根目录 的问题,我已经上图了!
  9. mysqli 语句和mysql语句一样吗_mysqli语句的用法
  10. 设计开发-API代付系统/自动代付系统
  11. 标准差、均方误差、均方根误差、平均绝对误差
  12. Word中批量进行中英文标点的转换
  13. B题Tomb Raider ---- 一 。启动emule客户端
  14. win10计算机入门使用教程,win10系统使用教程_windows10基本使用教程图文步骤
  15. ipq wifi校准
  16. 看天下网络资讯浏览器 下载
  17. Redis原子计数器incr，防止并发请求
  18. NOIP2016·洛谷·天天爱跑步
  19. python中if条件语句的代码实例
  20. 高校教材管理系统mysql_高校教材管理系统的设计与实现

热门文章

  1. 了解NVIDIA显卡驱动（包括CUDA、CUDA Driver、CUDA Toolkit、CUDNN、NCVV）
  2. 数据结构--01|逻辑结构和物理结构（存储结构）
  3. 推荐一些好用的视频处理工具
  4. 【应急基础】安全应急响应工具年末大放送（含下载）
  5. 小明有5本新书，要借给A、B、C三位小朋友，若每人每次只能借一本，  共有多少种借书的方案？
  6. 《精彩绝伦的CSS》——提示（二）无单位的行高值
  7. 网上打印室打印收费标准是什么？
  8. 用友成立用友畅捷通，专针对中小企业管理软件公司
  9. 【MySQL】MySQL 架构
  10. 笔记本键盘上输入0变成/怎么回事？