CN113034374A

CN113034374A - 一种oct图像自适应噪声估计方法及系统

Info

Publication number: CN113034374A
Application number: CN201911348675.4A
Authority: CN
Inventors: 秦嘉; 刘华宗
Original assignee: Guangdong Weiren Medical Technology Co ltd
Current assignee: Guangdong Weiren Medical Technology Co ltd
Priority date: 2019-12-24
Filing date: 2019-12-24
Publication date: 2021-06-25

Abstract

本发明公开了一种OCT图像自适应噪声估计方法及系统，根据系统采集图像的信号与噪声特点，结合双变量收缩复小波降噪理论，进行了相关数学模型的建立与实际OCT图像的降噪实验，多帧图像叠加降噪因其良好的降噪效果被广泛应用，但对特殊病例无法采集足够多帧图、病理诊断过程对特征细节要求高等问题无法很好解决；在无法提供足够帧图情况下，依然能提供高质量的成像效果图像，对于需要去除背景噪声前提下最大程度保留特征细节信号的场景，能在保留图像细节特征前提下，有效降低噪声对图像质量的影响。

Description

一种OCT图像自适应噪声估计方法及系统

技术领域

本公开涉及光学相干层析技术、图像处理领域，具体涉及一种OCT图像自适应噪声估计方法及系统。

背景技术

光学相干层析技术(Optical Coherence Tomography，OCT)^[1]于1991年问世，是一种新型的光学成像技术，该技术采用低相干光源作为激励，能够对探测样品进行层析成像。系统的深度方向分辨率由宽带光源的低相干长度决定，可达到微米级，是超声成像的10-100倍。随着近三十年的发展，从时域到频域的光谱域和扫频域，OCT系统的灵敏度和成像速度都有了很大的提升。由于OCT优越成像性能以及临床检查的需要，使得在眼科等临床检测中得到广泛的应用，其中，光谱域OCT(SD-OCT)相对时域拥有更高的成像速度以及更好的灵敏度，对比扫频域则拥有更成熟的计算以及更便宜的系统造价，是商用OCT系统中常用的一种。

虽然OCT系统能提供高分辨率的图像，但是由于系统特性，成像图像会带有许多噪声，造成图像特征模糊，质量下降等影响，其中白噪声影响尤为严重。为消除这类白噪声，常用办法是采集多张同位置图像，然后进行平均降噪处理，此方法能够有效消除白噪声影响。但是这意味着系统要采集更多数据，从而使得系统需要更长的采集时间以及较大的内存。因此，对单帧成像图像进行预处理降噪十分必要。

目前，针对OCT图像单帧降噪的方法主要集中在时域和频域分析上。Rogowska,J等研究人员^[2]将超声图像降噪技术--旋转核变换(Rotating Kernel Transformation，RKT)应用于OCT图像,该方法能在抑制背景噪声的同时，识别图像边界进行特征区域增强。但该算法的噪声抑制效果差，无法有效的区分图像中信号与噪声混叠区域。Alexander，Wong等人^[3]考虑OCT成像图像噪声的贝叶斯先验分布，成功分离出噪声和信号，降低了OCT噪声对图像质量的影响。Farzana,Z等人^[4]将OCT成像图转换到小波变换域分析，使用自适应小波阈值降噪，有效的减少背景白噪声，但同时也降低微小结构细节的可见度。

发明内容

现有的OCT图像处理方法，对图像帧数要求较高，需要较多帧图像才能得到降噪效果良好的图像数据，而且对细节处理不够精细，不能满足诊断需求，尤其是，对一些特殊病例无法采集足够多帧图、而病理诊断过程对特征细节要求又高时，不能满足需求。

为此，本公开基于OCT图像噪声的分布特性，结合双变量收缩复小波降噪理论，进行了相关数学模型的建立与实际OCT图像的降噪实验，提出一种OCT图像自适应噪声估计方法图像，具体包括以下步骤：

S100：采集OCT图像数据；

S200：通过双树复小波变换将图像数据的实际观测信号值转换到小波域；

S300：估算小波域中的各个子带的全局噪声的标准差；

S400：估算小波域中的各个子带的局部噪声的标准差；

S500：计算各个子带中修正后的小波系数；

S600：根据修正后的小波系数重构复原降噪后的图像。

进一步地，在S300中，估算小波域中的各个子带的全局噪声的标准差的方法为：

式中，w_s代表当前子带的小波系数，T为自适应参数，N表示子带中小波系数的总数，median函数是返回给定数值的中值的函数。

进一步地，在S400中，估算小波域中的各个子带的局部噪声的标准差的方法为：

式中,σ_y为领域内小波系数的方差，M代表当前邻域N(k)的大小。

进一步地，在S500中，计算各个子带中修正后的小波系数的方法为：

通过结合双变量收缩贝叶斯阈值模型可以得到联合收缩公式：

式中，

表示为处在第j层小波尺度s位置修正后的小波系数，

表示为处在第j层小波尺度s位置的小波系数,

表示去噪后的小波系数,σ_n为全局估计的噪声值，σ为位置s邻域估计的噪声值，g()+为收缩函数，定义为：

本发明还提供了一种OCT图像自适应噪声估计系统，所述系统包括：存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序运行在以下系统的单元中：

图像采集单元，用于采集OCT图像数据；

小波域转换单元，用于通过双树复小波变换将图像数据的实际观测信号值转换到小波域；

全局噪声估算单元，用于估算小波域中的各个子带的全局噪声的标准差；

局部噪声估算单元，用于估算小波域中的各个子带的局部噪声的标准差；

小波系数修正单元，用于计算各个子带中修正后的小波系数；

图像重构单元，用于根据修正后的小波系数重构复原降噪后的图像。

本公开的有益效果为：本发明公开的OCT图像自适应噪声估计方法，在无法提供足够帧图情况下，依然能提供高质量的成像效果图像，且能在去除背景噪声前提下最大程度保留特征细节信号的场景，图像从而可以适应更广泛的应用场景。

附图说明

通过对结合附图所示出的实施方式进行详细说明，本公开的上述以及其他特征将更加明显，本公开附图中相同的参考标号表示相同或相似的元素，显而易见地，下面描述中的附图仅仅是本公开的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，在附图中：

图1所示为本公开的一种OCT图像自适应噪声估计方法的流程图；

图2为实验系统采集到的人眼后节OCT成像图像和进行模型降噪后的人眼后节OCT图像对比图；

图3为降噪前后含有眼部血管区域图像和降噪前后脉络膜后层含有微血管区域图像感兴趣区域对比图；

图4为OCT原始图像的截取感兴趣区域、自适应中值滤波后截取的感兴趣区域、双树复小波降噪后截取的感兴趣区域、本公开方法降噪后截取的感兴趣区域对比图；

图5为使用100帧未做处理图像所得到的多帧平均图、使用8帧图像进行本公开方法先验处理所得到的多帧平均对比图；

图6所示为本公开实施方式的一种OCT图像自适应噪声估计系统结构图。

具体实施方式

以下将结合实施例和附图对本公开的构思、具体结构及产生的技术效果进行清楚、完整的描述，以充分地理解本公开的目的、方案和效果。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。

如图1所示为根据本公开的一种OCT图像自适应噪声估计方法的流程图，下面结合图1来阐述根据本公开的实施方式的方法。

本公开提出一种OCT图像自适应噪声估计方法，具体包括以下步骤：

S100：采集OCT图像数据；

S300：估算小波域中的各个子带的全局噪声的标准差；

S400：估算小波域中的各个子带的局部噪声的标准差；

S500：计算各个子带中修正后的小波系数；

S600：根据修正后的小波系数重构复原降噪后的图像。

其中，采集的OCT图像是被方差为

的高斯白噪声污染图像，其数学公式可以描述为：

m_(s)＝a_(s)+n_(s) (1)

式中s为图像素的空间位置坐标，m为图像的实际观测信号值，a为理论分析得到的信号值，n为噪声的估计值。

采用双树复小波变换将m转换到小波域w中，结合双变量收缩贝叶斯阈值模型^[5]可以得到联合收缩公式：

式中

表示为处在第j层小波尺度s位置的小波系数,

表示去噪后的小波系数,σ_n为全局估计的噪声值，σ为位置s邻域估计的噪声值，双变量收缩贝叶斯阈值模型见参考文献[5]。

其中，g(X)₊为收缩函数其定义为

由公式(2)可知，得到去噪后的小波系数

需要计算σ_n和σ。其中σ_n传统估计方法为文献^[6]所提出的，将当前子带的所有小波系数排序，并取其中值除以一个鲁棒系数得到，如公式(3)。

median函数是返回给定数值的中值的函数。

该噪声估计使用在大多数图像中都具鲁棒性，但对于OCT图像，其实际效果不是很理想，因此，本申请融合OCT图像噪声特性，提出的自适应噪声方差估计如公式(4)：

式中，w_s代表当前子带的小波系数，T为自适应参数，N表示子带中小波系数的总数。对于σ的估计引入邻域范围内的均值进行计算，因为模型建立在公式(1)中，故σ可以通过

得出，如公式(5)：

实施例

将本公开方法应用于一组由SD-OCT系统采集到的人眼后节图像。从多名志愿者眼睛部位采集二维成像的OCT眼底断层结构图，如图2的(a)所示。所采集部位摄于视网膜视盘区域，图像具有清晰的分层结构和良好的形态学特征，图中可看到主干血管的阴影部分，脉络膜后层微毛细血管。但因为OCT成像系统引入的背景噪声，直接影响了图2的(a)中各细节部分的展现。

在图2中的(a)中实验系统采集到的人眼后节OCT成像图像。(b)进行上述方法降噪后的人眼后节OCT图像。输入图像分辨率为1024*1000，依据本公开所提出方法，先对图2的(a)使用双树复小波变换(Dual Tree-Complex Wavelet Transform，DT-CWT)，其分解尺度Level＝5。根据文献^[7]介绍复小波分解所具有的方向选择性与多分辨率特性，复小波域中的图像在每一层分辨率上会分出6个方向子带。对每一个方向子带应用本公开所建立的数学模型，然后进行小波重构恢复降噪后图像。最终效果如图2中的(b)所示，通过对比可知，本公开算法能在保留OCT图像形态学特征的前提下去除背景含有的高斯白噪声。

针对细节的比较，可以从图3中的a、b区域对比中发现，在去除噪声的同时，图4中b中各部分层次结构清晰，血管部分并未因为降噪算法而模糊。图3的c与d区域对比可以看到，降噪前脉络膜血管信号与噪声混叠区域在降噪算法处理后被显示出来。

图3中感兴趣区域比较，图3中的a、b分别表示降噪前后含有眼部血管区域图像，图3中c、d分别表示降噪前后脉络膜后层含有微血管区域图像。

性能评价

为了评价本公开所构建的方法与同类相关算法在降低OCT图像中高斯白噪声的性能，选用了其他高性能的小波域滤波方法和基于时域的图像降噪算法，如T.Loupas等人^[8]提出的自适应中值滤波，I.W.Selesnick等人^[9]提出的双树复小波滤波，使用相关研究文献中提出的参数，高效的实现以上测试方法。对于本公开所提出的方法，经过超过100组数据的验证，选用邻域大小为17*17，双树复小波分解选用qshift_a与sym_a作为基，自适应参数T＝0.45，能更好的降低图像受到高斯白噪声的影响。控制单一变量，以上所有算法，都将使用python(3.6版本)语言进行实现，并运行在一台搭载了(i5-4570cpu 3.3Ghz 8G内存)的计算器上进行测试。为了合理的估计各算法的性能，选取了如图4中显示了4个感兴趣区域，聚焦于血管的形态学特征和更加深层次的信号分析。

图4的(A)为OCT原始图像的截取感兴趣区域，图4的(B)为自适应中值滤波后截取的感兴趣区域，图4的(C)为双树复小波降噪后截取的感兴趣区域图4的(D)为本公开方法降噪后截取的感兴趣区域，对处理后的图像进行主观视觉评价，发现自适应中值滤波(Adaptive Median Filtering，AMF)^[8]保留了图像细节，但滤除噪声效果不佳，无法有效的分离出背景噪声和特征信号。双树复小波降噪(DT-CWT)^[9]能保留大部分细节的同时滤除背景噪声，但在内部层次结构分层和脉络膜后端信号成像上模糊。对比本公开方法，在保留细节的同时滤除背景噪声，同时保留了完好的层次结构。

同时，为了定量比较各检测方法的性能，计算每种方法在人眼OCT图像上的峰值信噪比(SNR)，对比度信噪比(CNR)，等效视数(ENL)。使用的评价指标与文献^[3,10,11]相同，具体定义如下：

SNR＝10log₁₀[max(A²)/σ²] (6)

在SNR公式中，A代表信号区域灰度值，σ为选定背景区域噪声的标准差大小，SNR越大，代表图像品质越好。在ENL公式与CNR公式中，h为所选取感兴趣区域编号，H为感兴趣区域数，μ为感兴趣区域的平均值，σ为该区域的标准差大小，b表示背景噪声区域编号，其中ENL值越大说明图中所选区域平滑度越高，CNR值越大说明图像所选区域对比度越高。

表1给出了原始图像测量数据和各滤波器数据的定量比较，其中可以看到本公开所提出的算法在SNR值呈现明显优势，即对微小细节的处理提升明显，相比于原始图像提升为77.35％。对比于DT-CWT算法，CNR和ENL略微逊色，SNR呈现略微提高。实验结果表明，本公开方法能在保留图像细节特征前提下，有效降低噪声对图像质量的影响，信噪比提升约为77.35％。

表1.原始图像测量数据和各滤波器数据的定量比较

同时，还通过与现有眼科OCT临床降噪图像的多帧图像平均算法^[12]对比，评价了本公开算法的数据处理对图像帧数的要求。

现有的多帧图像平均算法核心思想为对图像进行多帧信息互补，采取的是用时间换图像质量的思想，该方法针对不同情况图像的降噪，具有一定鲁棒性。但对于人眼OCT图像采集有一个不可避免的因素，需要考虑人的眼动所造成帧图间特征位置不对应问题，故而需要引入配准机制^[13]，匹配算法的性能将会直接决定最终输出图像的成像质量，帧图越多，因配准算法导致的误差将会越大。

对比效果如表2所示，表2所示为平均降噪和本公开方法，相同效果下所需要帧数对比。

表2.平均降噪和本公开方法，相同效果下所需要帧数对比

算法	SNR	CNR	ENL	帧数
					多帧平均降噪	38.306	0.0413	3.3287	100
本公开方法	43.814	0.0410	3.3032	8

图5中的(a)为使用100帧未做处理图像所得到的多帧平均图。图5中的(b)为使用8帧图像进行本公开方法先验处理所得到的多帧平均图。

从主观视觉上分析，图5的(a)视觉显示更好，图像更加平滑，但在100帧的叠加过程中，丢失了部分细节特征，如图中红色圈区域，其中血管部分的图像已经被稀释掉。图5的(b)，因为算法过程中考虑了邻域平均，导致降噪后图像视觉效果下降，但得益于复小波变换对高频细节的保护，绝大多数细节特征都被保留下来。而在临床眼科OCT中，眼底层析成像的局部细节是临床诊断疾病的一个重要基准。

而且，图像在临床检测中，多帧图像平均算法需要较多的帧图数，而对于某些特殊病例，无法接受长时间的优质图像采集。若考虑本公开算法，在进行多帧平均前对每帧图像按照本公开算法进行降噪预先处理，达到相同降噪效果下，所需要帧图数量大幅减少。

综上，本公开提出了一种针对OCT图像去除背景噪声的算法，其对单帧图像经过DT-CWT变换后的小波系数进行增益处理，增益方法为考虑上下层级之间的贝叶斯先验关系。不同于OCT处理中常见的多帧平均方法，该方法使用单帧图像信号作为输入。

定量评价结果显示，本公开方法去除OCT图像背景噪声的效果优于自适应中值滤波或传统的小波降噪方法。同时在与多帧降噪算法比较中，只用8帧图像就能达到100帧叠加降噪图像相同的CNR和ENL。在眼科OCT临床检测中，保证某些疾病患者无法提供足够帧图情况下，依然能提供高质量的成像效果图像，对于需要去除背景噪声前提下最大程度保留特征细节信号的场景也有应用价值。

另外，本公开的实施例也提供的一种OCT图像自适应噪声估计系统，如图6所示为本公开的一种OCT图像自适应噪声估计系统结构图，该实施例的一种OCT图像自适应噪声估计系统包括：处理器、存储器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述一种OCT图像自适应噪声估计系统实施例中的步骤。

所述系统包括：存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序运行在以下系统的单元中：

图像采集单元，用于采集OCT图像数据；

所述一种OCT图像自适应噪声估计系统可以运行于桌上型计算机、笔记本、掌上电脑及云端服务器等计算设备中。所述一种OCT图像自适应噪声估计系统可运行的系统可包括，但不仅限于，处理器、存储器。本领域技术人员可以理解，所述例子仅仅是一种OCT图像自适应噪声估计系统的示例，并不构成对一种OCT图像自适应噪声估计系统的限定，可以包括比例子更多或更少的部件，或者组合某些部件，或者不同的部件，例如所述一种OCT图像自适应噪声估计系统还可以包括输入输出设备、网络接入设备、总线等。所称处理器可以是中央处理单元(Central Processing Unit，CPU)，还可以是其他通用处理器、数字信号处理器(Digital Signal Processor，DSP)、专用集成电路(Application Specific IntegratedCircuit，ASIC)、现场可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等，所述处理器是所述一种OCT图像自适应噪声估计系统运行系统的控制中心，利用各种接口和线路连接整个一种OCT图像自适应噪声估计系统可运行系统的各个部分。

所述存储器可用于存储所述计算机程序和/或模块，所述处理器通过运行或执行存储在所述存储器内的计算机程序和/或模块，以及调用存储在存储器内的数据，实现所述一种OCT图像自适应噪声估计系统的各种功能。所述存储器可主要包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需的应用程序(比如声音播放功能、图像播放功能等)等；存储数据区可存储根据手机的使用所创建的数据(比如音频数据、电话本等)等。此外，存储器可以包括高速随机存取存储器，还可以包括非易失性存储器，例如硬盘、内存、插接式硬盘，智能存储卡(Smart Media Card,SMC)，安全数字(SecureDigital,SD)卡，闪存卡(Flash Card)、至少一个磁盘存储器件、闪存器件、或其他易失性固态存储器件。

参考文献:

[1]Huang D,Swanson E,Lin C,et al.Optical coherence tomography[J].Science,1991,254(5035):1178-1181.

[2]Rogowska J,Brezinski M E.Evaluation of the Adaptive SpeckleSuppression Filter for Coronary Optical Coherence Tomography Imaging[J].IEEETransactions on Medical Imaging,2000,19(12):1261-1266.

[3]Wong A,Mishra A,Bizheva K,et al.General Bayesian estimation forspeckle noise reduction in optical coherence tomography retinal imagery[J].Opt.Express,2010,18:8338-8352.

[4]Aki F a Z,Ang Y a W,Iu X U L.Noise adaptive wavelet thresholdingfor speckle noise removal in optical coherence tomography[J].BiomedicalOptics Express,2017,8(5):2720-2731.

[5]Sendur L,Selesnick I W.Bivariate shrinkage with localvarianceestimation[J].IEEE Signal Processing Letters,2002,9(12):438-441.

[6]L.D,Donoho,Johnstone M,et al.Ideal Spatial Adaptation by WaveletShrinkage[J].Oxford University Press on behalf of Biometrika Trust,1994,81(3):425-455.

[7]Selesnick I W,Baraniuk R G,Kingsbury N C.The dual-tree complexwavelet transform[J].IEEE Signal Processing Magazine,2005,22(6):123-151.

[8]Loupas T,Mcdicken W N,Allan P L.An adaptive weighted median filterfor speckle suppression in medical ultrasonic images[J].IEEE Transactions onCircuits and Systems,1989,36(1):129-135.

[9]Selesnick I W,Baraniuk R G,Kingsbury N G.The Dual-Tree ComplexWavelet Transform[J].IEEE Signal Processing Magazine,2005,22:no.6.

[10]Puvanathasan P,Bizheva K.Interval type-II fuzzy anisotropicdiffusion algorithm for speckle noise reduction in optical coherencetomography images[J].Optics Express,2009,17(2)(2):733.

[11]Puvanathasan P,Bizheva K.Speckle noise reduction algorithm foroptical coherence tomography based on interval type II fuzzy set[J].Opt.Express,2007,15:15747-15758.

[12]Sander B,Larsen M,Thrane L,et al.Enhanced optical coherencetomography imaging by multiple scan averaging[J].British Journal ofOphthalmology,2005,28(2):207-212.

[13]Miri M S,Abràmoff M D,Kwon Y H,et al.Multimodal registration ofSD-OCT volumes and fundus photographs using histograms of oriented gradients[J].Biomed Opt Express,2016,7(12):5252-5267.

尽管本公开的描述已经相当详尽且特别对几个所述实施例进行了描述，但其并非旨在局限于任何这些细节或实施例或任何特殊实施例，而是应当将其视作是通过参考所附权利要求考虑到现有技术为这些权利要求提供广义的可能性解释，从而有效地涵盖本公开的预定范围。此外，上文以发明人可预见的实施例对本公开进行描述，其目的是为了提供有用的描述，而那些目前尚未预见的对本公开的非实质性改动仍可代表本公开的等效改动。