CN111028258B - 一种大尺度灰度图像自适应阈值提取方法 - Google Patents
一种大尺度灰度图像自适应阈值提取方法 Download PDFInfo
- Publication number
- CN111028258B CN111028258B CN201911112512.6A CN201911112512A CN111028258B CN 111028258 B CN111028258 B CN 111028258B CN 201911112512 A CN201911112512 A CN 201911112512A CN 111028258 B CN111028258 B CN 111028258B
- Authority
- CN
- China
- Prior art keywords
- value
- gray
- peak
- image
- 1stpk
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/10—Segmentation; Edge detection
- G06T7/136—Segmentation; Edge detection involving thresholding
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Image Analysis (AREA)
- Image Processing (AREA)
Abstract
本发明实施例公开了一种大尺度灰度图像自适应阈值提取方法，先矩阵化读入系列大尺度灰度图像，对符合要求的图像依次进行维纳滤波和高斯滤波处理，得到图像样品；然后统计图像样品的灰度值信息，得到图像样品的灰度直方分布图，读取系列图像的灰度图进行直方图的累加；再根据累加的直方分布进行求其梯度，得到图像样品不同的梯度分布数据；接着根据图像样品的灰度直方分布图及其梯度分布数据结果，进行灰度分布特征值求解，并由经验公式计算得到分割阈值；最后计算分割阈值的可调节范围，提供进行手动调节的范围限制。本发明的方法具有良好的自适应能力，全过程自动化进行，极大地提高了识别提取工作效率和识别准确度。
Description
技术领域
本发明实施例涉及数字图像处理技术领域，具体涉及一种具有显著四组分特征的大尺度灰度图像自适应阈值提取方法。
背景技术
图像阈值分割技术指的是根据不同阈值，将图像中不同组分进行区分的技术，传统的阈值分割主要有两类，一类为手动调节方法，通过目测识别，对图像进行不同分割阈值的试验性调节，得到最优分割阈值；另一类是以大津法为代表的自动阈值分割方法，基于最大类间差思想，进行图像的二组分分割。
手动调节方法有时候能够得到相对可信的阈值，但是需要反复分割试验的过程，这一过程存在太多的主观判断，也就带来了不准确性；另一方面，分割试验往往耗费较长时间，不具有经济性和操作性；同时，手动调节方法对于小尺度单张图像具有操作性，对于大尺度、具有系列图片的处理对象，根本无法执行准确调节。
基于最大类间差判断的大津法，具有相对较好的稳定性和可操作性，是进行二组分识别的重要方法，但是很难用其进行大尺度、四组分(抽象地指在图像中，具有四类可肉眼分类的显著成份的特征)的分析，主要原因是，随着矩阵增大，带来的计算量是非线性增长的，同时用二组分方法计算四组分特征的阈值，将带来反复迭代的过程，进一步增大的计算量，同时带来一定的不准确性。阈值分割仅仅是进行图像识别与提取的第一步，目前方法带来的计算量，甚至会超过识别提取本身。因此需要针对四组分图像，设计快速分割方法。
大尺度图像可以有两种形式，一种是单张图像，但是像素矩阵较大；另一种是单张图像像素矩阵很小，但是具有数目较多的系列图像。本方法对以上两种形式的处理，均能达到快速效果。
发明内容
为此，本发明实施例提供一种具有显著四组分特征的大尺度灰度图像自适应阈值提取方法，以解决现有技术中的问题。
为了实现上述目的，本发明的实施方式提供如下技术方案：一种大尺度灰度图像自适应阈值计算方法，包括如下步骤：
步骤100、图像前处理：矩阵化读入系列大尺度图像，并裁剪去除不符合要求的图像，然后对符合要求的图像依次进行维纳滤波和高斯滤波处理，得到图像样品；
步骤200、图像灰度直方分布提取：统计图像样品的灰度值信息，得到图像样品的灰度直方分布数据，再批量读取系列图像样品，并计算灰度直方分布数据的累加结果；
步骤300、根据累加的直方分布数据，通过求数据梯度的方法，得到图像样品对应的一阶梯度分布、二阶梯度分布和三阶梯度分布数据；
步骤400、特征点提取：根据图像样品的灰度直方分布数据及其一阶梯度分布、二阶梯度分布和三阶梯度分布数据，并进行特征值求解；
步骤500、分割阈值预测：根据求解的特征值，由经验公式计算得到分割阈值；
步骤600、计算分割阈值的可调节范围，提供进行手动调节的范围限制。
作为本发明一种优选的方案，所述图像样品为页岩样品，所述特征值包括页岩样品的有机质、无机质和黄铁矿的峰值对应灰度值，及其各自左侧、右侧的一阶极值、二阶极值和三阶极值，加上有机质、无机质之间的谷值，一共22个特征点。
作为本发明一种优选的方案，所述维纳滤波的方法包括：
生成一个矩阵大小为3×3或5×5的滤波模板，对目标点周围8或24个像素求取平均值和方差，目标点处于滤波模板中心位置，使用所述滤波模板内像素灰度的平均值和所述方差创建一个像素矩阵维纳滤波器，用所述维纳滤波器求滤波模板像素矩阵内灰度平均值作为滤波模板矩阵的中心值，也即得到目标点的滤波后值；
平均值为：
方差为：
滤波模板值为：
其中a是模板中(n1,n2)位置的像素灰度值，v2是噪声的标准偏差，或采用模板内估计的局部偏差值，b(n1,n2)。
作为本发明一种优选的方案，所述高斯滤波的方法包括：
根据高斯分布特征，生成大小为3×3或5×5的滤波模板，与维纳滤波相同，用滤波模板去乘以中心点近邻的与滤波模板相同大小区域的原图像灰度值，将结果作为中心点滤波后的值，例如当该值取0.8时，生成的模板为：
中心值由高斯平均值代替，对图像所有像素点依次求邻域的高斯平均值，则得到滤波后的图像。
作为本发明一种优选的方案，所述步骤200提取图像灰度直方分布特征的方法包括：
依次读入批量图像中的每一张图像，统计处于0-255各个灰度值的像素的个数，并将批量图像得到的灰度分布数据按照灰度值进行数量累加，得到样品图像总灰度直方分布数据；对总灰度直方分布分别求一阶梯度、二阶梯度和三阶梯度；提取总灰度直方分布及其一阶梯度分布数据、二阶梯度分布数据和三阶梯度分布数据的极大值和极小值，进行组合分析运算，自适应地得到进行灰度阈值分割三个特征阈值，从而进行显著四组分的划分，并根据图像自身特征，得到相应分割灰度的可调节上下限值。
作为本发明一种优选的方案，所述步骤300提取特征分布值的方法包括：
提取灰度直方分布的峰值，根据最大峰值，找到左侧、右侧的各两个极大值标记，共得到五个极大值，将其中最大峰值标记为为maxPk，右侧两个峰值标记为RPk1,RPk2，左侧两个峰值标记为LPk1,LPk2；
峰值整理，从提取得到的五个峰值中，取最大三个峰值所对应峰值点，峰值点从小至大，分别认为是第一峰值1stPk、第二峰值2ndPk和第三峰值3rdPk；判断如果最小的峰值过于小，接近零，那么只保留最大的两个峰值；
对于每个峰值，分别找到左、右侧分别对应的一阶导的极大值和极小值，进行分析；第一峰值的左、右两个极值分别标注为1stLSp,1stRSp；第二峰值的左右俩极值分别标注为2ndLSp，2ndRSp，第三峰值的左右两个极值分别标注为3rdLSp，3rdRSp，如果最小峰值对应为零，同样其左右极值也设定为零；
对于每个峰值，分别找左、右侧二阶梯度的极大值，作为拐点TF，每个峰值的TF分别标注为：2ndRTF，2ndLTF，1stRTF,1stLTF,3rdRTF,3rdLTF；
对于每个点的TF，分别找到左、右侧分别对应的三阶导的极大值和极小值，进行分析，分别标注为：2ndRKt，2ndLKt，1stRKt,1stLKt,3rdRKt,3rdLKt；
在1stPk和2ndPk之间找到灰度直方分布的谷值，标记为medvally；
当灰度图像中，没有第三个峰值时，根据已有数据构建一个第三峰的虚峰，以用于下一步计算，虚峰的峰点及其相关的点用经验公式给出：
3rdRPk＝floor(1stPk^(-2)*2+2ndLSp^(-2)*3+2ndPk^0.71*5.79+2ndRSp^(-2)*3)；
3rdLTF＝3rdRPk-14；
3rdLSp＝3rdRPk-10；
3rdRSp＝3rdRPk+10；
3rdRTF＝3rdRPk+14；
当灰度图像中，没有第三峰值时，根据已有的数据计算得到一个第三峰的虚峰，以用于下一步计算，第三峰的峰点可以用经验公式给：
1stPk＝round(2ndLSp^0.13*119.56-2ndPk^(-5.09)*0.31-2ndRSp^(-0.81)*8541.39)；
1stLSp＝round(1stPk^0.54*9.30+medvally^(-13.97)*0.20+2ndLSp^(0.75)*0.99-2ndPk^(-3.91)*2.89-2ndRSp^(-0.04)*86.37)；
1stRSp＝round(1stLSp^(-0.06)*0.50+1stPk^0.99*1.29+medvally^(-0.31)*0.21-2ndLSp^0.24*10.84+2ndPk^(-0.054)*50.83-2ndRSp^(-3.74)*26.14-3rdPk^1.09*0.034)；
1stLTF＝round(1stLSp^0.94*1.26-1stPk^(-6.81)*0.88-1stRSp^0.45*2.56-medvally^0.24*2.02+2ndPk^0.31*1.39+3rdPk^1.31*0.013)；
1stLKt＝round(1stLTF^1.15*0.23-1stLSp^(-0.0048)*0.77+1stPk^0.76*3.58-1stRSp^0.70*3.50+medvally^(-7.62)*0.63+2ndPk^0.85*0.26)；
1stRTF＝floor(1stPk+(1stPk-1stLTF)*(1stRSp-1stPk)^0.6/(1stPk-1stLSp)^0.6)。
作为本发明一种优选的方案，五个所述的极大值和五个所述的峰值均包含零值。
作为本发明一种优选的方案，所述步骤400中经验公式：
3rdlevel＝floor(-2ndLSp^0.85*3.78+2ndPk^0.76*6.83-2ndRSp^0.81*0.55+2ndRTF^0.82*3.22+3rdLTF^0.93*0.41-2ndRKt^1.00*0.53+3rdLSp^0.77*2.63-3rdPk^0.98*0.88)；
2ndlevel＝floor(1stPk^1.03*1.50-medvally^1.02*0.15+2ndLKt^1.09*0.73-2ndLTF^0.92*0.48-2ndLSp^0.98*3.32+2ndPk^0.94*6.50-2ndRSp^0.94*6.62+2ndRTF^0.98*2.39+2ndRKt^0.91*0.026)；
1stlevel＝floor(-1stLKt^(-1.48)*5.10+1stLTF^(2.37)*0.00069-1stLSp^0.30*54.92+1stPk^0.30*58.25-1stRSp^1.77*0.0032+medvally^0.31*10.41+2ndPk^(-0.40)*0.78)。
作为本发明一种优选的方案，所述步骤500中分割阈值的可调节范围的计算方法包括：
第一阈值的上下限，分别标记为u1stlevel和d1stlevel，计算经验公式为：
u1stlevel＝1stLTF；
u1stlevel(u1stlevel<＝1stlevel)＝1stlevel+5；
d1stlevel＝1stLKt-4；
d1stlevel(d1stlevel>＝1stlevel)＝1stlevel-5；
第二阈值的上下限制区间，分别标记为u2ndlevel和d2ndlevel，计算经验公式为：
u2ndlevel＝floor(medvally*0.5+2ndLKt*0.5)；
u2ndlevel(u2ndlevel<＝2ndlevel)＝2ndlevel+5；
d2ndlevel＝floor(0.5*1stRSp+0.5*1stRTF)；
d2ndlevel(d2ndlevel>＝2ndlevel)＝2ndlevel-5；
第三阈值的上下变化空间，分别标记为u3rdlevel和d3rdlevel，其计算经验公式为：
u3rdlevel＝floor(2ndRKt*0.6+3rdLSp*0.4)；
u3rdlevel(u3rdlevel<＝3rdlevel)＝3rdlevel+7；
d3rdlevel＝floor(2ndRSp*0.7+3rdLTF*0.3)；d3rdlevel(d3rdlevel>＝3rdlevel)＝3rdlevel-7。
本发明的实施方式具有如下优点：
本发明的方法具有良好的自适应能力，能够根据图像的特征，自动识别图像的灰度分布特征，从而根据具有四组分图像的本身特性，进行特征参数的提取，全过程自动化进行，避免了手动设置带来的低效率和不确定性，极大地提高了识别提取工作效率和识别准确度。
本发明能够对批量文件进行处理提取，能够对大尺度图像进行处理，得到这一系列图像(大尺度图像)所具有的共同灰度阈值，这比对单张图像进行处理结果具有更高的可信度。
附图说明
为了更清楚地说明本发明的实施方式或现有技术中的技术方案，下面将对实施方式或现有技术描述中所需要使用的附图作简单地介绍。显而易见地，下面描述中的附图仅仅是示例性的，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图引伸获得其它的实施附图。
本说明书所绘示的结构、比例、大小等，均仅用以配合说明书所揭示的内容，以供熟悉此技术的人士了解与阅读，并非用以限定本发明可实施的限定条件，故不具技术上的实质意义，任何结构的修饰、比例关系的改变或大小的调整，在不影响本发明所能产生的功效及所能达成的目的下，均应仍落在本发明所揭示的技术内容得能涵盖的范围内。
图1为本发明实施的算法流程图；
图2为本发明实施例页岩扫描电镜灰度图；
图3为本发明实施例页岩灰度直方分布图；
图4为本发明实施例页岩灰度直方分布一阶梯度图；
图5为本发明实施例页岩灰度直方分布二阶梯度图；
图6为本发明实施例页岩灰度直方分布三阶梯度图；
图7为本发明实施例页岩各组分提取结果图。
具体实施方式
以下由特定的具体实施例说明本发明的实施方式，熟悉此技术的人士可由本说明书所揭露的内容轻易地了解本发明的其他优点及功效，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1至图7所示，本发明提供了一种大尺度灰度图像自适应阈值计算方法，主要应用于具有显著四组分的图像组分识别与处理过程中，采用数字图像处理技术，具体为灰度直方分布统计的方法对图像特征进行自动分析与判别，进而自适应地得到图像的分割阈值。包括如下步骤：
步骤100、图像前处理：矩阵化读入系列大尺度图像，并裁剪去除不符合要求的图像，然后对符合要求的图像依次进行维纳滤波和高斯滤波处理，得到图像样品；
步骤200、图像灰度直方分布提取：统计图像样品的灰度值信息，得到图像样品的灰度直方分布数据，再批量读取系列图像样品，并计算灰度直方分布数据的累加结果；
步骤300、根据累加的直方分布数据，通过求数据梯度的方法，得到图像样品对应的一阶梯度分布、二阶梯度分布和三阶梯度分布数据；
步骤400、特征点提取：根据图像样品的灰度直方分布数据及其一阶梯度分布、二阶梯度分布和三阶梯度分布数据，并进行特征值求解；
步骤500、分割阈值预测：根据求解的特征值，由经验公式计算得到分割阈值；
步骤600、计算分割阈值的可调节范围，提供进行手动调节的范围限制。
在本实施方式中，图像样品可以为页岩样品，特征值包括页岩样品的有机质、无机质和黄铁矿的峰值对应灰度值，及其各自左侧、右侧的一阶极值、二阶极值和三阶极值，加上有机质、无机质之间的谷值，一共22个特征点。
在步骤100中，对样品图像进行滤波处理，先利用维纳(Wiener)自适应滤波方法进行处理，再用高斯(Gaussian)滤波进行。Wiener滤波是一种最佳滤波系统，可以用于提取被平稳噪声污染的信号。维纳滤波的方法包括：
生成一个矩阵大小为3×3或5×5的滤波模板，对目标点周围8或24个像素求取平均值和方差，目标点处于滤波模板中心位置，使用所述滤波模板内像素灰度的平均值和所述方差创建一个像素矩阵维纳滤波器，用所述维纳滤波器求滤波模板像素矩阵内灰度平均值作为滤波模板矩阵的中心值，也即得到目标点的滤波后值；
平均值为：
方差为：
滤波模板值为：
其中a是模板中(n1,n2)位置的像素灰度值，v2是噪声的标准偏差，或采用模板内估计的局部偏差值，b(n1,n2)。
在本步骤中，高斯滤波的方法具体为：
根据高斯分布特征，生成大小为3×3或5×5的滤波模板，与维纳滤波相同，用滤波模板去乘以中心点近邻的与滤波模板相同大小区域的原图像灰度值，将结果作为中心点滤波后的值，例如当该值取0.8时，生成的模板为：
中心值由高斯平均值代替，对图像所有像素点依次求邻域的高斯平均值，则得到滤波后的图像。
在步骤200中，提取图像灰度直方分布特征的具体方法如下：
依次读入批量图像中的每一张图像，统计处于0-255各个灰度值的像素的个数，并将批量图像得到的灰度分布数据按照灰度值进行数量累加，得到样品图像总灰度直方分布数据；对总灰度直方分布分别求一阶梯度、二阶梯度和三阶梯度；提取总灰度直方分布及其一阶梯度分布数据、二阶梯度分布数据和三阶梯度分布数据的极大值和极小值，进行组合分析运算，自适应地得到进行灰度阈值分割三个特征阈值，从而进行显著四组分的划分，并根据图像自身特征，得到相应分割灰度的可调节上下限值。
在步骤300提取特征分布值的方法具体为：
提取灰度直方分布的峰值，根据最大峰值，找到左侧、右侧的各两个极大值标记，共得到五个极大值(可以为零值)，将其中最大峰值标记为为maxPk，右侧两个峰值标记为RPk1,RPk2，左侧两个峰值标记为LPk1,LPk2；
峰值整理，从提取得到的五个峰值(可以为零值)中，取最大三个峰值所对应峰值点，峰值点从小至大，分别认为是第一峰值1stPk、第二峰值2ndPk和第三峰值3rdPk；判断如果最小的峰值过于小，接近零，那么只保留最大的两个峰值；
对于每个峰值，分别找到左、右侧分别对应的一阶导的极大值和极小值，进行分析；第一峰值的左、右两个极值分别标注为1stLSp,1stRSp；第二峰值的左右俩极值分别标注为2ndLSp，2ndRSp，第三峰值的左右两个极值分别标注为3rdLSp，3rdRSp，如果最小峰值对应为零，同样其左右极值也设定为零；
对于每个峰值，分别找左、右侧二阶梯度的极大值，作为拐点TF，每个峰值的TF分别标注为：2ndRTF，2ndLTF，1stRTF,1stLTF,3rdRTF,3rdLTF；
对于每个点的TF，分别找到左、右侧分别对应的三阶导的极大值和极小值，进行分析，分别标注为：2ndRKt，2ndLKt，1stRKt,1stLKt,3rdRKt,3rdLKt；
在1stPk和2ndPk之间找到灰度直方分布的谷值，标记为medvally；
当灰度图像中，没有第三个峰值时，根据已有数据构建一个第三峰的虚峰，以用于下一步计算，虚峰的峰点及其相关的点用经验公式给出：
3rdRPk＝floor(1stPk^(-2)*2+2ndLSp^(-2)*3+2ndPk^0.71*5.79+2ndRSp^(-2)*3)；
3rdLTF＝3rdRPk-14；
3rdLSp＝3rdRPk-10；
3rdRSp＝3rdRPk+10；
3rdRTF＝3rdRPk+14；
当灰度图像中，没有第三峰值时，根据已有的数据计算得到一个第三峰的虚峰，以用于下一步计算，第三峰的峰点可以用经验公式给：
1stPk＝round(2ndLSp^0.13*119.56-2ndPk^(-5.09)*0.31-2ndRSp^(-0.81)*8541.39)；
1stLSp＝round(1stPk^0.54*9.30+medvally^(-13.97)*0.20+2ndLSp^(0.75)*0.99-2ndPk^(-3.91)*2.89-2ndRSp^(-0.04)*86.37)；
1stRSp＝round(1stLSp^(-0.06)*0.50+1stPk^0.99*1.29+medvally^(-0.31)*0.21-2ndLSp^0.24*10.84+2ndPk^(-0.054)*50.83-2ndRSp^(-3.74)*26.14-3rdPk^1.09*0.034)；
1stLTF＝round(1stLSp^0.94*1.26-1stPk^(-6.81)*0.88-1stRSp^0.45*2.56-medvally^0.24*2.02+2ndPk^0.31*1.39+3rdPk^1.31*0.013)；
1stLKt＝round(1stLTF^1.15*0.23-1stLSp^(-0.0048)*0.77+1stPk^0.76*3.58-1stRSp^0.70*3.50+medvally^(-7.62)*0.63+2ndPk^0.85*0.26)；
1stRTF＝floor(1stPk+(1stPk-1stLTF)*(1stRSp-1stPk)^0.6/(1stPk-1stLSp)^0.6)。
步骤400中，由于已经得到了计算灰度分割阈值所需要的22个基本特征值，因此采用一定的经验公式进行图像灰度值预测，经验公式的获取，来自对3类不同样品图像，具有较大的代表性，具体经验公式如下：
3rdlevel＝floor(-2ndLSp^0.85*3.78+2ndPk^0.76*6.83-2ndRSp^0.81*0.55+2ndRTF^0.82*3.22+3rdLTF^0.93*0.41-2ndRKt^1.00*0.53+3rdLSp^0.77*2.63-3rdPk^0.98*0.88)；
2ndlevel＝floor(1stPk^1.03*1.50-medvally^1.02*0.15+2ndLKt^1.09*0.73-2ndLTF^0.92*0.48-2ndLSp^0.98*3.32+2ndPk^0.94*6.50-2ndRSp^0.94*6.62+2ndRTF^0.98*2.39+2ndRKt^0.91*0.026)；
1stlevel＝floor(-1stLKt^(-1.48)*5.10+1stLTF^(2.37)*0.00069-1stLSp^0.30*54.92+1stPk^0.30*58.25-1stRSp^1.77*0.0032+medvally^0.31*10.41+2ndPk^(-0.40)*0.78)。
步骤500中，根据上一步得到的基本分割阈值，计算得到分割阈值的可调节范围，也即推荐的准确度空间，具体根据经验公式得到，具体如下：
第一阈值的上下限，分别标记为u1stlevel和d1stlevel，计算经验公式为：
u1stlevel＝1stLTF；
u1stlevel(u1stlevel<＝1stlevel)＝1stlevel+5；
d1stlevel＝1stLKt-4；
d1stlevel(d1stlevel>＝1stlevel)＝1stlevel-5；
第二阈值的上下限制区间，分别标记为u2ndlevel和d2ndlevel，计算经验公式为：
u2ndlevel＝floor(medvally*0.5+2ndLKt*0.5)；
u2ndlevel(u2ndlevel<＝2ndlevel)＝2ndlevel+5；
d2ndlevel＝floor(0.5*1stRSp+0.5*1stRTF)；
d2ndlevel(d2ndlevel>＝2ndlevel)＝2ndlevel-5；
第三阈值的上下变化空间，分别标记为u3rdlevel和d3rdlevel，其计算经验公式为：
u3rdlevel＝floor(2ndRKt*0.6+3rdLSp*0.4)；
u3rdlevel(u3rdlevel<＝3rdlevel)＝3rdlevel+7；
d3rdlevel＝floor(2ndRSp*0.7+3rdLTF*0.3)；d3rdlevel(d3rdlevel>＝3rdlevel)＝3rdlevel-7。
在本实施方式中，对于图像处理的软件可以是MATLAB，也可以利用其他软件实现。
本发明的方法具有良好的自适应能力，能够根据图像的特征，自动识别图像的灰度分布特征，从而根据具有四组分图像的本身特性，进行特征参数的提取，全过程自动化进行，避免了手动设置带来的低效率和不确定性，极大地提高了识别提取工作效率和识别准确度。还具有良好的普适性，能够对岩石扫描电镜图像等典型的四组分特征的图像样品进行分割阈值提取。
另一方面，本发明能够对批量文件进行处理提取，能够对大尺度图像进行处理，得到这一系列图像(大尺度图像)所具有的共同灰度阈值，这比对单张图像进行处理结果具有更高的可信度。
本发明的方法还具有极快的计算速度，在基于MatLab的统计与计算中，主要时间消耗在于提取图像的灰度直方分布，对于一共50张图像的样本，识别时间大约是2.7秒；对于一共150张图像的样本，识别提取时间大约是6.8秒，对于一共305张图像的样本，识别提取时间大约是13.5秒。
本发明的方法具有较好的鲁棒性，对于四组分分布比较清晰可见的图像，能够较为准确地得三个分割阈值，并给出各个结果可调节的阈值范围。
虽然，上文中已经用一般性说明及具体实施例对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。
Claims (7)
1.一种大尺度灰度图像自适应阈值计算方法，其特征在于，包括如下步骤：
步骤100、图像前处理：矩阵化读入系列大尺度图像，并裁剪去除不符合要求的图像，然后对符合要求的图像依次进行维纳滤波和高斯滤波处理，得到图像样品；
步骤200、图像灰度直方分布提取：统计图像样品的灰度值信息，得到图像样品的灰度直方分布数据，再批量读取系列图像样品，并计算灰度直方分布数据的累加结果；
步骤300、根据累加的直方分布数据，通过求数据梯度的方法，得到图像样品对应的一阶梯度分布、二阶梯度分布和三阶梯度分布数据；
所述步骤300提取特征分布值的方法包括：
提取灰度直方分布的峰值，根据最大峰值，找到左侧、右侧的各两个极大值标记，共得到五个极大值，将其中最大峰值标记为为maxPk，右侧两个峰值标记为RPk1,RPk2，左侧两个峰值标记为LPk1,LPk2；
峰值整理，从提取得到的五个峰值中，取最大三个峰值所对应峰值点，峰值点从小至大，分别认为是第一峰值1stPk、第二峰值2ndPk和第三峰值3rdPk；判断如果最小的峰值过于小，接近零，那么只保留最大的两个峰值；
对于每个峰值，分别找到左、右侧分别对应的一阶导的极大值和极小值，进行分析；第一峰值的左、右两个极值分别标注为1stLSp,1stRSp；第二峰值的左右俩极值分别标注为2ndLSp，2ndRSp，第三峰值的左右两个极值分别标注为3rdLSp，3rdRSp，如果最小峰值对应为零，同样其左右极值也设定为零；
对于每个峰值，分别找左、右侧二阶梯度的极大值，作为拐点TF，每个峰值的TF分别标注为：2ndRTF，2ndLTF，1stRTF,1stLTF,3rdRTF ,3rdLTF；
对于每个点的TF，分别找到左、右侧分别对应的三阶导的极大值和极小值，进行分析，分别标注为：2ndRKt，2ndLKt，1stRKt, 1stLKt, 3rdRKt ,3rdLKt；
在1stPk和2ndPk之间找到灰度直方分布的谷值，标记为medvally；
当灰度图像中，没有第三个峰值时，根据已有数据构建一个第三峰的虚峰，以用于下一步计算，虚峰的峰点及其相关的点用经验公式给出：
3rdRPk=floor(1stPk^(-2)*2+2ndLSp^(-2)*3+2ndPk^0.71*5.79+2ndRSp^(-2)*3)；
3rdLTF=3rdRPk-14；
3rdLSp=3rdRPk-10；
3rdRSp=3rdRPk+10；
3rdRTF=3rdRPk+14；
1stPk=round(2ndLSp^0.13*119.56-2ndPk^(-5.09)*0.31-2ndRSp^(-0.81)*8541.39)；
1stLSp=round(1stPk^0.54*9.30+medvally^(-13.97)*0.20+2ndLSp^(0.75)*0.99-2ndPk^(-3.91)*2.89-2ndRSp^(-0.04)*86.37)；
1stRSp=round(1stLSp^(-0.06)*0.50+1stPk^0.99*1.29+medvally^(-0.31)*0.21-2ndLSp^0.24*10.84+2ndPk^(-0.054)*50.83-2ndRSp^(-3.74)*26.14-3rdPk^1.09*0.034)；
1stLTF=round(1stLSp^0.94*1.26-1stPk^(-6.81)*0.88-1stRSp^0.45*2.56-medvally^0.24*2.02+2ndPk^0.31*1.39+3rdPk^1.31*0.013);
1stLKt=round(1stLTF^1.15*0.23-1stLSp^(-0.0048)*0.77+1stPk^0.76*3.58-1stRSp^0.70*3.50+medvally^(-7.62)*0.63+2ndPk^0.85*0.26);
1stRTF=floor(1stPk+(1stPk-1stLTF)*(1stRSp-1stPk)^0.6/(1stPk-1stLSp)^0.6)；
步骤400、特征点提取：根据图像样品的灰度直方分布数据及其一阶梯度分布、二阶梯度分布和三阶梯度分布数据，进行特征值求解；
步骤500、分割阈值预测：根据求解的特征值，由经验公式计算得到分割阈值；
所述步骤500中经验公式：
3rdlevel=floor(-2ndLSp^0.85*3.78+2ndPk^0.76*6.83-2ndRSp^0.81*0.55+2ndRTF^0.82*3.22+3rdLTF^0.93*0.41-2ndRKt^1.00*0.53+3rdLSp^0.77*2.63-3rdPk^0.98*0.88);
2ndlevel=floor(1stPk^1.03*1.50-medvally^1.02*0.15+2ndLKt^1.09*0.73-2ndLTF^0.92*0.48-2ndLSp^0.98*3.32+2ndPk^0.94*6.50-2ndRSp^0.94*6.62+2ndRTF^0.98*2.39+2ndRKt^0.91*0.026);
1stlevel=floor(-1stLKt^(-1.48)*5.10+1stLTF^(2.37)*0.00069-1stLSp^0.30*54.92+1stPk^0.30*58.25-1stRSp^1.77*0.0032+medvally^0.31*10.41+2ndPk^(-0.40)*0.78)；
步骤600、计算分割阈值的可调节范围，提供进行手动调节的范围限制。
2.根据权利要求1所述的一种大尺度灰度图像自适应阈值提取方法，其特征在于，所述图像样品为页岩样品，所述特征值包括页岩样品的有机质、无机质和黄铁矿的峰值对应灰度值，及其各自左侧、右侧的一阶极值、二阶极值和三阶极值，加上有机质、无机质之间的谷值，一共22个特征点。
3.根据权利要求1所述的一种大尺度灰度图像自适应阈值提取方法，其特征在于，所述维纳滤波的方法包括：
生成一个矩阵大小为3×3或5×5的滤波模板，对目标点周围8或24个像素求取平均值和方差，目标点处于滤波模板中心位置，使用所述滤波模板内像素灰度的平均值和所述方差创建一个像素矩阵维纳滤波器，用所述维纳滤波器求滤波模板像素矩阵内灰度平均值作为滤波模板矩阵的中心值，也即得到目标点的滤波后值；
平均值为：；
方差为： ；
滤波模板值为： ；
其中是模板中位置的像素灰度值，是噪声的标准偏差，或采用模板内估计的局部偏差值，。
4.根据权利要求1所述的一种大尺度灰度图像自适应阈值提取方法，其特征在于，所述高斯滤波的方法包括：
根据高斯分布特征，生成大小为3×3或5×5的滤波模板，与维纳滤波相同，用滤波模板去乘以中心点近邻的与滤波模板相同大小区域的原图像灰度值，将结果作为中心点滤波后的值，当该值取0.8时，生成的模板为：；
中心值由高斯平均值代替，对图像所有像素点依次求邻域的高斯平均值，则得到滤波后的图像。
5.根据权利要求1所述的一种大尺度灰度图像自动阈值提取方法，其特征在于，所述步骤200提取图像灰度直方分布特征的方法包括：
依次读入批量图像中的每一张图像，统计处于0-255各个灰度值的像素的个数，并将批量图像得到的灰度分布数据按照灰度值进行数量累加，得到样品图像总灰度直方分布数据；对总灰度直方分布分别求一阶梯度、二阶梯度和三阶梯度；提取总灰度直方分布及其一阶梯度分布数据、二阶梯度分布数据和三阶梯度分布数据的极大值和极小值，进行组合分析运算，自适应地得到进行灰度阈值分割三个特征阈值，从而进行显著四组分的划分，并根据图像自身特征，得到相应分割灰度的可调节上下限值。
6.根据权利要求1所述的一种大尺度灰度图像自适应阈值提取方法，其特征在于，五个所述的极大值和五个所述的峰值均包含零值。
7.根据权利要求2所述的一种大尺度灰度图像自适应阈值提取方法，其特征在于，所述步骤600中分割阈值的可调节范围的计算方法包括：
第一阈值的上下限，分别标记为u1stlevel和d1stlevel，计算经验公式为：
u1stlevel=1stLTF；
u1stlevel(u1stlevel<=1stlevel)=1stlevel+5；
d1stlevel=1stLKt-4；
d1stlevel(d1stlevel>=1stlevel)=1stlevel-5；
第二阈值的上下限制区间，分别标记为u2ndlevel 和d2ndlevel，计算经验公式为：
u2ndlevel=floor(medvally*0.5+2ndLKt*0.5)；
u2ndlevel(u2ndlevel<=2ndlevel)=2ndlevel+5；
d2ndlevel=floor(0.5*1stRSp+0.5*1stRTF)；
d2ndlevel(d2ndlevel>=2ndlevel)=2ndlevel-5；
第三阈值的上下变化空间，分别标记为u3rdlevel和d3rdlevel，其计算经验公式为：
u3rdlevel=floor(2ndRKt*0.6+3rdLSp*0.4)；
u3rdlevel(u3rdlevel<=3rdlevel)=3rdlevel+7；
d3rdlevel=floor(2ndRSp*0.7+3rdLTF*0.3)； d3rdlevel(d3rdlevel>=3rdlevel)=3rdlevel-7。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911112512.6A CN111028258B (zh) | 2019-11-14 | 2019-11-14 | 一种大尺度灰度图像自适应阈值提取方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911112512.6A CN111028258B (zh) | 2019-11-14 | 2019-11-14 | 一种大尺度灰度图像自适应阈值提取方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111028258A CN111028258A (zh) | 2020-04-17 |
| CN111028258B true CN111028258B (zh) | 2023-05-16 |
Family
ID=70200138
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911112512.6A Active CN111028258B (zh) | 2019-11-14 | 2019-11-14 | 一种大尺度灰度图像自适应阈值提取方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111028258B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115661135B (zh) * | 2022-12-09 | 2023-05-05 | 山东第一医科大学附属省立医院(山东省立医院) | 一种心脑血管造影的病灶区域分割方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6335980B1 (en) * | 1997-07-25 | 2002-01-01 | Arch Development Corporation | Method and system for the segmentation of lung regions in lateral chest radiographs |
| CN105261017A (zh) * | 2015-10-14 | 2016-01-20 | 长春工业大学 | 基于路面约束的图像分割法提取行人感兴趣区域的方法 |
| CN109509199A (zh) * | 2018-10-10 | 2019-03-22 | 华南理工大学 | 一种基于三维重建的医学影像组织智能分割方法 |
| CN110033458A (zh) * | 2019-03-12 | 2019-07-19 | 中国矿业大学 | 一种基于像素梯度分布的图像阈值确定方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN1213592C (zh) * | 2001-07-31 | 2005-08-03 | 佳能株式会社 | 采用自适应二值化的图象处理方法和设备 |
| US8861814B2 (en) * | 2010-12-22 | 2014-10-14 | Chevron U.S.A. Inc. | System and method for multi-phase segmentation of density images representing porous media |
-
2019
- 2019-11-14 CN CN201911112512.6A patent/CN111028258B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6335980B1 (en) * | 1997-07-25 | 2002-01-01 | Arch Development Corporation | Method and system for the segmentation of lung regions in lateral chest radiographs |
| CN105261017A (zh) * | 2015-10-14 | 2016-01-20 | 长春工业大学 | 基于路面约束的图像分割法提取行人感兴趣区域的方法 |
| CN109509199A (zh) * | 2018-10-10 | 2019-03-22 | 华南理工大学 | 一种基于三维重建的医学影像组织智能分割方法 |
| CN110033458A (zh) * | 2019-03-12 | 2019-07-19 | 中国矿业大学 | 一种基于像素梯度分布的图像阈值确定方法 |
Non-Patent Citations (3)
| Title |
|---|
| 利用多重-大津阈值算法和扫描电镜分割CT图像;高衍武;吴伟;张虔;赵燕红;邵广辉;李国利;毛超杰;;长江大学学报(自然科学版);16(04);全文 * |
| 基于DTM阈值分割法的孔裂隙煤岩体瓦斯渗流数值模拟;王刚;杨鑫祥;张孝强;李文鑫;史林肯;;岩石力学与工程学报;35(01);全文 * |
| 基于图像像素间空间信息的加权模糊阈值分割算法;于勇;郑姣;郭希娟;;计算机应用与软件;30(03);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111028258A (zh) | 2020-04-17 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107507173B (zh) | 一种全切片图像的无参考清晰度评估方法及系统 | |
| CN107038416B (zh) | 一种基于二值图像改进型hog特征的行人检测方法 | |
| CN109460722B (zh) | 一种车牌智能识别方法 | |
| CN107481374B (zh) | 一种智能终端指纹解锁开门装置 | |
| CN111145209A (zh) | 一种医学图像分割方法、装置、设备及存储介质 | |
| CN103175844A (zh) | 一种金属零部件表面划痕缺陷检测方法 | |
| CN109815762B (zh) | 远距离识别二维码的方法、存储介质 | |
| CN112614062A (zh) | 菌落计数方法、装置及计算机存储介质 | |
| CN116630813B (zh) | 一种公路路面施工质量智能检测系统 | |
| CN110687122A (zh) | 一种陶瓦表面裂纹检测方法及系统 | |
| CN111429372A (zh) | 一种增强低对比度图像边缘检测效果的方法 | |
| CN110555373A (zh) | 一种基于图像识别的混凝土振捣质量实时检测方法 | |
| CN115731493A (zh) | 基于视频图像识别的降水微物理特征参量提取与分析方法 | |
| CN111028258B (zh) | 一种大尺度灰度图像自适应阈值提取方法 | |
| CN115601379A (zh) | 一种基于数字图像处理的表面裂纹精确检测技术 | |
| CN116958503B (zh) | 一种基于图像处理的污泥干化等级识别方法及系统 | |
| Shi et al. | Image enhancement for degraded binary document images | |
| CN107704864B (zh) | 基于图像对象性语义检测的显著目标检测方法 | |
| CN109829511B (zh) | 基于纹理分类的下视红外图像中云层区域检测方法 | |
| CN108319927B (zh) | 一种自动识别病害的方法 | |
| CN116433978A (zh) | 一种高质量瑕疵图像自动生成与自动标注方法及装置 | |
| CN114627463B (zh) | 一种基于机器识别的非接触式配电数据识别方法 | |
| CN113643290B (zh) | 一种基于图像处理的吸管计数方法、装置及存储介质 | |
| Bi et al. | Adaptive blind image restoration algorithm of degraded image | |
| CN114758139A (zh) | 基坑积水检测方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |