基于压缩感知的频域OCT图像稀疏重构

陈明惠，王 帆，张晨曦，李福刚，郑 刚

(1. 上海理工大学 生物医学光学与视光学研究所 教育部现代微创医疗器械及技术工程研究中心 上海介入医疗器械工程技术研究中心，上海 200093；2. 上海奥普生物医药有限公司，上海 201203)

1 引 言

光学相干断层扫描(Optical Coherence Tomography,OCT)成像技术利用宽带光源的低相干特性，以非侵入、无损伤的方式进行断层成像，已广泛应用于医学诊断和研究[1-2]。频域OCT技术因其高灵敏度和高速性渐渐取代传统时域OCT技术，但是拥有这些优点的同时意味着后续数据采集和处理系统的造价更高且负担更重[3]；同时，该成像技术需要对病变部位进行反复多角度多层次的深度扫描以获得大量高分辨率图像，高分辨率必然伴随着高数据量，同时在对眼底部位进行成像时，眼球的抖动会造成运动模糊，所以迫切需要一种能降低硬件系统存储和传输代价且花费时间较少并获得较高图像质量的方式。

2006年，Donoho[4]，E Candes[5-6]及华裔科学家Tao[7]等人提出了一种新的信息采集方式，即压缩感知(Compressive Sensing，CS)理论。CS技术表明，如果信号(图像)本身是稀疏的，或能在某个稀疏集上进行稀疏表示，则通过一定的测量矩阵进行测量，加以选取适当的优化算法，在获得少量测量值的基础上，能对原始信号进行有效地恢复重构[8]。研究高性能的恢复重构算法是CS领域的一个热点，目前优化算法主要分为贪婪算法、凸优化算法和非凸算法三类[9]。研究者们在匹配追踪算法(Matching Pursuit，MP)的基础上获得了更多的贪婪算法，Wei等人[10]提出的SAMP(Sparsity Adaptive Matching Pursuit)重构算法通过更新支持集并逐渐增加稀疏度来逼近原始信号，实现了稀疏度未知的信号重构，但重构精度不够高，重构后SNR最高仅为25 dB左右；Chen等人[11]利用小波变换获得图像的稀疏表示，使用改进的最佳正交匹配追踪算法完成MRI图像的重构，但算法计算效率不高，算法运行时间约为26 s；Ma等人[12]提出一种变步长SAMP算法，在一定程度上提高了重构精度，但引入分块思想后，图像重构出现明显的人为噪声，导致重构质量下降；石昊苏等人[13]采用共轭梯度算法替换OMP算法中的最小二乘法求取估计值进行超声图像的重构，提高了重构质量，但PSNR值仅达到13.75 dB，算法运行时间约为16.6 s。第二类凸优化算法通过将非凸问题转换成凸问题进行求解，Qu等人[14]将基追踪算法和非均匀快速傅立叶变换技术相结合，实现了观测场景的成像重构，提高了成像质量，但算法运行时间约为39 s，算法计算复杂度高，应用范围受限；郑万泽等人[15]基于Contourlet变换，采样改进梯度投影算法(Gradient Projection for Sparse Reconstruction，GPSR)恢复稀疏处理后的系数，实现了二维Lena和Barbara图像的重构，采样率为0.5时重构图像的PSNR值仅为30 dB左右且算法实时性较差。第三类非凸算法通过改变目标函数求解优化问题，该类算法克服了基追踪算法中的计算量大问题，同时可以采用比正交匹配追踪算法少的测量值达到同样的重建效果，主要以FOCUSS算法[16]和稀疏贝叶斯学习[17]为代表。

上述文献中改进的压缩感知重构算法在保证算法的实时性以及图像最终重构质量之间总是存在矛盾，算法运行速度最快约为17 s，PSNR值最高为30 dB。为了在一定程度上实现频域OCT图像重构效率和重构质量之间的平衡，本文以非凸算法中的FOCUSS算法为基础进行改进，改进后算法运行时间仅耗时约2 s，PSNR值均在37 dB以上，具体从两方面进行，一是根据频域OCT图像特点，结合分块思想降低图像的数据量，引入lp范数简化运算过程，提高算法运行速度；二是从提高图像重构质量进行考虑，对图像进行分块会造成严重的块效应，因此在每次迭代求解过程中嵌入各向异性平滑滤波算子，抑制图像块效应，提高图像质量，实验结果表明，压缩感知技术可成功应用于频域OCT图像，且改进的压缩感知重构算法在频域OCT图像重构中能花费较短的时间获得较高的重构质量。

2 算法原理

2.1 压缩感知

传统的信号获取需满足采样频率大于两倍信号最高频率的Nyquist采样定理。而CS利用图像(信号)的稀疏性，以远低于传统采样方法获取的数据量便可完成初始信号的完整恢复，突破了传统采样定理的限制[4-8]，因此本文将压缩感知技术引入应用到频域OCT图像的重构上，下面将介绍CS理论的三个组成部分。

(1)稀疏表示

信号的稀疏表示是应用压缩感知的前提，假设φi是RN空间的一组向量基，构成基矩阵Ψ=[φ1,φ2，...，φN]，则对于任意的频域OCT图像x∈RN均可表示为：

(1)

式中：α是x在Ψ域的稀疏表示，若α中非零元数目比x少得多，则称信号x是稀疏的。常见的稀疏变换基有DCT，DWT(Discrete Wavelet Transform)，DDWT(Discrete Dyadic Wavelet Transform)，Curvelet以及冗余字典基等，在下文仿真实验中选择DCT基获得频域OCT图像的稀疏表示。

(2)观测矩阵

压缩感知理论指出，若采用一个与稀疏变换基Ψ不相关的测量矩阵Φ∈RM×N，对频域OCT图像x进行投影得到观测值y，则可以通过求解优化问题高概率精确地重构出x，表达式为：

y=Φx=ΦΨα，

(2)

式中：A=ΦΨ为传感矩阵，要求测量矩阵Φ和稀疏变换基Ψ不相关，即在较高概率上满足约束等距性条件，常见的观测矩阵分为两类，确定性矩阵和随机矩阵，本文在仿真实验中选择第二类高斯随机矩阵对频域OCT图像进行观测。

(3)重构算法

对式(2)求解可转化为式(3)：

(3)

求解上述l0范数是NP难解问题，需要列出α中所有非零项位置的线性组合才有可能得到精确解，计算过程复杂且不易实现。因此常将该问题转换为其他替代模型进行求解。常用的重构算法有三类，本文利用第三类非凸算法中的改进FOCUSS重构算法进行求解，接下来将详细阐述该改进算法的原理及实现。

2.2 改进重构算法

利用CS可实现对传统二维图像的重构，但此方法是对整幅图像进行重构，观测矩阵所需的存储空间大，重构时间长。为了解决这个问题，本文提出分块FOCUSS重构算法，根据频域OCT图像特点，将图像分块理论和FOCUSS重构算法相结合。具体做法是先将一幅大小为N×N的OCT图像x分成n个B×B大小的图像块，设xi(i=1,2,...,n,n=(N/B)2)是第i个图形块的列向量形式，采用相同的高斯观测矩阵ΦB，则对第i块进行压缩采样可表示为：

yi=ΦB·xi.

(4)

本文中压缩观测矩阵ΦB大小为nB×B2。

传统FOCUSS算法中由于伪逆的存在，往往导致能量分散，无法获得精确的稀疏解[18]。将lp范数作为正则项引入到FOCUSS算法中，通过逐步迭代使解的能量局部化。对每一块分别应用FOCUSS算法，设置p=1，用以简化运算过程，即求解以下优化问题:

(5)

使用拉格朗日方法将上式中约束优化问题转化为非约束优化问题，其中，拉格朗日函数被定义为:

L(xi,λ)=f(xi)+λT(ΦBxi-yi)，

(6)

式中λ=(λ1,λ2,...,λmB)T构成拉格朗日乘数。

(7)

(8)

为了找到目标函数的梯度，将矩阵Ψ分解为其行向量φi得到:

(9)

分解之后，目标函数表示为：

(10)

目标函数相对于的xi的梯度为：

(11)

目标函数相对于元素xj的偏导数为:

(12)

使用因子形式表示目标函数的梯度向量即:

经过半年多的运转，浮选车间运行状况稳定。浮选精煤灰分控制在8%以内，三班工作制，每班精煤压滤机可卸3～4个循环，按每个循环8 t计，则每天可生产浮选精煤80 t左右。每月可多洗精煤约2 500 t，按当时精煤价格1 300元/t测算，每月可增加效益近350万元。经测算，实际运行4个月后，已收回投资成本。

(13)

式中Π(xψ)=diag(|xψ|-1)是大小为n×n的对角矩阵且xψ=Ψxi。

参考式(8)和式(13)，得出:

(14)

参考式(5)，得出:

(15)

综合式(7)，式(14)，式(15)，稳定点满足如下两个表达式:

(16)

(17)

(18)

将式(18)代入式(17)，得出:

(19)

将式(20)代入式(19)，得出:

(20)

定义Aψφ=(ΦBΨ-1)，经过简单的转换，将(20)写成:

(21)

(22)

(23)

通过修改拉格朗日函数来调整松弛方法以找到迭代解，得到如下方程:

(24)

(25)

对所有的小图像块进行相同的处理，再把得到的图像块加以组合就完成了整幅频域OCT图像的重构。本文的算法流程图如图1所示。

图1 改进的分块压缩感知FOCUSS重构算法流程图

3 仿真实验与结果

本文选择两种频域OCT图像进行仿真实验，其中指尖图像来自实验室Thorlabs公司VEGA系列的SS-OCT成像系统，光源中心波长为1 310 nm，穿透深度为12 mm，轴向分辨率为16 μm，扫频速率为100 kHz，系统灵敏度为101 dB。眼底图像来自杜克大学眼科实验数据采集中心[19]，成像仪器是德国海德堡公司的Spectralis SD-OCT系统，光源中心波长为870 nm，穿透深度为1.8 mm。纵向分辨率和横向分辨率分别为3.8 μm和14 μm。实验室采集到的SS-OCT指尖图像为455×512 pixel的RGB图像，图像细节较少，整体较平坦。杜克大学数据库中的SD-OCT眼底图像为496×512 pixel的灰度图像，图像细节较丰富。为了方便对其进行分块处理，均将图像大小调整为512×512 pixel。

3.1 评价指标

为了验证本文算法能否有效重构频域OCT图像，除了直接从人眼视觉效果评价重构图像质量，同时采用峰值信噪比(Peak Signal to Noise Ratio，PSNR)，用RPSNR表示，结构相似性(Structural Similarity Index，SSIM)，用MSSIM表示，取值范围0-1，以及重构时间t定量评价算法的重构性能。RPSNR值越大，表明与完整采样图像相比，重构图像与之越接近，重构效果越好，重构精度越高；MSSIM值越大，表明重构图像保留了越多的原始完整采样图像中的信息，重构图像质量越好；重构时间t越小，算法的重构效率越高。

3.2 实验结果分析

本文中所有实验在64位Windows10系统下利用软件MATLAB R2018a进行，采样率(Sampling Rate，SR)是重构过程中采样值大小z与原始信号x大小的比值，首先对实验室采集到的指尖OCT图像进行压缩感知重构处理，利用DCT基对原始图像稀疏表示，采用高斯随机矩阵对图像进行观测，并运用5种优化算法对指尖图像进行恢复重构。图2给出SR为0.3，分块大小为8×8时原始FOCUSS重构算法、当前的主流分块重构算法-平滑投影算法(SPL)以及改进FOCUSS重构算法的恢复效果以及局部细节对比图。

据图2可知，当图像采样率SR为0.3，分块大小为8×8时，改进FOCUSS重构算法对OCT指尖图像进行重构后能够获得更好的视觉重构效果，图像细节保持完整，图像边缘处容易分辨，与原图更接近。图2(b)表明直接运行原始FOCUSS重构算法得到的重构图像有一定程度的信息缺失，图像的边缘处不光滑；从图2(c)～图2(e)可以看出3种主流的图像压缩感知分块算法对图像进行重构后，重构图像的细节存在明显的边界效应；通过图2(f)可以得出改进的FOCUSS算法显著改善了图像块效应，获得的重构图像细节清晰，轮廓明显，满足人眼视觉要求，成功将压缩感知技术应用到扫频OCT系统产生的指尖图像上，本文算法重构细节较少，画面平坦的彩色图像效果较好。

为了便于定量比较实验结果，将采样率固定为0.3，重复试验10次，取平均值作为最终仿真结果值，表1给出原始FOCUSS算法、使用3种不同稀疏变换基(DCT，DWT，DDWT)的分块平滑投影重构算法以及改进FOCUSS算法对指尖OCT图像的压缩感知重构客观评价指标值。

图2 采样率为0.3时不同算法重构效果及局部细节图

表1SR为0.3时各算法对指尖图像的重构结果

Tab.1 Reconstruction results of fingertip images by each algorithm when the SR is 0.3

算法PSNR/dBSSIMt/s原始FOCUSS35.340.878 778.65分块SPL-DCT35.010.825 68.64分块SPL-DWT35.660.853 113.79分块SPL-DDWT36.040.865 710.72改进FOCUSS37.710.938 31.89

据表1可知，与原始FOCUSS算法相比，改进FOCUSS算法对图像进行重构后提高了图像质量，PSNR值约高出2 dB，SSIM值提高到0.938 3，并且对图像进行分块后大大降低了算法运行时间，算法运行仅花费1.89 s，提高了算法运行效率，算法实时性强；与同样结合了分块思想的SPL算法相比，改进FOCUSS算法不仅改善了图像的块效应，能获得更高的PSNR值和SSIM值，重构图像质量高，与原始完整采样OCT指尖图像更接近且保留了图像中的大部分信息，而且算法运行花费时间更短，说明该算法计算复杂度更低，运算效率高，更易推广使用。同时，为了说明本文算法的普遍性，另外进行了4组仿真实验，将图像分块大小分别固定为8×8，16×16，32×32，64×64 4种情况，仿真获取SR分别为0.1，0.2，0.3，0.4和0.5时5种算法对指尖图像的重构结果图和客观评价指标值，实验图和数据均表明，当图像分块大小为8×8，16×16，32×32时，采样率大于0.2时，改进的压缩感知重构算法较以往算法能获得重构速度和重构精度的平衡，是具有一定优势的改进算法。

上述仿真采用的指尖图像细节较少，为了进一步说明本文算法的有效性，同时具体分析采样率与OCT图像重构质量的关系，随机选取一张眼底OCT切片，切片图像的细节较丰富，将采样率分别设置为0.1，0.2，0.3，0.4和0.5，图3～图7给出了当采样率为0.1，0.3和0.5时各算法对完整采样眼底OCT图像的重构结果。据五组图可以看出，随着采样率的增加，重构图像的细节、轮廓越来越清晰，视觉效果越来越好，与原始完整采样图像越来越接近。据图3～图7(b)可知当采样率为0.1时，各算法都不能很好地重构原始完整采样图像，图像信息缺失严重，图像边缘位置无法确定，图像视觉效果差；据图3～图7(c)和图7(d)可知，随着采样率增加至0.3以上时，重构后的图像更加清晰，边缘细节保持完整，边缘处较为光滑，图像轮廓基本保持良好，能获得较高的图像质量，容易分辨各视网膜层，与原始完整采样图像更接近，从视觉效果看并无明显差别。所以采用压缩感知技术可以用较少的采样数据量精确重构出谱域OCT系统产生的原始眼底OCT图像，本文算法对细节丰富的灰度图像的重构效果也较好，适用范围较大。

图3 原始FOCUSS算法重构结果

图4 分块SPL-DCT算法重构结果

图5 分块SPL-DWT算法重构结果

图6 分块SPL-DDWT算法重构结果

图7 改进FOCUSS算法重构结果

为了直观反映重构算法的重构性能和采样率的关系，图8给出了不同采样率下各算法对眼底OCT切片图像重构后的客观评价指标值变化曲线图。可以看出随着采样率的增加，各算法对图像的重构效果越来越好，重构质量越来越高，重构时间并没有大幅度变化，重构性能越来越好。据图8(a)可知，在采样率为0.1和0.2时，各算法对图像重构后的PSNR值均在36 dB以下，重构质量不高，但是随着采样率增加至0.3以上时，本文改进的结合了分块思想的FOCUSS重构算法的PSNR值迅速高于其他4种对比算法，图像重构精度高；从图8(b)可以看出改进的FOCUSS算法的SSIM值始终高于其他算法，对图像进行重构后保留了图像中的大部分有用信息，图像重构质量高。图8(c)则可以看出原始FOCUSS重构算法的重构时间均在75 s以上，算法重构效率不高，而改进的算法则将重构时间控制在2 s以下，算法运行时间也远远低于3种SPL算法，优势显著，性能良好，综合来看，改进重构算法能同时获得较高的图像重构质量和重构效率。

本文改进的压缩感知FOCUSS重构算法结合了图像分块思想，为了进一步分析分块大小对频域OCT图像重构质量和重构效率的影响，本文将采样率固定为0.5，将图像块大小分别设置为8×8，16×16，32×32，64×64对眼底图像进行改进压缩感知重构算法的仿真实验，重复10次，取平均值，最终结果如表2所示。

图8 5种算法在不同采样率下的PSNR值、SSIM值以及重构时间

表2 改进FOCUSS算法在不同分块大小下的仿真结果

Tab.2 Simulation results of improved FOCUSS algorithm under different block sizes

PSNR/dBSSIMt/s8×842.470.964 11.9616×1642.290.961 52.1632×3241.840.956 84.6164×6440.990.946 553.45

据表2可知，随着OCT图像分块大小的增加，算法的运行时间加长，重构效率降低，算法实时性变差，这是由于高斯随机观测矩阵对各个块图像进行测量是同时进行的，前一个图像块对后一个图像块没有依赖性，图像分块越大，观测矩阵越大，系统的存储代价越大，算法复杂度越高。从PSNR值和SSIM值的变小可以看出图像重构后图像质量有小幅度降低，这是由于本文处理的OCT图像并不是标准测试图像，自身不可避免含有散斑噪声，并且在算法迭代过程中结合了各向异性平滑滤波算法，图像块尺寸越小，边缘平滑得越好，散斑也得到抑制。综合来看，当图像分块大小为8×8时，重构图像的PSNR值高达42.47 dB，SSIM值高达0.964 1，接近于1，算法运行仅耗时1.96 s。所以综合来看，频域OCT图像分块大小为8×8时，在图像进行重构后可以同时获得较高的重构质量和重构效率。

4 结 论

本文首先介绍了压缩感知的基本框架，通过分析频域OCT图像特点，成功将压缩感知技术引入到频域OCT图像稀疏重构上。通过对细节丰富程度不同的两种频域OCT图像在不同采样率下及不同分块大小的仿真实验得出，本文提出的改进FOCUSS重构算法通过结合分块思想、引入正则项lp范数以及在算法迭代求解过程中嵌入各向异性滤波算子，不仅降低了计算复杂度，同时有效抑制了重构图像的块效应，提高了图像质量，对OCT图像进行压缩感知重构后能获得更好的视觉效果以及更短的计算时间，本文算法适应范围较广。通过进一步仿真得到当图像分块大小为8×8时能同时获得较高的重构质量和重构速度，PSNR值高达42.47 dB，远高于文献14中的30 dB，SSIM值可达到0.964 1，重构时间仅花费1.96 s，比文献12中的运行时间16.6 s要快得多，本文改进的压缩感知重构算法在一定程度上实现了重构效率和重构精确度的平衡，具有一定的优势。