[go: up one dir, main page]

CN109255764A - 一种去除fib-sem图像幕帘噪声的方法及装置 - Google Patents

一种去除fib-sem图像幕帘噪声的方法及装置 Download PDF

Info

Publication number
CN109255764A
CN109255764A CN201811013158.7A CN201811013158A CN109255764A CN 109255764 A CN109255764 A CN 109255764A CN 201811013158 A CN201811013158 A CN 201811013158A CN 109255764 A CN109255764 A CN 109255764A
Authority
CN
China
Prior art keywords
image
noise
sectioning
tonal range
curtain
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.)
Pending
Application number
CN201811013158.7A
Other languages
English (en)
Inventor
孙朗秋
赵建国
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum Beijing
Original Assignee
China University of Petroleum Beijing
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China University of Petroleum Beijing filed Critical China University of Petroleum Beijing
Priority to CN201811013158.7A priority Critical patent/CN109255764A/zh
Publication of CN109255764A publication Critical patent/CN109255764A/zh
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00Image enhancement or restoration
    • G06T5/70Denoising; Smoothing
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/11Region-based segmentation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00Image analysis
    • G06T7/10Segmentation; Edge detection
    • G06T7/136Segmentation; Edge detection involving thresholding
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/10Image acquisition modality
    • G06T2207/10056Microscopic image
    • G06T2207/10061Microscopic image from scanning electron microscope
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20024Filtering details
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00Indexing scheme for image analysis or image enhancement
    • G06T2207/20Special algorithmic details
    • G06T2207/20048Transform domain processing
    • G06T2207/20056Discrete and fast Fourier transform, [DFT, FFT]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Image Processing (AREA)

Abstract

本发明提供了一种去除FIB‑SEM图像幕帘噪声的方法及装置,涉及图像处理技术领域。方法包括:将待处理的FIB‑SEM图像序列按主要矿物组分变化分为若干组,每组中选择关键切片图像划分不同灰度范围区域,形成每个灰度范围区域的带阻滤波器,并确定各带阻滤波器的参数;得到各带阻滤波器对应的滤波后的图像;形成关键切片图像对应的降噪处理图像;降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数,应用于关键切片图像所在的切片图像组的首尾两个切片图像,若降噪处理效果合格,将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;重复上述过程直至获得FIB‑SEM图像序列所有切片图像的降噪处理图像。

Description

一种去除FIB-SEM图像幕帘噪声的方法及装置
技术领域
本发明涉及图像处理技术领域,尤其涉及一种去除FIB-SEM图像幕帘噪声的方法及装置。
背景技术
当前,为满足人们对清洁能源的大量需求,非常规油气的勘探和开发已成为我国能源发展战略的重点。非常规油气储层的岩石具有结构复杂、致密、非均质性强的特点,因此,有必要利用聚焦离子束—扫描电子显微镜(Focused Ion Beam-ScanningElectronicMicroscope,简称FIB-SEM)对这类岩石样品进行数字岩心研究。FIB-SEM测试是用聚集离子束切削样品表面,形成平整的观察面,再用扫描电子显微镜拍摄纳米级分辨率的样品图像,如此反复,形成一套连续的切片(如图1所示,一组连续的FIB-SEM图像序列(示意图),被测试样品为富有机质的页岩。),可以组合成关于样品内部微观结构的三维图像。在FIB-SEM图像中(尤其是背散射图像),灰度与矿物的平均原子量正相关,即与密度正相关,高密度矿物的亮度高。因此,可利用灰度的变化来区分不同矿物。
非常规储层岩石样品多具有强非均质性。在FIB-SEM测试中,利用预设的固定功率的聚焦离子束切削时,样品表面的硬度的变化引起切削效率变化,在样品表面造成平行于切削方向(多为竖向)的划痕,从而在FIB-SEM图像中产生明暗变化的竖条纹噪声(如图2所示,一幅富有机质的页岩样品的FIB-SEM图像切片,有明显的幕帘噪声。中间部分的点状是利用软件自动提取的孔隙结构,受幕帘噪声影响,出现大量假的针状孔。)。由于形似悬挂的幕帘褶皱,通常称为“幕帘噪声(curtain noise)”。它由岩石样品本身的特性和FIB-SEM测试方式引起,不能完全避免,尤其是含有硬度偏低的矿物时。幕帘噪声对于岩石图像分析、基于图像分析的岩石性质提取和地质解释、岩石物理建模及岩石内部物理化学过程的模拟计算等都非常有害。例如,在图像的自动化处理中,幕帘噪声可能被识别为针状孔(如图2所示),导致对成岩过程的误判(例如误判为有孔隙水的淋滤作用),并在建模、模拟时产生虚假的各向异性。
在目前的实践中,现有技术中多采用均值滤波、基于二维傅立叶变换的低通滤波或带阻滤波来降低FIB-SEM图像的幕帘噪声,但效果较差。均值滤波和低通滤波相当于对图像进行模糊化,不能有效去除幕帘噪声,还会严重损失图像细节。带阻滤波可能存在参数冲突的问题,即对于同一幅FIB-SEM图像切片,同一套参数在部分区域滤波不足,在另一区域滤波过度并产生伪影。由于幕帘噪声形似遥感图像中的条带噪声,而针对后者已有多种成熟有效的方法,近年部分研究者试图借用遥感中的去条带(de-stripe)方法来去除FIB-SEM图像的幕帘噪声。应注意到,FIB-SEM图像的幕帘噪声与遥感中的条带噪声有本质区别:幕帘噪声是测试中实际存在过的划痕,其强度受聚焦离子束功率和矿物硬度的影响,难以预知;条带噪声是电磁信号干扰,空间分布上规律性较强,沿条带方向强度一定。因此,上述现有技术的这些方法用于去除幕帘噪声的效果多不理想,易产生严重的滤波伪影。个别方法(例如Münch et al,2009)即使有一些效果,但需结合傅立叶变换、小波变换等,且需要人工试验并设置复杂的滤波参数,计算量较大。一套FIB-SEM图像往往有几百幅切片,依靠人工逐张试验、设置复杂的参数在实际工作中缺乏可行性。可见,如何实现一种针对岩石的FIB-SEM图像的有效、快速、节省人力的幕帘噪声去除方法,并避免产生滤波伪影,获得细节清晰、结构真实的图像,成为了一个亟待解决的问题。
发明内容
本发明的实施例提供一种去除FIB-SEM图像幕帘噪声的方法及装置,以实现一种针对岩石的FIB-SEM图像的有效、快速、节省人力的幕帘噪声去除方法,并避免产生滤波伪影,获得细节清晰、结构真实的图像。
为达到上述目的,本发明采用如下技术方案:
一种去除FIB-SEM图像幕帘噪声的方法,包括:
步骤101、根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致。
步骤102、在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著。
步骤103、根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值。
步骤104、对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数。
步骤105、用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像。
步骤106、将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像。
步骤107、将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数。
步骤108、将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像。
步骤109、将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;
重复步骤101至109,直至步骤110、获得FIB-SEM图像序列所有切片图像的降噪处理图像。
具体的,所述对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数,包括:
根据公式:将部分图像数据f(x,y)进行二维傅里叶变换,确定部分图像数据f(x,y)对应的频谱F(u,v);其中,f(x,y)表示进行截取处理后获得的含有幕帘噪声的灰度范围区域的部分图像数据,f(x,y)的大小为M×N像素,x和y分别表示横坐标像素和纵坐标像素;i为虚数单位,u为频谱F(u,v)的横坐标空间频率,u=–M/2,–M/2+1,…,0,1,…,M/2–1;v为频谱F(u,v)的纵坐标空间频率,v=–N/2,–N/2+1,…,0,1,…,N/2–1;M和N为偶数;
根据频谱F(u,v)得到振幅谱|F(u,v)|,并在振幅谱|F(u,v)|中检测得到幕帘噪声所在的空间频率范围;
在幕帘噪声所在的空间频率范围内压制振幅,得到二维带阻滤波器H(u,v);
所述二维带阻滤波器H(u,v)的通带和阻带之间设置有过渡带,过渡带函数为其中,L为过渡带宽度;j=1,2,…,L,为相对于L的位置;该过渡带函数在下降区间时为在上升区间时为
将幕帘噪声所在的空间频率范围和过渡带宽度L作为二维带阻滤波器H(u,v)的参数。
具体的,所述用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像,包括:
根据公式:
对关键切片图像进行处理,得到各含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像gn(x,y);
其中,g(x,y)表示关键切片图像,g(x,y)的大小为X×Y像素,X和Y为偶数;gn(x,y)表示g(x,y)中的第n个含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像;Hn(u,v)表示第n个含有幕帘噪声的灰度范围区域对应的二维带阻滤波器。
具体的,所述将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像,包括:
根据平滑和灰度分割法,对关键切片图像进行平滑处理,然后利用灰度分割将各滤波器对应区域的位置提取出来;
在平滑后的灰度图像gs(x,y)中,获得二维带阻滤波器Hn(u,v)对应的含有幕帘噪声的灰度范围区域的灰度值区间:t1<gs(x,y)<t2;t1和t2为灰度阈值;
根据公式:将原图中二维带阻滤波器Hn(u,v)对应的幕帘噪声的灰度范围区域替换为Hn(u,v)对应的滤波后的图像gn(x,y)的相同区域,形成替换后的灰度图像g′n(x,y);
将各滤波器对应的灰度范围区域替换后,形成关键切片图像对应的降噪处理图像。
一种去除FIB-SEM图像幕帘噪声的装置,包括:
切片图像分组单元,用于根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致;
关键切片图像选择单元,用于在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著;
区域划分单元,用于根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值;
截取处理单元,用于对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数;
滤波处理单元,用于用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像;
图像替换单元,用于将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像;
区域及参数记录单元,用于将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数;
降噪处理效果测试单元,用于将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像;
降噪处理图像获得单元,用于将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;
重复切片图像分组单元至降噪处理图像获得单元的过程,直至获得FIB-SEM图像序列所有切片图像的降噪处理图像。
另外,所述截取处理单元,具体用于:
对每一个含有幕帘噪声的灰度范围区域进行截取处理,以获得每一个含有幕帘噪声的灰度范围区域的部分图像数据;
根据公式:将部分图像数据f(x,y)进行二维傅里叶变换,确定部分图像数据f(x,y)对应的频谱F(u,v);其中,f(x,y)表示进行截取处理后获得的含有幕帘噪声的灰度范围区域的部分图像数据,f(x,y)的大小为M×N像素,x和y分别表示横坐标像素和纵坐标像素;i为虚数单位,u为频谱F(u,v)的横坐标空间频率,u=–M/2,–M/2+1,…,0,1,…,M/2–1;v为频谱F(u,v)的纵坐标空间频率,v=–N/2,–N/2+1,…,0,1,…,N/2–1;M和N为偶数;
根据频谱F(u,v)得到振幅谱|F(u,v)|,并在振幅谱|F(u,v)|中检测得到幕帘噪声所在的空间频率范围;
在幕帘噪声所在的空间频率范围内压制振幅,得到二维带阻滤波器H(u,v);
所述二维带阻滤波器H(u,v)的通带和阻带之间设置有过渡带,过渡带函数为其中,L为过渡带宽度;j=1,2,…,L,为相对于L的位置;该过渡带函数在下降区间时为在上升区间时为
将幕帘噪声所在的空间频率范围和过渡带宽度L作为二维带阻滤波器H(u,v)的参数。
另外,所述滤波处理单元,具体用于:
根据公式:
对关键切片图像进行处理,得到各含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像gn(x,y);
其中,g(x,y)表示关键切片图像,g(x,y)的大小为X×Y像素,X和Y为偶数;gn(x,y)表示g(x,y)中的第n个含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像;Hn(u,v)表示第n个含有幕帘噪声的灰度范围区域对应的二维带阻滤波器。
另外,所述图像替换单元,具体用于:
根据平滑和灰度分割法,对关键切片图像进行平滑处理,然后利用灰度分割将各滤波器对应区域的位置提取出来;
在平滑后的灰度图像gs(x,y)中,获得二维带阻滤波器Hn(u,v)对应的含有幕帘噪声的灰度范围区域的灰度值区间:t1<gs(x,y)<t2;t1和t2为灰度阈值;
根据公式:将原图中二维带阻滤波器Hn(u,v)对应的幕帘噪声的灰度范围区域替换为Hn(u,v)对应的滤波后的图像gn(x,y)的相同区域,形成替换后的灰度图像g′n(x,y);
将各滤波器对应的灰度范围区域替换后,形成关键切片图像对应的降噪处理图像。
一种计算机可读存储介质,其上存储有计算机程序,该程序被处理器执行时实现以下步骤:
步骤101、根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致;
步骤102、在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著;
步骤103、根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值;
步骤104、对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数;
步骤105、用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像;
步骤106、将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像;
步骤107、将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数;
步骤108、将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像;
步骤109、将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;
重复步骤101至109,直至步骤110、获得FIB-SEM图像序列所有切片图像的降噪处理图像。
一种计算机设备,包括存储器、处理器及存储在存储上并可在处理器上运行的计算机程序,所述处理器执行所述程序时实现以下步骤:
步骤101、根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致;
步骤102、在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著;
步骤103、根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值;
步骤104、对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数;
步骤105、用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像;
步骤106、将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像;
步骤107、将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数;
步骤108、将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像;
步骤109、将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;
重复步骤101至109,直至步骤110、获得FIB-SEM图像序列所有切片图像的降噪处理图像。
本发明实施例提供的一种去除FIB-SEM图像幕帘噪声的方法及装置,方法包括:将待处理的FIB-SEM图像序列按主要矿物组分变化分为若干组,在每个切片图像组中,选择一幅关键切片图像;根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数;用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像;将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像;将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数;将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像;将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;重复上述过程直至获得FIB-SEM图像序列所有切片图像的降噪处理图像。可见,本发明可以实现一种针对岩石的FIB-SEM图像的有效、快速、节省人力的幕帘噪声去除方法,并避免产生滤波伪影,获得细节清晰、结构真实的图像。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动性的前提下,还可以根据这些附图获得其他的附图。
图1为现有技术中一组连续的FIB-SEM图像序列示意图;
图2为具有明显幕帘噪声的富有机质的页岩样品的FIB-SEM图像切片示意图;
图3为本发明实施例提供的一种去除FIB-SEM图像幕帘噪声的方法的流程图一;
图4为本发明实施例提供的一种去除FIB-SEM图像幕帘噪声的方法的流程图二;
图5为压制幕帘噪声的二维带阻滤波器设计示意图;
图6为在FIB-SEM图像序列中选取的第420幅切片的示意图;
图7为二维傅立叶变换的振幅谱示意图;
图8为第420幅切片经二维带阻滤波后的示意图;
图9为第420幅切片滤波前后的差值示意图;
图10为提取的需替换的区域,即原图中没有幕帘噪声的高亮度黄铁矿区域的示意图;
图11为第420幅切片的最终降噪结果示意图;
图12为第20、220、620幅切片的去幕帘噪声前后的示意图;
图13为本发明实施例提供的一种去除FIB-SEM图像幕帘噪声的装置的结构示意图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
如图3所示,本发明实施例提供一种去除FIB-SEM图像幕帘噪声的方法,包括:
步骤101、根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致。
步骤102、在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著。
步骤103、根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值。
步骤104、对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数。
步骤105、用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像。
步骤106、将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像。
步骤107、将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数。
步骤108、将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像。
另外,若测试灰度范围区域和带阻滤波器的参数的降噪处理效果不合格,则需要重新观察、划分该组切片,修正该组降噪参数的适用范围。
另外,如果灰度范围区域和带阻滤波器的参数对该组切片的首切片降噪效果合格,而对尾切片的降噪效果不合格,则说明其适用范围是从该组首切片到关键切片图像之后的某一幅。这时可从尾切片向前查看,在适当位置(例如矿物组分有较大变化处)将原切片组一分为二,并在此试用灰度范围区域和带阻滤波器的参数。如果合格,则此处为该组切片新的尾切片,从这里向前到首切片为该灰度范围区域和带阻滤波器的参数的适用范围;从这里向后到原来的尾切片成为新的切片组,需重复步骤101至107确定其灰度范围区域和带阻滤波器的参数。如果不合格,则继续向前寻找,直到关键切片图像处。同理,如果对首切片降噪效果不合格,而对尾切片降噪效果合格,或首、尾切片都不合格,处理方法与此类似,此处不再一一赘述。
步骤109、将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像。
重复步骤101至109,直至步骤110、获得FIB-SEM图像序列所有切片图像的降噪处理图像。
本发明实施例提供的一种去除FIB-SEM图像幕帘噪声的方法可以实现针对岩石的FIB-SEM图像的有效、快速、节省人力的幕帘噪声去除方法,并避免产生滤波伪影,获得细节清晰、结构真实的图像。
为了使本领域的技术人员更好的了解本发明,下面列举一个更为详细的实施例,如图4所示,本发明实施例提供一种去除FIB-SEM图像幕帘噪声的方法,包括:
步骤201、根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致。
步骤202、在每个切片图像组中,选择一幅关键切片图像。
其中,所述关键切片图像存在明显幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著。
步骤203、根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值。
此处对于关键切片图像,可以定性分析幕帘噪声强度与局部灰度的关系。一般情况下,不同矿物往往具有不同硬度,因此,不同灰度区域的幕帘噪声可能具有不同特征。进而可以将关键切片图像划分为若干个灰度范围区域,同一灰度范围区域内的幕帘噪声特征相似,而不同灰度范围区域之间的幕帘噪声特征则差别较大。
例如,可将关键切片图像划分为有幕帘噪声的软矿物区(有机质以及黏土等)和基本无幕帘噪声的硬矿物区。
步骤204、对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据,根据公式:将部分图像数据f(x,y)进行二维傅里叶变换,确定部分图像数据f(x,y)对应的频谱F(u,v)。
其中,f(x,y)表示进行截取处理后获得的含有幕帘噪声的灰度范围区域的部分图像数据,f(x,y)的大小为M×N像素,x和y分别表示横坐标像素和纵坐标像素;i为虚数单位,u为频谱F(u,v)的横坐标空间频率,u=–M/2,–M/2+1,…,0,1,…,M/2–1;v为频谱F(u,v)的纵坐标空间频率,v=–N/2,–N/2+1,…,0,1,…,N/2–1;M和N为偶数。
由上式计算得到的频谱F(u,v)中,横、纵坐标轴为空间频率,其只与空间采样率相关,纵坐标空间频率在FIB-SEM图像中可表示每像素的宽度,通常为几至几十纳米。同一幅二维图像在横、纵方向上通常具有相同的像素宽度Δl,故空间频率范围是–1/2Δl~1/2Δl。频域中,空间频率的间隔Δk与采样点数有关。横、纵坐标的Δk分别为1/MΔl和1/NΔl。若M≠N,则横、纵坐标的Δk不等,这在设计和使用滤波器时要特别注意。
步骤205、根据频谱F(u,v)得到振幅谱|F(u,v)|,并在振幅谱|F(u,v)|中检测得到幕帘噪声所在的空间频率范围。
步骤206、在幕帘噪声所在的空间频率范围内压制振幅,得到二维带阻滤波器H(u,v)。
应注意的是,基于傅立叶变换的滤波器可能在突变边界处引入伪波。为了减少这类干扰,通带和阻带之间一般需要设置过渡带。
其中,所述二维带阻滤波器H(u,v)的通带和阻带之间设置有过渡带,过渡带函数可以为余弦型,即为其中,L为过渡带宽度;j=1,2,…,L,为相对于L的位置;该过渡带函数在下降区间时为在上升区间时为当然,也可以应用其他过渡带函数,此处不再一一列举。如图5所示,是压制幕帘噪声的二维带阻滤波器设计示意图,包括通带(频域大部分)、阻带(幕帘噪声的空间频率范围)和过渡带三类区域。图5中,中心点为频域原点(空间频谱为0);阴影区是宽度为L的过渡带;阴影区包围的区域是需要压制的幕帘噪声的频谱范围;阴影区之外是通带。
步骤207、将幕帘噪声所在的空间频率范围和过渡带宽度L作为二维带阻滤波器H(u,v)的参数。
通过二维带阻滤波器H(u,v)的参数,后续可计算各含有幕帘噪声的灰度范围区域对应的滤波后的图像,进而可检测、改进二维带阻滤波器的效果。最终,应无明显的幕帘噪声或伪波。
步骤208、根据公式:
对关键切片图像进行处理,得到各含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像gn(x,y)。
其中,g(x,y)表示关键切片图像,g(x,y)的大小为X×Y像素,X和Y为偶数;gn(x,y)表示g(x,y)中的第n个含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像;Hn(u,v)表示第n个含有幕帘噪声的灰度范围区域对应的二维带阻滤波器。
应注意的是,不同滤波器对应的图像区域具有地质意义,针对不同的矿物大类,它们对应的灰度范围不重叠。还需注意,由于采用g(x,y)的局部图像来设计滤波器,g(x,y)频谱的空间频率间隔小于局部图像频谱的间隔,需对滤波器插值,以便应用于g(x,y)。
步骤209、根据平滑和灰度分割法,对关键切片图像进行平滑处理,然后利用灰度分割将各滤波器对应区域的位置提取出来。
例如可采用均值滤波法,以减少图像噪声对边缘形态的影响。
步骤210、在平滑后的灰度图像gs(x,y)中,获得二维带阻滤波器Hn(u,v)对应的含有幕帘噪声的灰度范围区域的灰度值区间:t1<gs(x,y)<t2
其中,t1和t2为灰度阈值。
步骤211、根据公式:将原图中二维带阻滤波器Hn(u,v)对应的幕帘噪声的灰度范围区域替换为Hn(u,v)对应的滤波后的图像gn(x,y)的相同区域,形成替换后的灰度图像g′n(x,y)。
步骤212、将各滤波器对应的灰度范围区域替换后,形成关键切片图像对应的降噪处理图像。
步骤213、将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数。
具体的,确定降噪处理效果是否合格采用的是分析差值图像(滤波前后图像之差)。理想的差值图像特征为:在低振幅的随机噪声背景上,叠加条纹状的幕帘噪声,不同位置的幕帘噪声可能有不同空间频率和振幅,但应大致平行(这是由FIB-SEM的成像方式决定的);不包含其它具有相关性的信号。可通过自相关函数判断数组的随机性。差值图像上,与幕帘噪声方向、空间频率不一致的波纹状信号一般是带阻滤波产生的伪波,通常是过渡带设计不当造成的。表现在频域,理想的差值图像的振幅谱只在原幕帘噪声的空间频率范围内有明显信号;其它区域是随机噪声背景,其平均振幅一般不超过幕帘噪声平均振幅的十分之一。
步骤214、将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像。
步骤215、将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像。
重复步骤201至215,直至步骤216、获得FIB-SEM图像序列所有切片图像的降噪处理图像。
为了使本领域的技术人员更加了解上述步骤201至步骤216,下面给出一个具体实例:
以图1所示的FIB-SEM图像序列为例,具体说明本发明的实施方式。
该序列共有650幅连续切片,切片间距10纳米,每幅切片2048×1768像素,每像素5纳米。因此,所测区域大小为6.5μm×10.24μm×8.84μm。
样品主要包含以下类别的区域:
(1)高亮度的黄铁矿区域,基本无幕帘噪声。
(2)中等亮度的基质矿物区域,部分可见幕帘噪声。
(3)深灰色的有机质区域,可见明显的幕帘噪声。
(4)接近黑色的孤立小区域为孔隙。
选取第420幅切片为关键切片图像(如图6所示,灰度值23~244,平均值129,下部的线框表示用于分析幕帘噪声频谱特性的区域)。选用依据是:
(1)存在有明显的幕帘噪声。
(2)除幕帘噪声外,其它噪声不明显。
(3)包含上述全部4类区域。
(4)第420幅切片较为靠近FIB-SEM图像序列的中部。
定性分析认为该切片可分为2个灰度范围区域:无需滤波的高亮区和需要去除幕帘噪声的浅灰以下区域。在该切片中截取一部分(图6下部的线框内的部分),用于分析幕帘噪声的频谱特征,其二维傅立叶变换的振幅谱如图7所示,黑框内即为幕帘噪声的主要空间频率范围,横坐标和纵坐标分别表示X轴空间频率和Y轴空间频率。由此设计二维带阻滤波器(可以见上述图5所示),增加过渡带以减少滤波不当造成的伪波。
滤波效果如图8所示(灰度值25~248,平均值129),滤波前(图6)、滤波后(图8)的差值如图9所示。从步骤213下的判断依据可知,竖条纹状的幕帘噪声被较好地压制了,伪波较少,对样品原有特征的还原程度好。例如,图6线框中的基质矿物有较弱的密度(灰度)变化,但因受幕帘噪声影响,很难分辨。经滤波,图8中对应位置的灰度变化更清晰了,揭示出该位置的矿物正在经历某种地球化学变化(可能是矿物蚀变)。
尽管该滤波器是针对中等亮度的基质矿物(图6线框区域)设计的,它对深灰色的有机质区域的降噪效果也令人满意。因此,在这个算例中,不必再为有机质另外设计滤波器。也就是说,上述步骤中,只需设计和使用一个滤波器。
另外,需要注意的是,图8中箭头所指处。对比图6相同位置,原无幕帘噪声的高亮度区域出现了明暗变化,致使该处亮度、清晰度下降。相应地在图9中,箭头所指处出现了最大的差值,但此处原无幕帘噪声,不应出现差值。因此,按步骤209至步骤212的方法,可先将原图中高亮度区的位置提取出来,即采取10点均值滤波后,灰度值≥160的区域(图10所示,提取的需替换的区域,即原图中没有幕帘噪声的高亮度黄铁矿区域(仅显示包含所提取区域的局部。))。再将原图中该区域替换到滤波后的图像中。得到第420幅切片的最终降噪结果(图11所示),改善了图8中高亮区域的问题。
由此,可得到关于以上关键切片(第420幅)的参数组:
(1)需设计的二维带阻滤波器数量,在本例中为1个。
(2)需分区的数量,本例为2个,原有幕帘噪声,因而需要滤波的(用同一个滤波器),和原无幕帘噪声,故不需滤波的。
(3)二维带阻滤波器参数,包括在空间频域中的阻带范围、过渡带函数类型和参数,本例采用余弦型过渡带,阻带尽可能只包含幕帘噪声的主要空间频率范围。
(4)为了划分各滤波(或不滤波)区域,图像需要做的预处理及参数,以及提取各区域所采用的灰度阈值,本例为10点均值滤波后,提取灰度值≥160的区域为不需滤波区域。
如前所述,本套FIB-SEM测试结果包含650幅连续切片。由于测试区域小,岩性变化不大。依照步骤214,另选数幅切片进行试验,发现前述参数组适用于全部切片。因此,不必选用其它关键切片另行设计参数组。
将这些参数应用于全套切片,获得了完整的处理结果。图12展示第20、220、620幅,效果较好(左列三幅图像为去幕帘噪声前的图像,右列三幅图像为去幕帘噪声后的图像)。第420幅的详细处理过程和结果见上述步骤201至步骤214。
本发明实施例提供的一种去除FIB-SEM图像幕帘噪声的方法可以实现针对岩石的FIB-SEM图像的有效、快速、节省人力的幕帘噪声去除方法,并避免产生滤波伪影,获得细节清晰、结构真实的图像。
对应上述图3和图4所示的方法实施例,如图13所示,本发明实施例还提供一种去除FIB-SEM图像幕帘噪声的装置,包括:
切片图像分组单元31,用于根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致。
关键切片图像选择单元32,用于在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著。
区域划分单元33,用于根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值。
截取处理单元34,用于对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数。
滤波处理单元35,用于用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像。
图像替换单元36,用于将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像。
区域及参数记录单元37,用于将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数。
降噪处理效果测试单元38,用于将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像。
降噪处理图像获得单元39,用于将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像。
重复切片图像分组单元31至降噪处理图像获得单元39的过程,直至获得FIB-SEM图像序列所有切片图像的降噪处理图像。
另外,所述截取处理单元34,具体用于:
对每一个含有幕帘噪声的灰度范围区域进行截取处理,以获得每一个含有幕帘噪声的灰度范围区域的部分图像数据。
根据公式:将部分图像数据f(x,y)进行二维傅里叶变换,确定部分图像数据f(x,y)对应的频谱F(u,v);其中,f(x,y)表示进行截取处理后获得的含有幕帘噪声的灰度范围区域的部分图像数据,f(x,y)的大小为M×N像素,x和y分别表示横坐标像素和纵坐标像素;i为虚数单位,u为频谱F(u,v)的横坐标空间频率,u=–M/2,–M/2+1,…,0,1,…,M/2–1;v为频谱F(u,v)的纵坐标空间频率,v=–N/2,–N/2+1,…,0,1,…,N/2–1;M和N为偶数。
根据频谱F(u,v)得到振幅谱|F(u,v)|,并在振幅谱|F(u,v)|中检测得到幕帘噪声所在的空间频率范围。
在幕帘噪声所在的空间频率范围内压制振幅,得到二维带阻滤波器H(u,v)。
所述二维带阻滤波器H(u,v)的通带和阻带之间设置有过渡带,过渡带函数为其中,L为过渡带宽度;j=1,2,…,L,为相对于L的位置;该过渡带函数在下降区间时为在上升区间时为
将幕帘噪声所在的空间频率范围和过渡带宽度L作为二维带阻滤波器H(u,v)的参数。
另外,所述滤波处理单元35,具体用于:
根据公式:
对关键切片图像进行处理,得到各含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像gn(x,y);
其中,g(x,y)表示关键切片图像,g(x,y)的大小为X×Y像素,X和Y为偶数;gn(x,y)表示g(x,y)中的第n个含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像;Hn(u,v)表示第n个含有幕帘噪声的灰度范围区域对应的二维带阻滤波器。
另外,所述图像替换单元36,具体用于:
根据平滑和灰度分割法,对关键切片图像进行平滑处理,然后利用灰度分割将各滤波器对应区域的位置提取出来。
在平滑后的灰度图像gs(x,y)中,获得二维带阻滤波器Hn(u,v)对应的含有幕帘噪声的灰度范围区域的灰度值区间:t1<gs(x,y)<t2;t1和t2为灰度阈值。
根据公式:将原图中二维带阻滤波器Hn(u,v)对应的幕帘噪声的灰度范围区域替换为Hn(u,v)对应的滤波后的图像gn(x,y)的相同区域,形成替换后的灰度图像g′n(x,y)。
将各滤波器对应的灰度范围区域替换后,形成关键切片图像对应的降噪处理图像。
另外,本发明实施例还提供一种计算机可读存储介质,其上存储有计算机程序,该程序被处理器执行时实现以下步骤:
步骤101、根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致;
步骤102、在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著;
步骤103、根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值;
步骤104、对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数;
步骤105、用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像;
步骤106、将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像;
步骤107、将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数;
步骤108、将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像;
步骤109、将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;
重复步骤101至109,直至步骤110、获得FIB-SEM图像序列所有切片图像的降噪处理图像。
另外,本发明实施例还提供一种计算机设备,包括存储器、处理器及存储在存储上并可在处理器上运行的计算机程序,所述处理器执行所述程序时实现以下步骤:
步骤101、根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致;
步骤102、在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著;
步骤103、根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值;
步骤104、对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数;
步骤105、用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像;
步骤106、将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像;
步骤107、将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数;
步骤108、将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像;
步骤109、将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;
重复步骤101至109,直至步骤110、获得FIB-SEM图像序列所有切片图像的降噪处理图像。
本发明实施例提供的一种去除FIB-SEM图像幕帘噪声的装置,可以实现针对岩石的FIB-SEM图像的有效、快速、节省人力的幕帘噪声去除方法,并避免产生滤波伪影,获得细节清晰、结构真实的图像。
本领域内的技术人员应明白,本发明的实施例可提供为方法、系统、或计算机程序产品。因此,本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器,使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述,以上实施例的说明只是用于帮助理解本发明的方法及其核心思想;同时,对于本领域的一般技术人员,依据本发明的思想,在具体实施方式及应用范围上均会有改变之处,综上所述,本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种去除FIB-SEM图像幕帘噪声的方法,其特征在于,包括:
步骤101、根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致;
步骤102、在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著;
步骤103、根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值;
步骤104、对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数;
步骤105、用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像;
步骤106、将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像;
步骤107、将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数;
步骤108、将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像;
步骤109、将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;
重复步骤101至109,直至步骤110、获得FIB-SEM图像序列所有切片图像的降噪处理图像。
2.根据权利要求1所述的去除FIB-SEM图像幕帘噪声的方法,其特征在于,所述对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数,包括:
根据公式:将部分图像数据f(x,y)进行二维傅里叶变换,确定部分图像数据f(x,y)对应的频谱F(u,v);其中,f(x,y)表示进行截取处理后获得的含有幕帘噪声的灰度范围区域的部分图像数据,f(x,y)的大小为M×N像素,x和y分别表示横坐标像素和纵坐标像素;i为虚数单位,u为频谱F(u,v)的横坐标空间频率,u=–M/2,–M/2+1,…,0,1,…,M/2–1;v为频谱F(u,v)的纵坐标空间频率,v=–N/2,–N/2+1,…,0,1,…,N/2–1;M和N为偶数;
根据频谱F(u,v)得到振幅谱|F(u,v)|,并在振幅谱|F(u,v)|中检测得到幕帘噪声所在的空间频率范围;
在幕帘噪声所在的空间频率范围内压制振幅,得到二维带阻滤波器H(u,v);
所述二维带阻滤波器H(u,v)的通带和阻带之间设置有过渡带,过渡带函数为其中,L为过渡带宽度;j=1,2,…,L,为相对于L的位置;该过渡带函数在下降区间时为在上升区间时为
将幕帘噪声所在的空间频率范围和过渡带宽度L作为二维带阻滤波器H(u,v)的参数。
3.根据权利要求2所述的去除FIB-SEM图像幕帘噪声的方法,其特征在于,所述用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像,包括:
根据公式:
对关键切片图像进行处理,得到各含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像gn(x,y);
其中,g(x,y)表示关键切片图像,g(x,y)的大小为X×Y像素,X和Y为偶数;gn(x,y)表示g(x,y)中的第n个含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像;Hn(u,v)表示第n个含有幕帘噪声的灰度范围区域对应的二维带阻滤波器。
4.根据权利要求3所述的去除FIB-SEM图像幕帘噪声的方法,其特征在于,所述将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像,包括:
根据平滑和灰度分割法,对关键切片图像进行平滑处理,然后利用灰度分割将各滤波器对应区域的位置提取出来;
在平滑后的灰度图像gs(x,y)中,获得二维带阻滤波器Hn(u,v)对应的含有幕帘噪声的灰度范围区域的灰度值区间:t1<gs(x,y)<t2;t1和t2为灰度阈值;
根据公式:将原图中二维带阻滤波器Hn(u,v)对应的幕帘噪声的灰度范围区域替换为Hn(u,v)对应的滤波后的图像gn(x,y)的相同区域,形成替换后的灰度图像g′n(x,y);
将各滤波器对应的灰度范围区域替换后,形成关键切片图像对应的降噪处理图像。
5.一种去除FIB-SEM图像幕帘噪声的装置,其特征在于,包括:
切片图像分组单元,用于根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致;
关键切片图像选择单元,用于在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著;
区域划分单元,用于根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值;
截取处理单元,用于对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数;
滤波处理单元,用于用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像;
图像替换单元,用于将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像;
区域及参数记录单元,用于将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数;
降噪处理效果测试单元,用于将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像;
降噪处理图像获得单元,用于将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;
重复切片图像分组单元至降噪处理图像获得单元的过程,直至获得FIB-SEM图像序列所有切片图像的降噪处理图像。
6.根据权利要求5所述的去除FIB-SEM图像幕帘噪声的装置,其特征在于,所述截取处理单元,具体用于:
对每一个含有幕帘噪声的灰度范围区域进行截取处理,以获得每一个含有幕帘噪声的灰度范围区域的部分图像数据;
根据公式:将部分图像数据f(x,y)进行二维傅里叶变换,确定部分图像数据f(x,y)对应的频谱F(u,v);其中,f(x,y)表示进行截取处理后获得的含有幕帘噪声的灰度范围区域的部分图像数据,f(x,y)的大小为M×N像素,x和y分别表示横坐标像素和纵坐标像素;i为虚数单位,u为频谱F(u,v)的横坐标空间频率,u=–M/2,–M/2+1,…,0,1,…,M/2–1;v为频谱F(u,v)的纵坐标空间频率,v=–N/2,–N/2+1,…,0,1,…,N/2–1;M和N为偶数;
根据频谱F(u,v)得到振幅谱|F(u,v)|,并在振幅谱|F(u,v)|中检测得到幕帘噪声所在的空间频率范围;
在幕帘噪声所在的空间频率范围内压制振幅,得到二维带阻滤波器H(u,v);
所述二维带阻滤波器H(u,v)的通带和阻带之间设置有过渡带,过渡带函数为其中,L为过渡带宽度;j=1,2,…,L,为相对于L的位置;该过渡带函数在下降区间时为在上升区间时为
将幕帘噪声所在的空间频率范围和过渡带宽度L作为二维带阻滤波器H(u,v)的参数。
7.根据权利要求6所述的去除FIB-SEM图像幕帘噪声的装置,其特征在于,所述滤波处理单元,具体用于:
根据公式:
对关键切片图像进行处理,得到各含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像gn(x,y);
其中,g(x,y)表示关键切片图像,g(x,y)的大小为X×Y像素,X和Y为偶数;gn(x,y)表示g(x,y)中的第n个含有幕帘噪声的灰度范围区域所对应滤波器的滤波后的图像;Hn(u,v)表示第n个含有幕帘噪声的灰度范围区域对应的二维带阻滤波器。
8.根据权利要求7所述的去除FIB-SEM图像幕帘噪声的装置,其特征在于,所述图像替换单元,具体用于:
根据平滑和灰度分割法,对关键切片图像进行平滑处理,然后利用灰度分割将各滤波器对应区域的位置提取出来;
在平滑后的灰度图像gs(x,y)中,获得二维带阻滤波器Hn(u,v)对应的含有幕帘噪声的灰度范围区域的灰度值区间:t1<gs(x,y)<t2;t1和t2为灰度阈值;
根据公式:将原图中二维带阻滤波器Hn(u,v)对应的幕帘噪声的灰度范围区域替换为Hn(u,v)对应的滤波后的图像gn(x,y)的相同区域,形成替换后的灰度图像g′n(x,y);
将各滤波器对应的灰度范围区域替换后,形成关键切片图像对应的降噪处理图像。
9.一种计算机可读存储介质,其上存储有计算机程序,其特征在于,该程序被处理器执行时实现以下步骤:
步骤101、根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致;
步骤102、在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著;
步骤103、根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值;
步骤104、对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数;
步骤105、用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像;
步骤106、将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像;
步骤107、将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数;
步骤108、将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像;
步骤109、将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;
重复步骤101至109,直至步骤110、获得FIB-SEM图像序列所有切片图像的降噪处理图像。
10.一种计算机设备,包括存储器、处理器及存储在存储上并可在处理器上运行的计算机程序,其特征在于,所述处理器执行所述程序时实现以下步骤:
步骤101、根据FIB-SEM图像序列中各切片图像所表现的主要矿物组分变化,将各切片图像分为一至多组,每个切片图像组的主要矿物组分一致;
步骤102、在每个切片图像组中,选择一幅关键切片图像;所述关键切片图像存在幕帘噪声,且包括其所在的切片图像组的主要矿物组分,且除幕帘噪声之外的其它噪声不显著;
步骤103、根据关键切片图像上各位置的幕帘噪声特征值,将所述关键切片图像划分为多个灰度范围区域;其中,同一灰度范围区域内各位置的幕帘噪声特征值小于等于预先设置的特征阈值,任意两个不同灰度范围区域内位置的幕帘噪声特征值大于所述特征阈值;
步骤104、对于每一个含有幕帘噪声的灰度范围区域,截取部分图像数据进行二维傅里叶变换,确定该部分图像数据包含的幕帘噪声的频谱特征,并形成每一个含有幕帘噪声的灰度范围区域对应的带阻滤波器,确定各带阻滤波器的参数;
步骤105、用各带阻滤波器处理关键切片图像,得到各带阻滤波器对应的滤波后的图像;
步骤106、将关键切片图像中的含有幕帘噪声的灰度范围区域替换成各自对应的滤波后的图像,形成关键切片图像对应的降噪处理图像;
步骤107、将降噪处理图像与所述关键切片图像进行比较,确定降噪处理效果是否合格;并在确定降噪处理效果合格时,记录灰度范围区域和带阻滤波器的参数;
步骤108、将灰度范围区域和带阻滤波器的参数应用于关键切片图像所在的切片图像组的首尾两个切片图像,测试灰度范围区域和带阻滤波器的参数的降噪处理效果是否合格;如果合格,则认为该灰度范围区域和带阻滤波器的参数适用于该切片图像组中的所有切片图像;
步骤109、将该灰度范围区域和带阻滤波器的参数应用于该切片图像组,获得该切片图像组的降噪处理图像;
重复步骤101至109,直至步骤110、获得FIB-SEM图像序列所有切片图像的降噪处理图像。
CN201811013158.7A 2018-08-31 2018-08-31 一种去除fib-sem图像幕帘噪声的方法及装置 Pending CN109255764A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811013158.7A CN109255764A (zh) 2018-08-31 2018-08-31 一种去除fib-sem图像幕帘噪声的方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811013158.7A CN109255764A (zh) 2018-08-31 2018-08-31 一种去除fib-sem图像幕帘噪声的方法及装置

Publications (1)

Publication Number Publication Date
CN109255764A true CN109255764A (zh) 2019-01-22

Family

ID=65050038

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811013158.7A Pending CN109255764A (zh) 2018-08-31 2018-08-31 一种去除fib-sem图像幕帘噪声的方法及装置

Country Status (1)

Country Link
CN (1) CN109255764A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113570505A (zh) * 2021-09-24 2021-10-29 中国石油大学(华东) 一种页岩三维超分辨率数字岩心分级重构方法及系统
CN113592741A (zh) * 2021-08-04 2021-11-02 西北工业大学 一种数字图像处理方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010177217A (ja) * 2010-05-21 2010-08-12 Sii Nanotechnology Inc Fib−sem複合装置
CN106057620A (zh) * 2015-04-15 2016-10-26 Fei 公司 在带电粒子显微镜中执行层析成像的方法
CN107067379A (zh) * 2017-03-16 2017-08-18 中国科学院地质与地球物理研究所 基于三维fib‑sem图像的页岩孔隙定量表征方法
CN108318491A (zh) * 2017-12-04 2018-07-24 华南理工大学 一种基于频谱曲率分析的织物缺陷检测方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2010177217A (ja) * 2010-05-21 2010-08-12 Sii Nanotechnology Inc Fib−sem複合装置
CN106057620A (zh) * 2015-04-15 2016-10-26 Fei 公司 在带电粒子显微镜中执行层析成像的方法
CN107067379A (zh) * 2017-03-16 2017-08-18 中国科学院地质与地球物理研究所 基于三维fib‑sem图像的页岩孔隙定量表征方法
CN108318491A (zh) * 2017-12-04 2018-07-24 华南理工大学 一种基于频谱曲率分析的织物缺陷检测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
S.LIU 等: "A fast curtain‐removal method for 3D FIB‐SEM images of heterogeneous minerals", 《JOURNAL OF MICROSCOPY-WILEY ONLINE LIBRARY》 *
SHUAI LIU 等: "A Staged Filtering Approach to Kill Curtain Noise in FIB-SEM Images", 《RESEARCHGATE》 *
章鲁 等: "《医学图像处理与分析》", 31 August 2006, 上海科学技术出版社 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113592741A (zh) * 2021-08-04 2021-11-02 西北工业大学 一种数字图像处理方法
CN113570505A (zh) * 2021-09-24 2021-10-29 中国石油大学(华东) 一种页岩三维超分辨率数字岩心分级重构方法及系统
CN113570505B (zh) * 2021-09-24 2022-01-04 中国石油大学(华东) 一种页岩三维超分辨率数字岩心分级重构方法及系统

Similar Documents

Publication Publication Date Title
Chandra et al. A critical review on pore to continuum scale imaging techniques for enhanced shale gas recovery
CN106597539B (zh) 针对黄土塬地区的曲波域Radon变换噪声压制方法
CN108549106B (zh) 混叠噪声压制方法及装置
CN104932010B (zh) 一种基于近道镶边稀疏Radon变换的绕射波分离方法
CN104007469A (zh) 一种基于曲波变换的弱地震信号重构方法
de Oliveira Lyrio et al. Efficient automatic denoising of gravity gradiometry data
CN107144879A (zh) 一种基于自适应滤波与小波变换结合的地震波降噪方法
CN104849757B (zh) 消除地震信号中随机噪声系统及方法
Wang et al. Robust vector median filtering with a structure-adaptive implementation
Liu et al. A fast curtain‐removal method for 3D FIB‐SEM images of heterogeneous minerals
CN109255764A (zh) 一种去除fib-sem图像幕帘噪声的方法及装置
AU2014342986B2 (en) Method for post-stack noise mitigation
CN108169795A (zh) 基于随机采样的数据规则化方法
CN113158830A (zh) 一种剩余重力异常场分离方法
Belila et al. Petrophysical characterization of coquinas from Morro do Chaves Formation (Sergipe-Alagoas Basin) by X-ray computed tomography
CN109917459A (zh) 一种压制地震噪声的方法、装置及系统
Li et al. Fracture extraction from FMI based on multiscale mathematical morphology
Ushizima et al. Material science image analysis using Quant-CT in ImageJ
Wen et al. Highlighting display of geologic bodies based on directivity filtering
CN109655916A (zh) 用于分离地震数据中有效波与多次波的方法及系统
Liu et al. A staged filtering approach to kill curtain noise in FIB-SEM images
CN113589384A (zh) 基于信号随偏移距变化特征的叠前道集保幅去噪方法
Ferreira et al. The importance of outcrop reservoir characterization in oil-industry facies modelling workflows-a case study from the Middle Jurassic of the Maciço Calcário Estremenho, Portugal
Cui et al. Fault feature enhancement in seismic data based on steerable pyramid tensor voting
CN117908092A (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
RJ01 Rejection of invention patent application after publication

Application publication date: 20190122

RJ01 Rejection of invention patent application after publication