【老文新发】Otsu大津法详解及python实现

原文:A Threshold Selection Method from Gray-Level Histograms
A Fast Algorithm for Multilevel Thresholding

前言

大津法包含两个重要的概念:类间方差(between-class variance)类内方差(within-class variance)

两者的详细关系推导可后文。

大津法又称为最大类间方差法是有原因的。因为这个算法的目的就是最大化类间方差,且这个最优阈值一定存在。

大津法作为阈值自动分割的经典算法,其思想很巧妙,值得学习。

大津法推导

在这里插入图片描述
如图所示,右边分割分明,其类间的差异大,区分明显,所以其类间方差更大。

大津法就是要实现这个过程。

我们先做图像的直方图,统计每个很小像素区间包含的像素个数。

即将图像像素值分为 [ 1 , 2 , 3 , …. , L] 个区间。用 n i n_{i} ni 表示各个水平像素值的像素个数,总像素个数与 n i n_{i} ni关系为:
N = n 1 + n 2 + … + n L N=n_1+n_2+\ldots+n_L N=n1+n2++nL
像元个数 n i ni ni 比上总像元数 N 即可得到某个像素区间出现的频率,定义 p i pi pi
p i = n i / N i , p i ≥ 0 , ∑ i = 1 L p i = 1 p_i=n_i / N_i, p_i \geq 0, \sum_{i=1}^L p_i=1 pi=ni/Ni,pi0,i=1Lpi=1
根据原文,是利用一个阈值 k 把图像分为两类 C 0 、 C 1 。 k ∈ [ 1 , 2 , 3 , … . , L ] C0 、C1。k∈[ 1 , 2 , 3 , …. , L] C0C1k[1,2,3,.,L]
我们分别求出这个阈值前后的局部频率之和,定义如下:
C 0 = [ 1 , k ] C 1 = [ k + 1 , L ] w 0 = Pr ⁡ ( C 0 ) = ∑ i = 1 k p i = w ( k ) w 1 = Pr ⁡ ( C 1 ) = ∑ i = k + 1 L p i = 1 − w ( k ) \begin{gathered} C_0=[1, k] \\ C_1=[k+1, L] \\ w_0=\operatorname{Pr}\left(C_0\right)=\sum_{i=1}^k p_i=w(k) \\ w_1=\operatorname{Pr}\left(C_1\right)=\sum_{i=k+1}^L p_i=1-w(k) \end{gathered} C0=[1,k]C1=[k+1,L]w0=Pr(C0)=i=1kpi=w(k)w1=Pr(C1)=i=k+1Lpi=1w(k)
则灰度图像频率直方图的总的数学期望和 C0 、C1的数学期望如下:
u 0 = ∑ i = 1 k i ∗ Pr ⁡ ( i ∣ C 0 ) = ∑ i = 1 k i ∗ p i / w 0 = u ( k ) w ( k ) u 1 = ∑ i = k + 1 L i ∗ Pr ⁡ ( i ∣ C 1 ) = ∑ i = k + 1 L i ∗ p i / w 1 = u T − u ( k ) 1 − w ( k ) u T = u ( L ) = ∑ i = 1 L i ∗ p i \begin{gathered} u_0=\sum_{i=1}^k i * \operatorname{Pr}\left(i \mid C_0\right)=\sum_{i=1}^k i * p_i / w_0=\frac{u(k)}{w(k)} \\ u_1=\sum_{i=k+1}^L i * \operatorname{Pr}\left(i \mid C_1\right)=\sum_{i=k+1}^L i * p_i / w_1=\frac{u_T-u(k)}{1-w(k)} \\ u_T=u(L)=\sum_{i=1}^L i * p_i \end{gathered} u0=i=1kiPr(iC0)=i=1kipi/w0=w(k)u(k)u1=i=k+1LiPr(iC1)=i=k+1Lipi/w1=1w(k)uTu(k)uT=u(L)=i=1Lipi
上式各个变量之间的关系如下:
w 0 u 0 + w 1 u 1 = u T w 0 + w 1 = 1 w_0 u_0+w_1 u_1=u_T \quad w_0+w_1=1 w0u0+w1u1=uTw0+w1=1
期望有了,计算一下对应的方差:
σ 0 2 = ∑ i = 1 k ( i − u 0 ) 2 Pr ⁡ ( i ∣ C 0 ) = ∑ i = 1 k ( i − u 0 ) 2 p i / w 0 σ 1 2 = ∑ i = k + 1 L ( i − u 1 ) 2 Pr ⁡ ( i ∣ C 1 ) = ∑ i = k + 1 L ( i − u 1 ) 2 p i / w 0 σ T 2 = ∑ i = 1 L ( i − u T ) 2 p i \begin{gathered} \sigma_0^2=\sum_{i=1}^k\left(i-u_0\right)^2 \operatorname{Pr}\left(i \mid C_0\right)=\sum_{i=1}^k\left(i-u_0\right)^2 p_i / w_0 \\ \sigma_1^2=\sum_{i=k+1}^L\left(i-u_1\right)^2 \operatorname{Pr}\left(i \mid C_1\right)=\sum_{i=k+1}^L\left(i-u_1\right)^2 p_i / w_0 \\ \sigma_T^2=\sum_{i=1}^L\left(i-u_T\right)^2 p_i \end{gathered} σ02=i=1k(iu0)2Pr(iC0)=i=1k(iu0)2pi/w0σ12=i=k+1L(iu1)2Pr(iC1)=i=k+1L(iu1)2pi/w0σT2=i=1L(iuT)2pi
根据文献:Introduction to statistical pattern recognition,260-267。类内误差、类间误差、总误差有如下关系:
类内误差 σ w 2 = w 0 σ 0 2 + w 1 σ 1 2 类内误差 \sigma_w^2=w_0 \sigma_0^2+w_1 \sigma_1^2 类内误差σw2=w0σ02+w1σ12
类间误差 σ b 2 = w 0 ( u 0 − u T ) 2 + w 1 ( u 1 − u T ) 2 类间误差\sigma_b^2=w_0\left(u_0-u_T\right)^2+w_1\left(u_1-u_T\right)^2 类间误差σb2=w0(u0uT)2+w1(u1uT)2

总误差 σ w 2 + σ b 2 = σ T 2 总误差\sigma_w^2+\sigma_b^2=\sigma_T^2 总误差σw2+σb2=σT2
注意:总误差是与阈值k无关的,但类间误差和类内误差是与阈值k相关的函数
σ b 2 = w 0 ( u 0 − u T ) 2 + w 1 ( u 1 − u T ) 2 = w 0 ( u 0 − ( w 0 u 0 + w 1 u 1 ) ) 2 + w 1 ( u 1 − ( w 0 u 0 + w 1 u 1 ) ) 2 = w 0 w 1 ( u 1 − u 0 ) 2 \begin{gathered} \sigma_b^2=w_0\left(u_0-u_T\right)^2+w_1\left(u_1-u_T\right)^2 \\ =w_0\left(u_0-\left(w_0 u_0+w_1 u_1\right)\right)^2+w_1\left(u_1-\left(w_0 u_0+w_1 u_1\right)\right)^2 \\ =w_0 w_1\left(u_1-u_0\right)^2 \end{gathered} σb2=w0(u0uT)2+w1(u1uT)2=w0(u0(w0u0+w1u1))2+w1(u1(w0u0+w1u1))2=w0w1(u1u0)2
然后再分别把w0,u1,u0带入,可得到:
σ b 2 = [ u T w ( k ) − u ( k ) ] 2 w ( k ) [ 1 − w ( k ) ] \sigma_b^2=\frac{\left[u_T w(k)-u(k)\right]^2}{w(k)[1-w(k)]} σb2=w(k)[1w(k)][uTw(k)u(k)]2

则求解最大类间误差为:
σ b 2 ( k ∗ ) = max ⁡ 1 ≤ k < L σ b 2 ( k ) \sigma_b^2\left(k^*\right)=\max _{1 \leq k<L} \sigma_b^2(k) σb2(k)=1k<Lmaxσb2(k)

由上述 σ b 2 \sigma_b^2 σb2 的分母可以发现, w ( k ) w(k) w(k) 可以取到 1 也可以取到0,因此在边界上 σ b 2 \sigma_b^2 σb2 可以无穷大,而在开基 ( 0 , 1 ) (0,1) (0,1) 则类间方差有限,因此在定义域
S ∗ = k : w 0 w 1 = w ( k ) [ 1 − w ( k ) ] > 0 S^*=k: w_0 w_1=w(k)[1-w(k)]>0 S=k:w0w1=w(k)[1w(k)]>0
因此必定存在一个阈值k使得两类类间方差最大。
以下是python代码实现:

def otsu(gray_img):n_count = gray_img.sizegray_img_array = gray_img.flatten()index = np.flatnonzero(gray_img_array)gray_img_data = gray_img_array[index]threshold_t = 0max_g = 0t = np.linspace(start=-1, stop=1, num=256)# 遍历每一个灰度层for i in range(len(t)):# 使用numpy直接对数组进行运算n0 = gray_img_data[np.where(gray_img_data < t[i])]n1 = gray_img_data[np.where(gray_img_data >= t[i])]w0 = len(n0) / n_countw1 = len(n1) / n_countu0 = np.mean(n0) if len(n0) > 0 else 0.u1 = np.mean(n1) if len(n0) > 0 else 0.g = w0 * w1 * (u0 - u1) ** 2if g > max_g:max_g = gthreshold_t = t[i]print('类间方差最大阈值:', threshold_t)gray_img[gray_img < threshold_t] = 0gray_img[gray_img >= threshold_t] = 1return gray_img

这个在opencv中已经有实现,可以直接调用

import cv2
t, otsu = cv2.threshold(img, 0, 255, cv2>THRESH_BINARY + cv2.THRESH_OTSU)

多分类最大类间方差法

根据以上公式类推到多分类的最大类间方差法,假设有 m − 1 \mathrm{m}-1 m1 个阈值 { t 1 , t 2 , … , t M − 1 } \{\mathrm{t} 1, \mathrm{t} 2, \ldots, \mathrm{tM}-1\} {t1,t2,,tM1} 将图像分为 M \mathrm{M} M 类, C 1 C_1 C1 C 2 … C M C_2 \ldots C_M C2CM 。 则存在一组阈值 { t 1 ∗ , t 2 ∗ , … , t M − 1 ∗ } \left\{t 1^*, \mathrm{t} 2^*, \ldots, \mathrm{tM}-1^*\right\} {t1,t2,,tM1} 使得
{ t 1 ∗ , t 2 ∗ , … , t M − 1 ∗ } = Arg ⁡ Max ⁡ { σ B 2 ( t 1 , t 2 , … , t M − 1 ) } , 1 ≤ t 1 < … < t M − 1 < L \begin{aligned} \left\{\mathrm{t}_1 *, \mathrm{t}_2 *, \ldots, \mathrm{t}_{\mathrm{M}-1} *\right\}= & \operatorname{Arg} \operatorname{Max}\left\{\sigma_{\mathrm{B}}{ }^2\left(\mathrm{t}_1, \mathrm{t}_2, \ldots, \mathrm{t}_{\mathrm{M}-1}\right)\right\}, \\ & 1 \leq \mathrm{t}_1<\ldots<\mathrm{t}_{\mathrm{M}-1}<\mathrm{L} \end{aligned} {t1,t2,,tM1}=ArgMax{σB2(t1,t2,,tM1)},1t1<<tM1<L

成立
其中:
σ B 2 = ∑ k = 1 M ω k ( μ k − μ T ) 2 ω k = ∑ i ∈ C k p i , μ k = ∑ i ∈ C k i p i / ω ( k ) . \begin{aligned} \sigma_{\mathrm{B}}{ }^2 & =\sum_{k=1}^{\mathrm{M}} \omega_{\mathrm{k}}\left(\mu_{\mathrm{k}}-\mu_{\mathrm{T}}\right)^2 \\ \omega_{\mathrm{k}} & =\sum_{\mathrm{i} \in \mathrm{Ck}} \mathrm{p}_{\mathrm{i}}, \\ \mu_{\mathrm{k}} & =\sum_{\mathrm{i} \in \mathrm{Ck}} \mathrm{i} \mathrm{p}_{\mathrm{i}} / \omega(\mathrm{k}) . \end{aligned} σB2ωkμk=k=1Mωk(μkμT)2=iCkpi,=iCkipi/ω(k).

因为
∑ k = 1 M ω k = 1 μ T = ∑ k = 1 M ω k μ k . \begin{gathered} \sum_{k=1}^{\mathrm{M}} \omega_{\mathrm{k}}=1 \\ \mu_{\mathrm{T}}=\sum_{k=1}^{\mathrm{M}} \omega_{\mathrm{k}} \mu_{\mathrm{k}} . \end{gathered} k=1Mωk=1μT=k=1Mωkμk.
因此
σ B 2 ( t 1 , t 2 , … , t M − 1 ) = ∑ k = 1 M ω k μ k 2 − μ T 2 \sigma_{\mathrm{B}}{ }^{2}\left(\mathrm{t}_{1}, \mathrm{t}_{2}, \ldots, \mathrm{t}_{\mathrm{M}-1}\right)=\sum_{k=1}^{\mathrm{M}} \omega_{\mathrm{k}} \mu_{\mathrm{k}}^{2}-\mu_{\mathrm{T}}{ }^{2} σB2(t1,t2,,tM1)=k=1Mωkμk2μT2

μ T \mu \mathrm{T} μT 与间值无关,因此求上式的最大值可转为:
{ t 1 ∗ , t 2 ∗ , … , t M − 1 ∗ } = Arg ⁡ Max ⁡ { ( σ B ′ ) 2 { { t 1 , t 2 , … , t M − 1 } } 1 ≤ t 1 < … < t M − 1 < L ( σ B ) 2 = ∑ k = 1 M ω k μ k 2 \begin{array}{c} \left\{\mathrm{t}_{1}^{*}, \mathrm{t}_{2}^{*}, \ldots, \mathrm{t}_{\mathrm{M}-1}^{*}\right\}=\operatorname{Arg} \operatorname{Max}\left\{\left(\sigma_{\mathrm{B}}{ }^{\prime}\right)^{2}\left\{\left\{\mathrm{t}_{1}, \mathrm{t}_{2}, \ldots, \mathrm{t}_{\mathrm{M}-1}\right\}\right\}\right. \\ 1 \leq \mathrm{t}_{1}<\ldots<\mathrm{t}_{\mathrm{M}-1}<\mathrm{L} \\ \left(\sigma_{\mathrm{B}}\right)^{2}=\sum_{k=1}^{\mathrm{M}} \omega_{\mathrm{k}} \mu_{\mathrm{k}}{ }^{2} \end{array} {t1,t2,,tM1}=ArgMax{(σB)2{{t1,t2,,tM1}}1t1<<tM1<L(σB)2=k=1Mωkμk2

类间方差、类内方差和总方差关系

https://blog.csdn.net/qq_42164483/article/details/119064535
https://blog.csdn.net/m0_38024332/article/details/104226806
以上这两篇博文讲的也很详细,内容有所参考,欢迎访问阅读。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.hqwc.cn/news/205373.html

如若内容造成侵权/违法违规/事实不符,请联系编程知识网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

CRM商机管理软件:构建客户为中心的管理理念

企业为什么选择CRM商机管理软件&#xff1f;1.CRM软件能够帮助企业建立以客户为中心的管理理念&#xff1b;2.CRM商机管理软件全面直观的展示客户数据&#xff1b;3.市场人员可以制订个性化的营销策略&#xff1b;4.移动应用为外出的销售带来的便利。 1.构建客户为中心的管理理…

2023年中国雷达设备市场规模及市场份额分析[图]

雷达设备行业是一种利用无线电波对目标进行探测和定位的技术&#xff0c;也被称为无线电探测和定位。雷达通过发射电磁波对目标进行照射并接收其回波&#xff0c;经波形处理后获取目标的位置和速度等信息。雷达具有探测距离远&#xff0c;测定精度高&#xff0c;不受天气和地形…

Motion Plan之搜索算法笔记

背景&#xff1a; 16-18年做过一阵子无人驾驶&#xff0c;那时候痴迷于移动规划&#xff1b;然而当时可学习的资料非常少&#xff0c;网上的论文也不算太多。基本就是Darpa的几十篇无人越野几次比赛的文章&#xff0c;基本没有成系统的文章和代码讲解实现。所以对移动规划的认…

迅为RK3568开发板学习之Linux驱动篇第十三期输入子系统

驱动视频全新升级&#xff0c;并持续更新~更全&#xff0c;思路更科学&#xff0c;入门更简单。 迅为基于iTOP-RK3568开发板进行讲解&#xff0c;本次更新内容为第十三期&#xff0c;主要讲解输入子系统&#xff0c;共计24 讲。 关注B站&#xff1a;北京迅为电子&#xff0c;在…

uboot移植-野火imx6ull

uboot移植及LCD和网络驱动设置 环境&#xff1a; Ubuntu 20.04 野火imx6ull pro开发板 tftp服务的搭建 IMX6ULL开发板参考的是NXP官方的I.MX6ULL EVK开发板做的硬件。 Linux的移植要复杂的多&#xff0c;在移植Linux之前我们需要先移植一个 bootloader 代码&#xff0c;这…

CPSC发布关于亚马逊含有纽扣电池或硬币电池产品的相关规则标准!UL4200A

2023年9月21日&#xff0c;美国消费品安全委员会&#xff08;CPSC&#xff09;在《联邦公报》上发布了纽扣及硬币电池及相关产品的最终规则&#xff08;DFR&#xff09;16 CFR 1263&#xff0c;以保护6岁以下儿童免受电池摄入危害。DFR将于2023年10月23日生效&#xff0c;除非消…

赛氪荣幸受邀参与中国联合国采购促进会第五次会员代表大会

11 月21 日 &#xff08;星期二&#xff09; 下午14:00&#xff0c;在北京市朝阳区定福庄东街1号中国传媒大学&#xff0c;赛氪荣幸参与中国联合国采购促进会第五次会员代表大会。 2022年以来&#xff0c;联合国采购杯全国大学生英语大赛已经走上了国际舞台&#xff0c;共有来自…

替换的DLL用户电脑报错加载失败

编译后混淆加签名的dll 远程下载下来有个选项&#xff1a; 在某用户电脑上出现加载失败的报错 右键dll 属性里勾选解除锁定后 加载运行正常 跟用户电脑安全策略有关系 有的会出现 大部分不会

铸就匠心,打造西部最具权威的行业商会组织

中国商报陕西报道&#xff08;记者 朱清平&#xff09;西安市五金机电商会(以下简称商会)第二届一次会员代表大会暨新任理事、监事就职典礼于11月17日在西安经开洲际酒店召开。 商会于2018年10月成立,在5年的发展中,依托“一带一路”发展的“快车道”,通过新丝路国际工业品数字…

前端性能优化总结

这里写目录标题 页面性能测试工具测试指标 前端页面性能常见的问题前端页面性能优化常见策略及方案dns优化------预解析域名&#xff08;异步进行&#xff09;http请求优化减少请求次数同时多开持久连接 前面面试url从输入到确认搜索发生了什么js介绍new一个对象的过程&#xf…

深度神经网络下的风格迁移模型

版权声明&#xff1a;本文为博主原创文章&#xff0c;转载请在显著位置标明本文出处以及作者网名&#xff0c;未经作者允许不得用于商业目的。 斯坦福大学李飞飞团队的风格迁移模型是一种基于深度学习的图像处理技术&#xff0c;可以将一张图像的风格转移到另一张图像上。该模型…

Spring框架学习 -- 核心思想

目录 (1) Spring是什么? (2) 什么是IOC容器? (3) 从传统开发认识spring (4) 这种传统开发的缺陷 (5)解决传统开发中的缺陷 (6) 对比总结规律 (7) 理解IOC 创作不易多多支持 (1) Spring是什么? 我们常说的Spring的全称是: Spring Framework(Spring框架), 它是一个开源…