基于OTSU的图像多阈值分割算法

目录

• 基于OTSU的图像多阈值分割算法
  • • 实验需求分析：
    • 实验基础知识：
      • 什么是Otsu阈值分割？
        • 类内方差（the within-class variance）：
        • 类间方差（the between-class variance）：
        • 灰度级总方差（the total variance of levels）
        • 小总结
    • 实验原理分析:
      • Otsu多阈值分割原理
      • 实际问题分析及python代码实现
        • 灰度级直方图
        • 二值化
        • 从理论到实践
          • 原图像
          • 处理后的分层图像
          • 实现最终代码
    • 感受
    • 后话

实验需求分析：

• 要求在以灰度图读入摄像头画面的同时对画面中的太阳进行分割，以获取我们所需的红外灯所成的正确图像；
• 了解otsu之后我们知道只需在其单阈值分类基础上加多一个阈值即可完成本实验；
• 用otsu分割时我们知道起码需要两个阈值k1、k2，分成三类：
  • 前景
    • 太阳
  • 中间区域
    • 红外灯
  • 背景
    • 除此上面两者之外的所有噪点及环境

实验基础知识：

什么是Otsu阈值分割？

1979年NOBUYUKI OTSU在他的文章"A Threshold Selection Method from Gray-Level Histograms"提出了A Threshold Selection Method from Gray-Level Histograms.

看了网上较多大佬写的博客，但为了避免理解偏差，还是去找原论文了解了一下，也觉得有必要自己写写加深理解（参考原论文）；

在此之前我们需要一些基本知识，以便后面的知识读起来变得更加易于接受，又由于时间关系，我们直接在用到的地方给出相关知识：

• 灰度级：一说为用来区分图像中不同灰度的最大数量，但本笔记中的涉及到的灰度级只是一系列级别，比如一幅图像有L个灰度级[1,2,…,L]，（正如这篇论文而言，我们也应该输出图片的Gray-Level Histograms看看）（这里的灰度级直方图是指[0,255]中256个灰度级为横坐标，对应像素数为纵坐标的一幅直方图）

我们将一幅图片分成L层灰度级[1,2,...,L],且我们规定位于第i个灰度级中的像素点数为n_i,则有总像素数N

\[N=\sum^{L}_{i=1}n_i \]

我们将这个灰度级直方图当作一个概率分布，记：

\[p_i=\frac{n_i}{N},\quad p_i\ge0,\sum^{L}_{i=1}p_i=1 \]

（原论文做的是二值化）现在我们打算把众多的像素点用一个阈值k分成两类C_0和C_1，而且C_0类像素点的灰度级都处在[1,...,k]，C_1像素点的灰度级都处在[k+1,...,L]；

类似地，我们用概率分布的思想，计算两类出现的概率

\[\omega_0=Pr(C_0)=\sum^{k}_{i=1}p_i\\ 记\omega(k)=\omega_0\\ \omega_1=Pr(C_1)=\sum^{L}_{i=k+1}p_i\\ 则\omega_1=1-\omega(k) \]

为了更好的理解下一步的过程，我们可以继续学习额外知识， 参考这篇文章零阶矩、一阶矩、二阶矩

则两类中的平均灰度级（一阶矩、数学期望）分别是，（此时我们分出两类中灰度级期望，我们得知道我们的目标是分类，所以到目前为止还是在寻找最能分开两类的指标，或者说是特征值，其中平均灰度值可能是其中一个指标）

\[\mu_0=\sum^{k}_{i=1}i\;Pr(i\mid C_0)=\sum^{k}_{i=1}\frac{ip_i}{\omega_0}=\frac{\mu(k)}{\omega(k)} \quad(1) \\ \mu_1=\sum^{L}_{i=k+1}i\;Pr(i\mid C_1)=\sum^{L}_{i=k+1}\frac{ip_i}{\omega_1}=\frac{\mu_T-\mu(k)}{1-\omega(k)} \quad(2)\\ (上式需要注意的是\mu_T指的是\mu_{Total},\mu(k)、\omega(k)会在下面的式子解释)\\ (一个小细节:我们需要区分\mu_0和\mu(k),前者是类内的条件概率，而后者是针对全局而言) \]

其中

\[\omega(k)=\sum^{k}_{i=1}p_i\\ \mu(k)=\sum^{k}_{i=1}ip_i\quad(第一层灰度级平均值)\\ (表述为与阈值有关的样式更有利于我们下面讨论)\\ 则\mu_T=\mu(L)=\sum^{L}_{i=1}ip_i\quad(原图片灰度级平均值) \]

由(1)(2)式我们容易得到，也容易知道

\[\omega_0\mu_0+\omega_1\mu_1=\mu_T,\quad\omega_0+\omega_1=1 \]

两类的方差分别是（类内方差、二阶矩）

\[\sigma^2_0=\sum^{k}_{i=1}(i-\mu_0)^{2}Pr(i\mid C_0)=\sum^{k}_{i=1}(i-\mu_0)^{2}\frac{p_i}{\omega_0}\\ \sigma^2_1=\sum^{L}_{i=k+1}(i-\mu_1)^{2}Pr(i\mid C_1)=\sum^{L}_{i=k+1}(i-\mu_1)^{2}\frac{p_i}{\omega_1}\\ \]

下面的做法是为了评估灰度级k时的分类好坏程度，这里论文作者用了一种判别分析方法,为此引入三个因子（类似于归一化矩但又不太像）

\[\lambda=\frac{\sigma_B^2}{\sigma_W^2}\\ \kappa =\frac{\sigma_T^2}{\sigma_W^2}\\ \eta =\frac{\sigma_B^2}{\sigma_T^2}\\ \]

其中

类内方差（the within-class variance）：

\[\sigma_W^2=\omega_0\sigma_0^2+\omega_1\sigma_1^2 \quad (方差与类出现概率积的和)\\ \quad \]

类间方差（the between-class variance）：

\[\begin{align} \sigma_B^2 & = \omega_0(\mu_0-\mu_T)^2+\omega_1(\mu_1-\mu_T)^2\\ & = \omega_0\mu_0^{2}+ \omega_1\mu_1^{2} +(\omega_0+\omega_1)\mu_T^{2} -2(\omega_0\mu_0+\omega_1\mu_1)\mu_T\\ &=\omega_0\mu_0^{2}+ \omega_1\mu_1^{2} +\mu_T^{2} -2\mu_T^2\\ &=\omega_0\mu_0^{2}+ \omega_1\mu_1^{2}-\mu_T^2\\ &=\omega_0\mu_0^{2}+ \omega_1\mu_1^{2}-(\omega_0\mu_0+\omega_1\mu_1)^2\\ &=\omega_0\omega_1\mu_0^{2}+\omega_0\omega_1\mu_1^{2}-2\omega_0\omega_1\mu_0\mu_1\\ &=\omega_0\omega_1(\mu_0-\mu_1)^2 \end{align} \]

灰度级总方差（the total variance of levels）

\[\sigma_T^2=\sum_{i=1}^{L}(i-\mu_T)^2p_i \]

则目前我们的问题回到了最大化刚刚定义的一个判别因子即可：（我们想达到最好的分类效果，即类内方差最小同时还想类间方差最大，所以用比值因子来定义最合适）

\[我们知道可知类内方差+类间方差=总方差：\\ \sigma_W^2+\sigma_B^2=\sigma_T^2 \]

则有

\[\kappa =\frac{\sigma_T^2}{\sigma_W^2}=\frac{\sigma_W^2+\sigma_B^2}{\sigma_W^2}=1+\frac{\sigma_B^2}{\sigma_W^2}=1+\lambda\\ \eta=\frac{\sigma_T^2-\sigma_W^2}{\sigma_T^2} =1-\frac{1}{\kappa}=1-\frac{1}{1+\lambda}=\frac{\lambda}{1+\lambda}\\ (至此我们知道每个因子都是彼此正相关的，所以上面所说的只要最大化其中一个因子即可) \]

又前面的类内方差和类间方差可知都是关于阈值k的函数，但显然总方差是独立于k的，而且同样需要注意的是类内方差是基于二阶统计数据(类的方差)，而类间方差是基于一阶统计数据(类的平均值)，所以类间方差与总方差之比是应该最简单可以代表k的好坏程度的，所以我们只需找到一个k来最大化类间方差与总方差之比即可，又之前我们刚好定义了这个因子

\[\eta(k)=\sigma_B^2(k)/\sigma_T^2\quad*\\ \sigma_B^2(k)=\omega_0\omega_1(\mu_0-\mu_1)^2\\(我们显然希望上式表示为k的函数形式) \\我们上面定义过 \omega_0=\omega(k)\quad\omega_1=1-\omega(k)\\ \mu_0=\frac{\mu(k)}{\omega(k)}\quad\mu_1=\frac{\mu_T-\mu(k)}{1-\omega(k)} \\代入得 \sigma_B^2(k)=\frac{(\mu_T\omega(k)-\mu(k))^2}{\omega(k)(1-\omega(k))} \]

对于*式而言，我们知道因为总方差不受k影响是个定值，我们只需找到一个k值令类间方差最大即可（或者令类内方差最小即可），而找k值的方法便足够丰富多样了（如何降低复杂度等等）

小总结

提取上面重要的步骤其实我们只需几步即可

\[两类出现的概率\\ \omega_0=\omega(k),\;\omega_1=1-\omega(k)\\ 两类分别的期望\\ \mu_0=\frac{\mu(k)}{\omega(k)},\;\mu_1=\frac{\mu_T-\mu(k)}{1-\omega(k)}\\ 令\sigma^{2}_B=\omega_0(\mu_0-\mu_T)^2+\omega_1(\mu_1-\mu_T)^2取得最大值即可\\ \]

或者更简单(写成这样易于我们去拓展)

\[令\omega(k)=\omega_0 \\\omega_0\mu_0+\omega_1\mu_1=\mu_T,\quad\omega_0+\omega_1=1\\ 令\sigma^{2}_B=\omega_0(\mu_0-\mu_T)^2+\omega_1(\mu_1-\mu_T)^2取得最大值即可\\ \]

不要忘记我们其中一些重要参数

\[\sigma_B^2(k)=\frac{(\mu_T\omega(k)-\mu(k))^2}{\omega(k)(1-\omega(k))}\\ \omega(k)=\sum^{k}_{i=1}p_i\;(w_0=\omega(k))\\ \mu(k)=\sum^{k}_{i=1}ip_i\quad(第一类灰度级平均值,\mu(k)=w_0\mu_0)\\ 则\mu_T=\mu(L)=\sum^{L}_{i=1}ip_i\quad(原图片灰度级平均值) \]

实验原理分析:

Otsu多阈值分割原理

\[令\omega(k_1)=\omega_0,\omega(k_2)=\omega_1 \\\omega_0\mu_0+\omega_1\mu_1+\omega_2\mu_2=\mu_T,\quad\omega_0+\omega_1+\omega_2=1\\ 令\sigma^{2}_B=\omega_0(\mu_0-\mu_T)^2+\omega_1(\mu_1-\mu_T)^2+\omega_2(\mu_2-\mu_T)^2取得最大值即可\\ \]

\[其中\\ \omega_0=\omega(k_1)=\sum^{k_1}_{i=1}p_i\\ \omega_1=\omega(k_2)=\sum^{k_2}_{i=k_1+1}p_i\\ (\omega_2容易知道)\\ \mu(k_1)=\sum^{k_1}_{i=1}ip_i\quad(第一类灰度级平均值)\\ \mu(k_2)=\sum^{k_2}_{i=k_1+1}ip_i\quad(第二类灰度级平均值)\\ 则\mu_T=\mu(L)=\sum^{L}_{i=1}ip_i\quad(原图片灰度级平均值) \]

\[\begin{align} \sigma^{2}_B & = \omega_0(\mu_0-\mu_T)^2+\omega_1(\mu_1-\mu_T)^2+\omega_2(\mu_2-\mu_T)^2 \\&=\omega(k_1)(\frac{\mu(k_1)}{\omega(k_1)} -\mu_T)^2 +\omega(k_2)(\frac{\mu(k_2)}{\omega(k_2)} -\mu_T)^2 +(1-\omega(k_1)-\omega(k_2))(\frac{1-\mu(k_1)-\mu(k_2)}{(1-\omega(k_1)-\omega(k_2))} -\mu_T)^2 \\&=\frac{(\mu(k_1) -\omega(k_1)\mu_T)^2}{\omega(k_1)} +\frac{(\mu(k_2) -\omega(k_2)\mu_T)^2}{\omega(k_2)} +\frac{(1-\mu(k_1)-\mu(k_2) -(1-\omega(k_1)-\omega(k_2))\mu_T)^2}{1-\omega(k_1)-\omega(k_2)} \end{align} \]

最终就是一个二元函数（但不连续）求极值问题，

为了统一对称且公式好看，笔者这样拓展了一下上述公式

\[\sigma^{2}_B=\sum^{layers-1}_{i=0}\frac{(U_i-W_iU_T)^2}{W_i}\\ where\\ U_i是第i层的灰度级平均值\\ U_0=\sum^{k_1}_{j=1}jp_j\quad W_0=\sum^{k_1}_{j=1}p_j\quad(第0层:背景)\\ 我们不妨令k_0=0则有: U_0=\sum^{k_1}_{k_0+1}jp_j\quad W_0=\sum^{k_1}_{k_0+1}p_j\\ U_i=\sum^{k_{i+1}}_{k_i+1}jp_j\quad W_i=\sum^{k_{i+1}}_{k_i+1}p_j\quad(第i层)\\ 则到最后时有\\ U_{layers-1}=\sum^{k_{layers}}_{k_{layers-1}+1}jp_j\quad W_{layers-1}=\sum^{k_{layers}}_{k_{layers-1}+1}p_j\quad(第layers-1层)\\ 另外U_T=\mu_T=\sum^{L}_{j=1}jp_j=\sum^{k_{layers}}_{k_0+1}jp_j,\\ (细节:k_0=0，k_{layers}=L，真正的阈值是[k_1,\cdots,k_{layers-1}]共layers-1个) \]

实际问题分析及python代码实现

灰度级直方图

我们需要在[0,255]这256个灰度级中找出其中两个值k1,k2来分离太阳，红外灯，背景共三层：其中太阳应该是灰度级较高且在直方图中比较集中的发布在最右端且像素点占比不多，而我们的红外灯次之，应该位于中间区段，我们用python画灰度级直方图：

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


def showimg_pro_eachgray(imagepath):  # 显示每一个灰度级的像素数概率
    img = cv.imread(imagepath, 0)  # 以灰度图形式读取图片
    pixels = cv.calcHist([img], [0], None, [256], [0, 256])  # 计算每个灰度级中所含像素数，返回的是一个（256,1）的数组
    p = pixels / (img.shape[0] * img.shape[1])  # 获得每个灰度级中像素数占总像素数比例则我们获得了pi的一个向量
    x = np.linspace(0, 255, 256)  # 横坐标灰度级别
    plt.subplot(2,1,1)
    plt.imshow(img,'gray') # 如果不加gray就会出现以前hsv与rgb的问题
    plt.title('origin image')
    plt.colorbar(orientation="horizontal") #水平放置颜色条
    plt.subplot(2,1,2)
    plt.bar(x, p.ravel(), 0.9, alpha=1, color='b')
    plt.title('Histogram of the probablity for each gray level',y=-0.4) #调节标题上下移动的方法
    plt.show()


if __name__ == '__main__':
    # 显示直方图
    showimg_pro_eachgray("sun.jpg")

二值化

我们尝试对这张自然界太阳图片进行处理：用127的灰度级为阈值进行二值化；直接用otsu确定阈值进行二值化；两者结合的方法进行二值化，最后输出三种类型图片及每种方法的阈值：

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt
img = cv.imread('sun.jpg',0)
# 全局阈值
ret1,th1 = cv.threshold(img,127,255,cv.THRESH_BINARY) # 第二个参数是我们自己手动设置
# Otsu阈值
ret2,th2 = cv.threshold(img,0,255,cv.THRESH_OTSU)
# 混合方法
ret3,th3 = cv.threshold(img,0,255,cv.THRESH_BINARY+cv.THRESH_OTSU)
plt.subplot(2,2,1)
plt.imshow(img,'gray')
plt.title('origin image')
plt.subplot(2,2,2)
plt.imshow(th1,'gray')
plt.title('BINARY:'+str(ret1))
plt.subplot(2,2,3)
plt.imshow(th2,'gray')
plt.title('OTSU:'+str(ret2),y=-0.4)
plt.subplot(2,2,4)
plt.imshow(th2,'gray')
plt.title('BINARY+OTSU:'+str(ret3),y=-0.4)
plt.show()

但是显然单纯的二值化不是我们的目的，因为我们会有红外灯的光，还要进行进一步的分割:

从理论到实践

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


def process_pro_eachgray(imagepath):  # 显示每一个灰度级的像素数概率
    img = cv.imread(imagepath, 0)  # 以灰度图形式读取图片
    pixels = cv.calcHist([img], [0], None, [256], [0, 256])  # 计算每个灰度级中所含像素数，返回的是一个（256,1）的数组
    p = pixels / (img.shape[0] * img.shape[1])  # 获得每个灰度级中像素数占总像素数比例则我们获得了p_i的一个向量
    x=np.linspace(1,256,256) #灰度级我们定义为从0到255好像不合适这样做平均灰度级时会忽略第一个数据，所以我们定义为从1到256
    x=x.reshape(256,1)  # 这一步是因为我们需要的不是秩为1的向量，可以不加这一句输出看看shape的区别，但这样我们得到了（256,1）的灰度级向量
    maxvar=0
    th1=0
    th2=0
    cnt=0
    for k1 in range(1,257):
        for k2 in range(k1+1,257):
            w0 = np.sum(p[0:k1])
            w1 = np.sum(p[k1:k2])
            w2 = np.sum(p[k2:256])  # 注意python是左闭右开np.sum(p[0:256])=1
            u0 = np.dot(x[0:k1].T, p[0:k1])
            u1 = np.dot(x[k1:k2].T, p[k1:k2])
            u2 = np.dot(x[k2:256].T, p[k2:256])
            W = np.array([[w0, w1, w2]]).T  # 注意是两层中括号，用np.array组成的向量可以转置,输出矩阵大小等等操作
            U = np.array([[int(u0), int(u1), int(u2)]]).T
            UT = np.sum(U)
            temp = (U - UT * W) * (U - UT * W) / W
            varbetween = np.sum(temp)
            if varbetween>maxvar:
                cnt=cnt+1
                maxvar=varbetween
                th1=k1
                th2=k2
                print("第"+str(cnt)+"次找到最佳阈值："+str(th1)+","+str(th2)+"此时类间方差为"+str(maxvar))
    return th1,th2      

if __name__ == '__main__':
    process_pro_eachgray("sun.jpg")

本代码是用了最简单的暴力枚举方法（复杂度为O(n^2)以后有待改进），通过一个简单的450*600像素的图片找出了最大类间方差值，同时返回两个阈值，最终结果如下：

第120次找到最佳阈值：76,140此时类间方差为2430.652302720893

这里我们便通过了两个阈值将图片灰度值分成了三层，再稍微改进代码将不同层区分并显示出来，三层：白色，灰色，黑色，最终结果：

原图像

处理后的分层图像

实现最终代码

import cv2.cv2 as cv
import numpy as np
np.seterr(invalid='ignore') #之所以加这一句是因为下面除法的时候会除零报错
import matplotlib.pyplot as plt


def otsu3layers(imagepath):  # 显示每一个灰度级的像素数概率
    img = cv.imread(imagepath, 0)  # 以灰度图形式读取图片

    high, wide = img.shape
    pixels = cv.calcHist([img], [0], None, [256], [0, 256])  # 计算每个灰度级中所含像素数，返回的是一个（256,1）的数组
    p = pixels / (high * wide)  # 获得每个灰度级中像素数占总像素数比例则我们获得了p_i的一个向量
    x=np.linspace(1,256,256) #灰度级我们定义为从0到255好像不合适这样做平均灰度级时会忽略第一个数据，所以我们定义为从1到256
    x=x.reshape(256,1)  # 这一步是因为我们需要的不是秩为1的向量，可以不加这一句输出看看shape的区别，但这样我们得到了（256,1）的灰度级向量

    maxvar=0
    th1=0
    th2=0
    cnt=0
    for k1 in range(1,257):
        for k2 in range(k1+1,257):
            w0 = np.sum(p[0:k1])
            w1 = np.sum(p[k1:k2])
            w2 = np.sum(p[k2:256])  # 注意python是左闭右开np.sum(p[0:256])=1
            u0 = np.dot(x[0:k1].T, p[0:k1])
            u1 = np.dot(x[k1:k2].T, p[k1:k2])
            u2 = np.dot(x[k2:256].T, p[k2:256])
            W = np.array([[w0, w1, w2]]).T  # 注意是两层中括号，用np.array组成的向量可以转置,输出矩阵大小等等操作
            U = np.array([[int(u0), int(u1), int(u2)]]).T
            UT = np.sum(U)
            temp = (U - UT * W) * (U - UT * W) / W
            varbetween = np.sum(temp)
            if varbetween>maxvar:
                cnt=cnt+1
                maxvar=varbetween
                th1=k1
                th2=k2
                # print("第"+str(cnt)+"次找到最佳阈值："+str(th1)+","+str(th2)+"此时类间方差为"+str(maxvar))
    # mean1 = int(np.dot(x[0:th1].T, p[0:th1]))
    # mean2 = int(np.dot(x[th1:th2].T, p[th1:th2]))
    # mean3 = int(np.dot(x[th2:256].T, p[th2:256]))
    for row in range(0,high):
        for col in range(0,wide):
            if img[row, col]=th1 and img[row, col]<=th2:
                img[row, col] =127
            else:
                img[row, col] = 255
    cv.imshow("2th3gray",img)
    cv.waitKey()
    cv.destroyAllWindows()


if __name__ == '__main__':
    otsu3layers("sun.jpg")

感受

至此笔者一共花了一个星期五下午和星期六上下午来写这篇文章，真的收获了很多很多东西，python的plt画图、numpy的基础操作、向量化处理减少计算量、阅读原英文论文学习、思考模型简化模型、从理论到实践到敲出代码到看到图片，每个任务都是一个小小的挑战，做成这样真的非常开心??????????感谢！

后话

由于笔者没有实验器材（红外光源和滤光片等等器材），所以真正地识别到实际的红外光源可能还得添加层数，复杂度还将继续上升，这时可能要简化算法等等操作（网上另外有大佬用遗传算法求得），期待后续的研究。

to be continued...