当前位置: 首页 > news >正文

OTSU算法及其Python实现

文章目录

    • 原理
    • 实现和验证
    • 分析和优化

原理

OTSU算法是大津展之提出的阈值分割方法,又叫最大类间方差法。OTSU并不是一个英文缩写,而是日语假名,是其提出者的姓氏“大津”。

假设存在阈值 T T T可以将图像分为两部分,记二者均值为 m 1 , m 2 m_1, m_2 m1,m2,图像总均值为 m m m,像素被分入这两部分的比例分别为 p 1 , p 2 p_1, p_2 p1,p2。从而

p 1 + p 2 = 1 p 1 m 1 + p 2 m 2 = m \begin{aligned} p_1+p_2&=1\\ p_1m_1+p_2m_2&=m\\ \end{aligned} p1+p2p1m1+p2m2=1=m

则类间方差可表示为

σ 2 = p 1 ( m 1 − m ) 2 + p 2 ( m 2 − m ) 2 = p 1 m 1 2 + p 2 m 2 2 − m 2 \begin{aligned} \sigma^2&=p_1(m_1-m)^2+p_2(m_2-m)^2\\ &=p_1m_1^2+p_2m_2^2-m^2 \end{aligned} σ2=p1(m1m)2+p2(m2m)2=p1m12+p2m22m2

实现和验证

由于OTSU算法的逻辑和公式都非常清晰,所以实现起来也及其方便

import numpy as np
import matplotlib.pyplot as plt# 计算类间方差
def getVar(img, th):slct = img>thp1 = np.sum(slct)/img.sizeif p1 in [0, 1]:return 0p2 = 1-p1m1 = np.mean(img[slct])m2 = np.mean(img[~slct])m = p1*m1+p2*m2return p1*m1**2+p2*m2**2-m**2def otsu(img):vs = [getVar(img, th) for th in range(256)]return np.argmax(vs)

其中,getVar用于计算类间方差,后面的otsu则返回分割阈值。下面创建函数用于测试otsu算法

def drawOne(fig, index, img):fig.add_subplot(index)plt.imshow(img, cmap="gray")plt.axis('off')def test():img = plt.imread(r"D:\Code\NotePL\python\lena.jpg").astype(float)img = np.mean(img, axis=2)th = otsu(img)print(th)fig = plt.figure(figsize=(8,3))drawOne(fig, 131, img)drawOne(fig, 132, img>th)drawOne(fig, 133, img<=th)plt.tight_layout()plt.show()

效果如下

在这里插入图片描述

分析和优化

由于图像的像素值是八位整型,所以迭代256次就可以得到所有的类间方差。

img = plt.imread(r"D:\Code\NotePL\python\lena.jpg").astype(float)
img = np.mean(img, axis=2)
vs = [getVar(img, th) for th in range(256)]
print(np.argmax(vs))    # 121
plt.plot(vs)
plt.show()

类间方差分布如下,当阈值是121时,得到最大类间方差。

在这里插入图片描述

对于精度更高的16位图像,或者其他非图像的数值,遍历的方案效率太低了,为此可进行做一个步长二分的爬山算法,代码如下,最终得到的结果位121.7,由于在阈值分割时采用的是大于号,所以效果与121相同。

def climb(img, step, st=0, err=0.1):vSt = getVar(img, st)while abs(step)>err:ed = st+stepvEd = getVar(img, ed)if vEd < vSt:step = -step/2st, vSt = ed, vEdreturn edclimb(img, 40)
# 121.71875

相关文章:

  • python单引号怎么打,两种输入方法
  • 高德地图海量点标记图片(智慧交通道路事件情况)
  • 推荐10款App安全测试工具
  • [Firefly-Linux] RK3568 pca9555芯片驱动详解
  • AI跨界学习,不再是梦!
  • 宝塔是可以切换mongodb版本的
  • 用23种设计模式打造一个cocos creator的游戏框架----(四)装饰器模式
  • 用C语言实现链栈的基本操作
  • Day03 linux高级系统编程--进程
  • 电气火灾监控系统
  • 正则表达式(5):常用符号
  • 51单片机程序
  • 10_企业架构NOSQL数据库之MongoDB
  • Amazon CodeWhisperer 使用体验
  • 绝地求生:NH究极天命圈惊险吃鸡,17斩获单日积分第一,4AM梦游暂居倒数
  • 03Go 类型总结
  • 8年软件测试工程师感悟——写给还在迷茫中的朋友
  • CSS进阶篇--用CSS开启硬件加速来提高网站性能
  • DOM的那些事
  • Eureka 2.0 开源流产,真的对你影响很大吗?
  • Spring-boot 启动时碰到的错误
  • SQLServer插入数据
  • SSH 免密登录
  • thinkphp5.1 easywechat4 微信第三方开放平台
  • 后端_ThinkPHP5
  • 机器学习中为什么要做归一化normalization
  • 目录与文件属性:编写ls
  • 软件开发学习的5大技巧,你知道吗?
  • 深入体验bash on windows,在windows上搭建原生的linux开发环境,酷!
  • 一道面试题引发的“血案”
  • Mac 上flink的安装与启动
  • 正则表达式-基础知识Review
  • (2015)JS ES6 必知的十个 特性
  • (3)Dubbo启动时qos-server can not bind localhost22222错误解决
  • (WSI分类)WSI分类文献小综述 2024
  • (二十三)Flask之高频面试点
  • (七)Java对象在Hibernate持久化层的状态
  • (十八)三元表达式和列表解析
  • (详细版)Vary: Scaling up the Vision Vocabulary for Large Vision-Language Models
  • (一)Spring Cloud 直击微服务作用、架构应用、hystrix降级
  • (转)JVM内存分配 -Xms128m -Xmx512m -XX:PermSize=128m -XX:MaxPermSize=512m
  • (转)Oracle存储过程编写经验和优化措施
  • (最完美)小米手机6X的Usb调试模式在哪里打开的流程
  • .bat批处理(六):替换字符串中匹配的子串
  • .mkp勒索病毒解密方法|勒索病毒解决|勒索病毒恢复|数据库修复
  • .net 生成二级域名
  • .NET的微型Web框架 Nancy
  • ??eclipse的安装配置问题!??
  • @Autowired和@Resource的区别
  • @Import注解详解
  • [BZOJ1008][HNOI2008]越狱
  • [BZOJ2281][SDOI2011]黑白棋(K-Nim博弈)
  • [C#]获取指定文件夹下的所有文件名(递归)
  • [C和指针].(美)Kenneth.A.Reek(ED2000.COM)pdf
  • [flask]http请求//获取请求头信息+客户端信息