[Python知识库]大律法OTSU【详细】（python实现）

一、 前言

图像的灰度直方图

灰度直方图是关于灰度级分布的函数，是对图像中灰度级分布的统计。灰度直方图是将数字图像中的所有像素，按照灰度值的大小，统计其出现的频率。灰度直方图是灰度级的函数，它表示图像中具有某种灰度级的像素的个数，反映了图像中某种灰度出现的频率。【from百度百科】

计算直方图

1.使用opencv的函数

cv2.calcHist(images, channels, mask, histSize, ranges)

• 参数1：要计算的原图，以方括号的传入，如：[img]。
• 参数2：类似前面提到的dims，灰度图写[0]就行，彩色图B/G/R分别传入[0]/[1]/[2]。
• 参数3：要计算的区域ROI，计算整幅图的话，写None。
• 参数4：也叫bins,子区段数目，如果我们统计0-255每个像素值，bins=256；如果划分区间，比如0-15,16-31…240-255这样16个区间，bins=16。
• 参数5：range,要计算的像素值范围，一般为[0,256)。

eg:

import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt
img = cv.imread('test.jpg')
# hist是256x1数组，每个值对应于该图像中具有相应像素值的像素数
hist = cv.calcHist([img],[0],None,[256],[0,256])
# 绘制直方图
plt.plot(hist)
plt.show()

2. 使用Numpy函数

bincount(x, weights=None, minlength=None)

相对应的函数详解可参考：numpy.bincount详解
eg：

import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt
img = cv.imread('test.jpg')
# ravel()函数将二维矩阵展平变成一维数组
hist = np.bincount(img.ravel(), minlength=256) 
# 绘制直方图
plt.plot(hist)
plt.show()

二、 OTSU算法简介

OTSU算法是一种图像二值化算法。其算法假设图像中存在一个阈值 T ，判断图像中每个像素与T的大小关系，可以将所有像素分为背景 C₀ 和 前景 C₁ 两类，且当选取到最佳T阈值时，背景部分与前景部分的差别最大，而OTSU 算法利用最大类间方差来衡量这一差别。由于OTSU算法采用了最大类间方法的思想，因此该算法也被OTSU最大类间方差法

缺点：对图像噪声敏感；只能针对单一目标分割；当目标和背景大小比例悬殊、类间方差函数可能呈现双峰或者多峰，这个时候效果不好

三、 数学原理

该算法的原理主要涉及 均值、方差等概念和部分公式的推导

1. 将图片的像素值分为 [1, 2, …, L] 个水平，用n_i 表示各个水平像素值的像素个数，那么很容易得到总像素个数N为:N = n₁ + n₂ + … + n_L

2. 利用像素值对应个数与总数的商作为某个像素值出现的概率，定义p_i
   $$p_i=n_i/N,p_i\ge 0,\sum _{i=1}^L{p}_i=1\text{\hspace{1em}(1)}$$

3. 定义两个量w₀ , w₁为C₀,C₁的局部概率之和，并且得到二者的关系
   $$w_0=Pr(C_0)=\sum _{i=1}^k{p}_i=w(k)\text{\hspace{1em}(2)}$$
   $$w_1=Pr(C_1)=\sum _{i=k+1}^L{p}_i=1?w(k)\text{\hspace{1em}(3)}$$

4. 由此我们得到总的数学期望【期望是平均数随样本趋于无穷的极限，所以可以将期望 ≈ 均值】和C₀ , C₁各自的数学期望并指出三者的关系，式中 i 代表像素值，除于各自概率的和用于进行归一。
   $$u_0=\sum _{i=1}^{k}Pr(i|C_0)=\sum _{i=1}^{k}i?p_i/w_0=\frac{u(k)}{w(k)}\text{\hspace{1em}(4)}$$
   $$u_1=\sum _{i=k+1}^{L}Pr(i|C_1)=\sum _{i=k+1}^{L}i?p_i/w_1=\frac{u(T)?u(k)}{1?w(k)}\text{\hspace{1em}(5)}$$
   $$w(k)=\sum _{i=1}^k{p}_i\text{\hspace{1em}(6)}$$
   $$u(k)=\sum _{i=1}^{k}i?p_i\text{\hspace{1em}(7)}$$
   $$u_T=u(L)=\sum _{i=1}^{L}i?p_i\text{\hspace{1em}(8)}$$
   上述公式间的关系为：
   $$w_0{u}_0+w_1{u}_1=u_T,w_0+w_1=1\text{\hspace{1em}(9)}$$

5. 各自局部方差和总方差为
   $${\sigma }_0^2=\sum _{i=1}^k(i?u_0{)}^{2}Pr(i|C_0)=\sum _{i=1}^k(i?u_0{)}^2{p}_i/w_0\text{\hspace{1em}(10)}$$
   $${\sigma }_1^2=\sum _{i=k+1}^L(i?u_1{)}^{2}Pr(i|C_1)=\sum _{i=k+1}^L(i?u_1{)}^2{p}_i/w_1\text{\hspace{1em}(11)}$$
   $${\sigma }_T^2=\sum _{i=1}^L(i?u_T{)}^2{p}_i\text{\hspace{1em}(12)}$$

6. within_class variance 类内差异
   $${\sigma }_w^2=w_0{\sigma }_0^2+w_1{\sigma }_1^2\text{\hspace{1em}(13)}$$
   between_class variance 类间差异
   $${\sigma }_b^2=w_0(u_0?u_T{)}^2+w_1(u_1?u_T{)}^2\text{\hspace{1em}(14)}$$

7. 当最大化类间差距时，即最小化类内差距，因为 二者的和为定值。
   $${\sigma }_w^2+{\sigma }_b^2={\sigma }_T^2\text{\hspace{1em}(15)}$$

四、 python实现

假设一图像为灰度图，大小为M*N，背景部分较暗，OTSU算法步骤如下：

1. 选取一个分割阈值T，统计计算属于前景的像素点数（像素值小于T）在全部像素点的占比 $w_0$ 以及其平均灰度 $u_0$ ; 同理统计属于背景部分的 $w_1$ 和 $u_1$ ;
2. 计算图像全部像素的平均灰度 $u_T=w_0?u_0+w_1?u_1$ 以及类间方差

$$g=w_0?(u_0?u_T{)}^2+w_1?(u_1?u_T{)}^2=w_0?w_1?(u_0?u_1{)}^2\text{\hspace{1em}(16)}$$

3. 遍历全部分割阈值T，重复1-2步骤，得到使类间方差 g 最大的阈值T；

import numpy as np
import cv2 as cv

#将图片转为灰度图
# https://blog.csdn.net/weixin_43635647/article/details/99625105
img = cv.imread('test.jpg', 0)
cv.imshow("img", img)
cv.waitKey()

def OTSU(img_gray, GrayScale):
    assert img_gray.ndim == 2, "must input a gary_img"  #shape有几个数字, ndim就是多少
    img_gray = np.array(img_gray).ravel().astype(np.uint8)
    u1=0.0#背景像素的平均灰度值
    u2=0.0#前景像素的平均灰度值
    th=0.0

    #总的像素数目
    PixSum=img_gray.size
    #各个灰度值的像素数目
    PixCount=np.zeros(GrayScale)
    #各灰度值所占总像素数的比例
    PixRate=np.zeros(GrayScale)
    #统计各个灰度值的像素个数
    for i in range(PixSum):
        #默认灰度图像的像素值范围为GrayScale
        Pixvalue=img_gray[i]
        PixCount[Pixvalue]=PixCount[Pixvalue]+1
    
    #确定各个灰度值对应的像素点的个数在所有的像素点中的比例。
    for j in range(GrayScale):
        PixRate[j]=PixCount[j]*1.0/PixSum
    Max_var = 0
    #确定最大类间方差对应的阈值
    for i in range(1,GrayScale):#从1开始是为了避免w1为0.
        u1_tem=0.0
        u2_tem=0.0
        #背景像素的比列
        w1=np.sum(PixRate[:i])
        #前景像素的比例
        w2=1.0-w1
        if w1==0 or w2==0:
            pass
        else:#背景像素的平均灰度值
            for m in range(i):
                u1_tem=u1_tem+PixRate[m]*m
            u1 = u1_tem * 1.0 / w1
             #前景像素的平均灰度值
            for n in range(i,GrayScale):
                u2_tem = u2_tem + PixRate[n]*n
            u2 = u2_tem / w2
            #print(u1)
            #类间方差公式：G=w1*w2*(u1-u2)**2
            tem_var=w1*w2*np.power((u1-u2),2)
            #print(tem_var)
            #判断当前类间方差是否为最大值。
            if Max_var<tem_var:
                Max_var=tem_var#深拷贝，Max_var与tem_var占用不同的内存空间。
                th=i
    return th 
th = OTSU(img,256)
print("使用numpy的方法：" + str(th))		# 结果为 135

直接使用 opencv的函数

import numpy as np
import cv2 as cv
#该函数返回的第一个值就是输入的thresh值，第二个就是处理后的图像
img = cv.imread('test.jpg',0)
retVal, a_img = cv.threshold(img, 0, 255, cv.THRESH_OTSU)
print("使用opencv函数的方法：" + str(retVal))	# 结果为 134
cv.imshow("a_img",a_img)
cv.waitKey()

五、 参考文章

• 图像处理基础(8)：图像的灰度直方图、直方图均衡化、直方图规定化（匹配）
• OpenCV图像处理-直方图
• A Threshold Selection Method from Gray-Level Histograms
• 大津法（OTSU 最大类间方差法）详细数学推导
• OTSU（大津律法）数学原理及源码解析
• 大津法Python实现