基于低秩矩阵填充的背景荧光噪声抑制方法

王晓东; 耿国华; 易黄建; 何雪磊; 贺小伟

doi:doi:10.3788/AOS201838.1017003

光学学报, 2018, 38 (10): 1017003, 网络出版: 2019-05-09

基于低秩矩阵填充的背景荧光噪声抑制方法下载： 967次

Low-Rank-Matrix-Completion-Based Method for Suppressing Background Fluorescence

论文大纲

王晓东 ^**耿国华 ^*易黄建何雪磊贺小伟

作者单位

西北大学信息科学与技术学院, 陕西西安 710127

AI 词云图 AI一句话精读 AI短摘要

注：本部分内容由 AI 自动生成，请您知悉。

摘要

在面向精准医疗的分子影像领域,荧光分子断层成像(FMT)是当前的研究热点之一。由于FMT逆问题严重的病态性,背景荧光噪声会对重建结果产生严重的负面影响。在深入研究基于有限元的FMT重建方法的基础上,提出利用低秩矩阵填充技术克服背景荧光的方法。该方法将不同激发节点形成的外表面观测组成一个有元素缺失的观测矩阵,利用低秩矩阵填充算法恢复该矩阵的缺失元素,同时抑制观测矩阵含有的背景荧光噪声。利用去噪后的观测矩阵建立了新的FMT逆问题模型,并利用其对荧光目标进行重建。单荧光和双荧光目标重建实验表明:基于去噪后FMT逆问题模型的重建结果获得了显著改善。

Abstract

Fluorescence molecular tomography (FMT) is a hot research topic in molecular imaging applied to precision medicine. Due to the seriously ill posed FMT inverse problem, background fluorescence noise often degrades FMT reconstruction results greatly. After analyzing FMT reconstruction methods based on the finite element method, we propose a method for suppressing background fluorescence using low rank matrix completion technology. The proposed method builds an incomplete boundary observation matrix which columns correspond to different excitation sources. Then, a low rank matrix completion algorithm is employed to complete the matrix and suppress the background fluorescence. At last, a new FMT inverse problem is formed by the denoised boundary observation matrix, with which the fluorescent targets are reconstructed. The new FMT inverse problem is applied to reconstruct single and double fluorescent targets, numerical experiments illustrate that the reconstruction results are improved greatly.

1 引言

荧光分子断层成像(FMT)在分子影像领域是一个重要的研究分支。与其他分子断层成像技术相比,FMT具有副作用小、成本低、精度高、稳定性强等特点,在肿瘤的早期检测、新药物研发等领域具有巨大的应用价值。FMT通过光源重建定位荧光分子探针,从而达到检测早期肿瘤的目的。而背景荧光(BF)是影响FMT光源重建效果的重要因素。这些背景荧光通常有两个来源^[1],第一个来源是那些分布在非目标区域的荧光探针,第二个来源是生物体自身产生的自体荧光。由于FMT逆问题的病态性,这些背景荧光通常对光源重建精度影响很大。因此,研究减小或消除这些背景荧光的方法就显得非常有价值。

很多研究者致力于减小或消除这些背景荧光,所使用的方法大致可以分为3类。

第一类方法是估计单幅视图的背景荧光能量,并将其从外表面观测中减去^[2-4]。Gao等^[2]和Soubret等^[3]利用物理模型估计背景荧光的能量,Ale等^[4]使用Gao等^[2]和Soubret等^[3]给出的物理模型,并为将要减去的背景荧光项增加一个比例因子,这个因子是依据当前外表面观测结果按照一定公式计算获得的。

第二类是采用谱解混的方法^[5-6],这类方法主要用于多光谱荧光图像的去噪,与FMT有紧密的关系,其中的方法值得借鉴。Qin等^[5]利用非负矩阵分解的方法将观测到的多光谱荧光图像分解成不同的端元谱矩阵和丰度分数矩阵的乘积,背景荧光可以被认为是含在某一个端元谱中。Xu等^[6]采用多变量曲线的方法来寻找多光谱荧光图像中的不同谱成分。当然,也有文献采用物理方法测量并形成端元谱矩阵^[7],该方法只需用最小二乘方法计算丰度分数矩阵即可。尽管计算简单,但Qin等^[5]指出,Plaza等^[7]使用的矩阵分解方法效果并不理想。也有些学者将成分分离的方法用于FMT光源重建后的处理^[8],如Pu等^[8]采用独立成分分析的方法将重建出的非目标区域的荧光能量去除。

第三类方法是利用重建算法的稳健性克服背景荧光噪声^[1]。Zhang等^[1]采用L₂-范数正则方法构建FMT逆问题,将前一次逆问题获得的重建结果用3D离散余弦滤波,然后用L₁-范数正则方法抑制目标函数的能量,重复这一过程若干次后,算法停止。Zhang等^[1]采用的方法是专门为克服背景荧光噪声而设计的算法,该算法的重建质量很高。此外,一些基于稀疏正则项的FMT逆问题也具有一定的抑制背景荧光噪声的作用,如基于L₁-范数正则的逆问题^[9]、基于L_p-范数(0<p<1)正则的逆问题^[10]等。

本文的目标是利用低秩矩阵填充方法抑制由分布在非目标区域的荧光探针产生的背景荧光噪声。FMT的光源重建是根据多激发光源点产生的多角度观测协同完成的。每一个激发光源产生的外表面荧光分布可以变形成一个列向量,不同激发光源产生的列向量排在一起可以形成一个外表面观测矩阵。由于采集设备的限制,只能观测到与激发光源相反方向的外表面,可视区域的范围一般为120°或150°。因此,外表面观测矩阵的大部分元素都是缺失的。通过推导可以证明外表面观测矩阵具有低秩属性,因此可以通过矩阵低秩填充的方法恢复,在恢复这些缺失元素的同时,背景荧光噪声将被大幅降低。利用去噪后的外表面观测重建的FMT光源质量将会获得极大改善。

从上一段的描述可以看出,本文提出的方法是利用矩阵这一数学工具,结合生物体外表面荧光分布的先验信息,通过优化方法来达到去除或抑制背景荧光噪声的目的。在本质上讲,本方法应该属于前面提到的第二类去噪方法。当然,本方法与前面提到的方法相比^[1-10],也有非常多的不同之处,概括起来主要为:

首先,在去噪原理方面,本研究证明外表面观测矩阵具有低秩属性,并依据这一属性进行背景荧光去噪。本方法的去噪原理与基于物理模型和端元谱的去除背景荧光方法不同。

其次,在去噪策略方面,本研究采用多视角外表面观测联合背景荧光去噪策略。本方法的去噪测量与单视角外表面观测图像去噪及稳健重建算法去噪有所不同。

再次,在去噪工具方面,本研究创新性地将矩阵低秩填充技术这一数学工具用于背景荧光去噪。本方法去噪所用数学工具与非负矩阵分解、独立成分分析、3D离散余弦滤波等技术不同。

本研究通过仿真实验对所提方法加以验证。实验结果表明,对于不同的噪声强度和不同的光源数量,光源的重建效果都获得了很大提高。

2 原理与方法

2.1 光辐射传输模型及其简化方程

在生物断层成像领域,通常使用扩散方程(DE)来近似模拟光在生物组织中的传输过程^[11]。而FMT有激发和发射两个过程,因此能完整描述这两个光传输过程的方程为

\begin{matrix} \begin{matrix} \nabla [D_{x} (r) \nabla Φ_{x} (r)] - μ_{ax} (r) Φ_{x} (r) = - s (r), (1) \\ \begin{matrix} \nabla [D_{m} (r) \nabla Φ_{m} (r)] - μ_{am} (r) Φ_{m} (r) = - Φ_{x} (r) η (r) μ_{af} (r), (2) \end{matrix} \end{matrix} \end{matrix}

式中:r为3D位置向量;μ_a为介质对光的吸收系数(下角标x表示激发过程,m表示发射过程);D为介质对光的扩散系数;s(r)为激发光源的分布函数;Φ为光子流密度函数;η为荧光产额的转换效率比例;μ_af为荧光团的吸收系数。(1)式和(2)式被统称为耦合扩散方程(CDEs),其中D的计算公式为

\begin{matrix} D = \frac{1}{3 [μ_{a} + (1 - g) μ_{s}]}, (3) \end{matrix}

式中:μ_s为介质对光的散射系数;g为各向异性参数。D、μ_a、μ_s都是r的函数。

直接求解 CDEs并不容易,因此常使用数值的方法求CDEs的近似解。其中,有限元方法是最常用的数值方法^[12-13],该类方法用很小的四面体网格将生物体划分成不同的单元,每一个单元中光子流密度函数用线性函数做近似处理,再结合Robin边界条件^[14-15]就可以将CDEs转换成两个线性方程组:

\begin{matrix} \begin{matrix} K_{x} φ_{x} = s_{x}, (4) \\ K_{m} φ_{m} = K_{F} φ_{x} x \Leftrightarrow φ_{m} = K_{m}^{- 1} K_{F} diag (φ_{x}) x, (5) \end{matrix} \end{matrix}

式中:K为有限元系统矩阵;φ为列向量,φ中的每个元素都是对应有限元节点的光子流密度;s为列向量,其反映激发光源在有限元节点上的分布;K_F为对φ_x有限元进行离散获得的系统矩阵;diag(φ_x)为以φ_x的分量为对角线元素的对角矩阵;x为一个列向量,它指示荧光团的位置,若某个节点在荧光团内部,则x与之对应的分量值为ημ_af,反之则x与之对应的分量值为0。在文献[ 15]中,K_Fdiag(φ_x)用矩阵F表示。对于FMT逆问题,x是待求解的未知量。有限元方法将对微分方程求解的复杂问题转化为对线性方程组的求解问题,大大降低了构建FMT逆问题的难度。

FMT逆问题的构建将在下面两小节详细介绍。

2.2 外表面观测矩阵的低秩属性

(4)式和(5)式分别是模拟单个激发光源s_x产生的激发光和发射光在生物体内的传输过程。在真实实验中,通常需要多个激发光源产生的外表面观测来共同定位生物体中的荧光团。设多个激发光源分别是s_x,1,…,s_x,_l,…,s_x,_L,其中l=1,2,…,L,它们对应的光传输近似方程分别为

\begin{matrix} \{\begin{matrix} K_{x} φ_{x, l} = s_{x, l} \\ φ_{m, l} = K_{m}^{- 1} K_{F} diag (φ_{x, l}) x \end{matrix}, (6) \end{matrix}

(6)式中的K_x、K_m、K_F、x都已知,且是与激发光源s_x,_l无关的量;φ_x,_l可以通过公式φ_x,_l= $\begin{matrix} K_{x}^{- 1} \end{matrix}$ s_x,_l计算获得。FMT光源重建就是已知φ_m,_l,然后利用(6)式的第二个方程计算荧光团分布向量x。然而,φ_m,_l是无法被全部观测到的,只有位于外表面的节点对应的φ_m,_l的分量才有可能被观测到。设这些外表面的节点对应的指标集为I₀,令 $\begin{matrix} {\bar{φ}}_{m, l} \end{matrix}$ = $\begin{matrix} (φ_{m, l})_{I_{0}} \end{matrix}$ ,B= $\begin{matrix} (K_{m}^{- 1} K_{F})_{I_{0} :} \end{matrix}$ ,其中 $\begin{matrix} (φ_{m, l})_{I_{0}} \end{matrix}$ 表示由指标集I₀对应的φ_m,_l的分量形成的子向量, $\begin{matrix} (K_{m}^{- 1} K_{F})_{I_{0} :} \end{matrix}$ 表示由指标集I₀对应的矩阵 $\begin{matrix} K_{m}^{- 1} \end{matrix}$ K_F的行组成的子矩阵。于是,可以得到对FMT光源重建有意义的方程组:

\begin{matrix} {\bar{φ}}_{m, l} = Bdiag (φ_{x, l}) x, l = 1,2, \dots, L 。 (7) \end{matrix}

将(7)式中的 $\begin{matrix} {\bar{φ}}_{m, l} \end{matrix}$ 按l的顺序排成一个矩阵,则有:

\begin{matrix} \begin{matrix} Ψ = [{\bar{φ}}_{m, 1} \dots {\bar{φ}}_{m, l} \dots {\bar{φ}}_{m, L}] = B [diag (φ_{x, 1}) x \dots diag (φ_{x, l}) x \dots diag (φ_{x, L}) x] 。 (8) \end{matrix} \end{matrix}

下面证明外表面观测矩阵Ψ具有低秩属性。

首先给出两个假设,这两个假设在后面的证明中将被用到。假设1:真实发射光光源(荧光团)近似于点光源,激发光在每个真实光源上的分布近似为常数。假设2:在真实发射光光源处的转换效率η远大于在伪发射光光源处的转换效率。

假设2显然是合理的,因为在目标区域(如肿瘤部位)荧光分子探针的密度会非常大,因此更多的激发光被转换为发射光。而在非目标区域,尽管可能存在荧光分子探针,但是其密度一定是非常低的,因此转换效率也必然很低。对于假设1,现实中可能有部分荧光团体积较大而不能近似成为点光源,但是可以通过拆分的形式将其等效为几个点光源,这不会影响后面的推导结果。

为了更直观地展示,这里使用一个2D仿真模型(图1)辅助该证明,其结论很容易被推广到一般的3D模型或真实老鼠数据。该模型含有2个真实荧光目标(图1红色方框含有的节点)和5个等效激发点(图1绿色小圆圈对应的节点),同时该仿体还设置了4个伪光源(图1蓝色椭圆含有的节点)作为背景荧光的来源。为了叙述方便,仿体中的真实荧光目标分别被编号为1号发射光源和2号发射光源,伪光源分别被编号为3~6号,激发点分别被编号为7~11号。

图 1. 2D仿真模型

Fig. 1. 2D simulation model

下载图片查看所有图片

设1号光源节点对应的x的子向量为x¹= $\begin{matrix} (x_{1} x_{2} \dots x_{n_{1}})^{T} \end{matrix}$ ,2号光源节点对应的x的子向量为x²= $\begin{matrix} (x_{n_{1} + 1} x_{n_{1} + 2} \dots x_{n_{2}})^{T} \end{matrix}$ ,3~6号伪光源节点对应的x的子向量分别为x³= $\begin{matrix} (x_{n_{2} + 1} x_{n_{2} + 2} \dots x_{n_{3}})^{T} \end{matrix}$ ,…, x⁶= $\begin{matrix} (x_{n_{5} + 1} x_{n_{5} + 2} \dots x_{n_{6}})^{T} \end{matrix}$ ,其他非光源节点对应的x的子向量为x⁷= $\begin{matrix} (x_{n_{6} + 1} x_{n_{6} + 2} \dots x_{n})^{T} \end{matrix}$ =0,这里0表示元素均为0的列向量。不失一般性,可以令x= $\begin{matrix} [(x^{1})^{T} (x^{2})^{T} \dots (x^{7})^{T}]^{T} \end{matrix}$ 。对x进行一个平凡的拆分,令 $\begin{matrix} {\tilde{x}}^{1} \end{matrix}$ = $\begin{matrix} [(x^{1})^{T} 0^{T} \dots 0^{T}]^{T} \end{matrix}$ , $\begin{matrix} {\tilde{x}}^{2} \end{matrix}$ = $\begin{matrix} [0^{T} (x^{2})^{T} \dots 0^{T}]^{T} \end{matrix}$ ,…, $\begin{matrix} {\tilde{x}}^{7} \end{matrix}$ = $\begin{matrix} [0^{T} \dots 0^{T} (x^{7})^{T}]^{T} \end{matrix}$ 。显然有:

\begin{matrix} x = {\tilde{x}}^{1} + {\tilde{x}}^{2} + \dots + {\tilde{x}}^{7} 。 (9) \end{matrix}

对于第l个激发光源,利用(7)式可得:

\begin{matrix} \begin{matrix} {\bar{φ}}_{m, l} = {\tilde{φ}}_{l}^{1} + {\tilde{φ}}_{l}^{2} + \dots + {\tilde{φ}}_{l}^{7} = B [diag (φ_{x, l}) {\tilde{x}}^{1} + diag (φ_{x, l}) {\tilde{x}}^{2} + \dots + diag (φ_{x, l}) {\tilde{x}}^{7}], (10) \end{matrix} \end{matrix}

式中: $\begin{matrix} {\tilde{φ}}_{l}^{i} \end{matrix}$ =Bdiag(φ_x,_l) $\begin{matrix} {\tilde{x}}^{i} \end{matrix}$ ,i=1,2,…,7。根据假设1可知,对角矩阵diag(φ_x,_l)与 $\begin{matrix} {\tilde{x}}^{i} \end{matrix}$ (i=1,2)非零位置对应的对角线元素为常数,设这个常数为 $\begin{matrix} p_{l}^{i} \end{matrix}$ ,于是有:

\begin{matrix} {\tilde{φ}}_{l}^{i} = p_{l}^{i} B {\tilde{x}}^{i}, i = 1,2, l = 1,2, \dots, L 。 (11) \end{matrix}

对于所有L个激发点,将(10)式代入(8)式,就可以得到外表面观测矩阵:

\begin{matrix} \begin{matrix} Ψ = [\begin{matrix} {\bar{φ}}_{m, 1} & \dots & {\bar{φ}}_{m, l} & \dots & {\bar{φ}}_{m, L} \end{matrix}] = [{\tilde{φ}}_{1}^{1} + {\tilde{φ}}_{1}^{2} + \dots + {\tilde{φ}}_{1}^{7} \dots {\tilde{φ}}_{l}^{1} + {\tilde{φ}}_{l}^{2} + \dots + {\tilde{φ}}_{l}^{7} \dots {\tilde{φ}}_{L}^{1} + {\tilde{φ}}_{L}^{2} + \dots + {\tilde{φ}}_{L}^{7}] = \\ [\begin{matrix} {\tilde{φ}}_{1}^{1} & \dots & {\tilde{φ}}_{l}^{1} & \dots & {\tilde{φ}}_{L}^{1} \end{matrix}] + [\begin{matrix} {\tilde{φ}}_{1}^{2} & \dots & {\tilde{φ}}_{l}^{2} & \dots & {\tilde{φ}}_{L}^{2} \end{matrix}] + \dots + [\begin{matrix} {\tilde{φ}}_{1}^{7} & \dots & {\tilde{φ}}_{l}^{7} & \dots & {\tilde{φ}}_{L}^{7} \end{matrix}] = \\ [\begin{matrix} {\tilde{φ}}_{1}^{1} & \dots & {\tilde{φ}}_{l}^{1} & \dots & {\tilde{φ}}_{L}^{1} \end{matrix}] + [\begin{matrix} {\tilde{φ}}_{1}^{2} & \dots & {\tilde{φ}}_{l}^{2} & \dots & {\tilde{φ}}_{L}^{2} \end{matrix}] + E, (12) \end{matrix} \end{matrix}

式中:E=[ $\begin{matrix} {\tilde{φ}}_{1}^{3} \end{matrix}$ … $\begin{matrix} {\tilde{φ}}_{l}^{3} \end{matrix}$ … $\begin{matrix} {\tilde{φ}}_{L}^{3} \end{matrix}$ ]+…+[ $\begin{matrix} {\tilde{φ}}_{1}^{7} \end{matrix}$ … $\begin{matrix} {\tilde{φ}}_{l}^{7} \end{matrix}$ … $\begin{matrix} {\tilde{φ}}_{L}^{7} \end{matrix}$ ]。由假设2可知,在伪光源位置 $\begin{matrix} {\tilde{x}}^{i} \end{matrix}$ 的值非常小,因此 $\begin{matrix} {\tilde{φ}}_{l}^{i} \end{matrix}$ =Bdiag(φ_x,_l) $\begin{matrix} {\tilde{x}}^{i} \end{matrix}$ (i=3,4,5,6,l=1,2,…,L)的每个分量也非常小。再考虑到 $\begin{matrix} {\tilde{φ}}_{l}^{7} \end{matrix}$ =Bdiag(φ_x,_l) $\begin{matrix} {\tilde{x}}^{7} \end{matrix}$ =Bdiag(φ_x,_l)0=0,因此背景荧光引起的噪声矩阵E对于Ψ的秩不会产生很大影响。最后将(11)式代公(12)式,可得:

\begin{matrix} \begin{matrix} Ψ = [p_{1}^{1} B {\tilde{x}}^{1} \dots p_{l}^{1} B {\tilde{x}}^{1} \dots p_{L}^{1} B {\tilde{x}}^{1}] + \\ [p_{1}^{2} B {\tilde{x}}^{2} \dots p_{l}^{2} B {\tilde{x}}^{2} \dots p_{L}^{2} B {\tilde{x}}^{2}] + E = \\ [B {\tilde{x}}^{1} B {\tilde{x}}^{2}] [\begin{matrix} p_{1}^{1} \dots p_{l}^{1} \dots p_{L}^{1} \\ p_{1}^{2} \dots p_{l}^{2} \dots p_{L}^{2} \end{matrix}] + E 。 (13) \end{matrix} \end{matrix}

从以上推导可以看出,对于图1所示的仿体,如果忽略噪声矩阵E对外表面观测矩阵Ψ的影响,Ψ的秩为2。该证明方法很容易被推广到荧光团具有较大体积的情况。如果一个荧光团的体积较大而不能等效成点光源,则可以将其拆分,获得可以等效成点光源的几个小荧光团,这样利用本节类似的证明方法,可以证明Ψ仍然是低秩的。这里需要说明的一点是,由于FMT的一个很重要的应用是疾病(如肿瘤)的早期检测,因此患病部位的荧光团通常应该是很小的。

2.3 提出方法

由于硬件的限制,外表面观测矩阵Ψ中的很多元素在医疗实践中是无法获得的。图2为一个数字鼠的截面图,该截面与xy平面平行。在图2中,黑色点是等效激发节点,红色的圆圈(O点)是该截面数字鼠中心,蓝色曲线为数字鼠的外表面。如果当前激发节点是绿色圆圈中的黑色点,与之反方向一定角度(对于图2是120°)内的外表面的观测是可靠的,如阴影部分区域对应的外表面,这部分区域被称为可视域(FOV),其他部分的外表面观测被认为是未知的。

图 2. 可视域示意图

Fig. 2. Sketch of FOV

下载图片查看所有图片

上述分析说明,Ψ的每一个列向量 $\begin{matrix} {\bar{φ}}_{m, l} \end{matrix}$ (l=1,2,…,L)都有很多元素缺失,通常应对这一状况的方法是,将缺失元素对应的方程去掉。原始方程组由(7)式表示,删掉缺失观测对应的方程组后为

\begin{matrix} {\hat{φ}}_{m, l} = A_{l} x, l = 1,2, \dots, L, (14) \end{matrix}

式中: $\begin{matrix} {\hat{φ}}_{m, l} \end{matrix}$ = $\begin{matrix} ({\bar{φ}}_{m, l})_{I_{l}} \end{matrix}$ ;A_l= $\begin{matrix} [Bdiag (φ_{x, l} {)]}_{I_{l} :} \end{matrix}$ 。I_l是第l个激发节点的可视域中的外表面节点对应的指标集, $\begin{matrix} ({\bar{φ}}_{m, l})_{I_{l}} \end{matrix}$ 是指标集I_l对应的 $\begin{matrix} {\bar{φ}}_{m, l} \end{matrix}$ 的分量形成的子向量, $\begin{matrix} [Bdiag (φ_{x, l} {)]}_{I_{l} :} \end{matrix}$ 是指标集I_l对应的矩阵Bdiag(φ_x,_l)的行组成的子矩阵。由于(14)式中的每一个方程组都有相同的解x,这L个方程组可以合并成一个方程组:

\begin{matrix} φ = Ax, (15) \end{matrix}

式中:φ= $\begin{matrix} [({\hat{φ}}_{m, 1})^{T} ({\hat{φ}}_{m, 2})^{T} \dots ({\hat{φ}}_{m, L})^{T}]^{T} \end{matrix}$ ;A= $\begin{matrix} [(A_{1})^{T} (A_{2})^{T} \dots (A_{L})^{T}]^{T} \end{matrix}$ 。FMT逆问题就是利用方程组(15)式求x。

从2.2节的讨论中可以看到,背景荧光将会导致噪声矩阵E。尽管E中元素的值可能很接近于0,但是由于FMT逆问题的病态性,它可能导致比较大的重建误差,可以考虑采用矩阵低秩分解的办法抑制或去除(13)式中的噪声矩阵E。

在人工智能和图像处理等很多领域,矩阵的低秩填充技术一直是研究的热点,如缺损图像的修复^[16]、电影用户的偏好分析与电影推荐^[17]等。一个相对简单的矩阵低秩填充算法OR1MP^[18]常被用于恢复外表面观测矩阵Ψ中缺失的数据,并抑制噪声矩阵E。由于本课题组使用的符号与Wang等^[18]的论文表述有所差异,为了方便阅读,本课题组给出基于新符号的OR1MP算法流程。

为了叙述简洁,首先给出一些记号。Ψ_NO为含有背景荧光噪声的外表面观测矩阵。Ω为掩模矩阵,Ω与Ψ_NO是同型矩阵。如果Ψ_NO中i行j列的元素缺失,则Ω的i行j列元素为0;如果Ψ_NO中i行j列的元素没有缺失,则Ω的i行j列元素为1。Ψ_NO中缺失元素用0替代。vec(·)为一个函数,它将一个矩阵按照一定顺序拉成一个列向量。如:Ψ_NO=[ $\begin{matrix} {\bar{φ}}_{m, 1} \end{matrix}$ … $\begin{matrix} {\bar{φ}}_{m, l} \end{matrix}$ … $\begin{matrix} {\bar{φ}}_{m, L} \end{matrix}$ ],则vec(Ψ_NO)= $\begin{matrix} [({\bar{φ}}_{m, 1})^{T} \dots ({\bar{φ}}_{m, l})^{T} \dots ({\bar{φ}}_{m, L})^{T}]^{T} \end{matrix}$ 。I_Ω是一个指标集,它记录vec(Ω)非0元素的位置。令a是一个列向量, $\begin{matrix} {(a)}_{I_{Ω}} \end{matrix}$ 表示由a对应指标集I_Ω的元素组成的子向量。令Ψ_LR是利用矩阵低秩填充方法去噪后获得的外表面观测矩阵,它与Ψ_NO是同型矩阵。

OR1MP算法^[18]的步骤如下:

步骤1:输入Ψ_NO、Ω及最大迭代次数K。令 $\begin{matrix} Ψ_{LR}^{0} \end{matrix}$ =0是与Ψ_NO同型的零矩阵。置k=0。

步骤2:对矩阵R^k=Ψ_NO- $\begin{matrix} {Ψ^{k}}_{LR} \end{matrix}$ 进行奇异值分解,并令最大奇异值对应的左右奇异向量分别为u^k和v^k,M^k=u^k $\begin{matrix} (v^{k})^{T} \end{matrix}$ 。

步骤3: 计算 $\begin{matrix} {\bar{M}}^{k} \end{matrix}$ =[ $\begin{matrix} (vec (M^{1} {))}_{I_{Ω}} \end{matrix}$ $\begin{matrix} (vec (M^{1} {))}_{I_{Ω}} \end{matrix}$ … $\begin{matrix} (vec (M^{k} {))}_{I_{Ω}} \end{matrix}$ ]。

步骤4:计算权向量θ^k= $\begin{matrix} [({\bar{M}}^{k})^{T} {\bar{M}}^{k}]^{- 1} \end{matrix} \begin{matrix} ({\bar{M}}^{k})^{T} \end{matrix}$ × $\begin{matrix} [vec (Ψ_{NO} {)]}_{I_{Ω}} \end{matrix}$ 。

步骤5:计算 $\begin{matrix} Ψ_{LR}^{k + 1} \end{matrix}$ = $\begin{matrix} \overset{k}{\sum_{i = 1}} \end{matrix} \begin{matrix} θ_{i}^{k} \end{matrix}$ Mⁱ,其中 $\begin{matrix} θ_{i}^{k} \end{matrix}$ 是θ^k的第i个分量。

步骤6:如果迭代次数未达到最大迭代次数K,则k=k+1,转步骤2;否则,算法停止,输出去噪填充后的矩阵Ψ_LR= $\begin{matrix} {Ψ^{k + 1}}_{LR} \end{matrix}$ 。

为了便于理解OR1MP算法,图3给出了该算法的流程图。

图 3. 算法1的流程图

Fig. 3. Flow chart of algorithm 1

下载图片查看所有图片

对于含背景荧光噪声且有元素缺失的外表面观测矩阵Ψ_NO,可以采用(15)式构建FMT逆问题。类似于(15)式,此时的逆问题记为

\begin{matrix} φ_{NO} = Ax 。 (16) \end{matrix}

对于去噪后的外表面观测矩阵Ψ_LR,也可以使用相同的方法,构建的FMT逆问题记为

\begin{matrix} φ_{LR} = Ax 。 (17) \end{matrix}

所提算法并没用采用填充后的完整外表面观测矩阵构建FMT逆问题,原因有两个。第一,由完整外表面观测矩阵构建的FMT逆问题规模比逆问题[(17)式]大很多,求解速度会很慢。第二,在仿真实验中,真实的外表面观测是可以获得的,不妨设可视域内的真实观测为φ_TR,利用(17)式作为逆问题,可以对比去噪前后的两个相对误差,即:

\begin{matrix} \begin{matrix} e_{NO} = {(φ_{NO} - φ_{TR})}_{2} / {(φ_{TR})}_{2}, (18) \\ e_{LR} = {(φ_{LR} - φ_{TR})}_{2} / {(φ_{TR})}_{2} 。 (19) \end{matrix} \end{matrix}

如果采用填充后的完整外表面观测矩阵构建FMT逆问题,则将没有办法获得相对应的含噪外表面的相对误差,因为可视域外的含噪观测是缺失的,是不能参与FMT光源重建的。

3 仿真实验

3.1 实验设计

采用两组非匀质数字鼠仿真实验验证所提方法的有效性。一组实验对单目标光源进行重建,另一组对双目标光源进行重建。实验采用的数字鼠模型由南加利福尼亚大学生物医学影像组提供^[19]。为了计算方便,实验中依据惯例^[20]对数字鼠模型进行了简化,沿z轴截取数字鼠躯干中部约34 mm作为实验对象,仅保留肌肉、心脏、胃、肝、肾、肺等器官,其余器官全部合并成肌肉。各器官的光学参数见表1,其中激发光的波长为670 nm,发射光波长为710 nm,μ'_sx=(1-g)μ_sx和μ'_sm=(1-g)μ_sm分别为激发阶段和发射阶段数字鼠器官的约化散射系数。

表 1. 数字鼠器官光学参数

Table 1. Optical parameters of the organs of digital mouse

Organ	$\begin{matrix} {μ_{ax}}_{} \end{matrix}$ /mm^-¹	μ' $\begin{matrix} {_{sx}}_{} \end{matrix}$ /mm^-¹	$\begin{matrix} {μ_{am}}_{} \end{matrix}$ /mm^-¹	μ' $\begin{matrix} {_{sm}}_{} \end{matrix}$ /mm^-¹	g
Muscle	0.075	0.412	0.043	0.350	0.90
Heart	0.051	0.944	0.030	0.870	0.85
Stomach	0.010	1.417	0.007	1.340	0.92
Liver	0.304	0.668	0.176	0.629	0.90
Kidneys	0.058	2.204	0.034	2.021	0.86
Lungs	0.170	2.157	0.097	2.093	0.90

查看所有表

实验采用的激发平面对应的z轴坐标为z=17 mm,位于截取躯干部位的正中间。18个等效激发节点等角度地散布在激发平面上,等效激发节点以及反方向的120°可视域如图2所示。逆问题采用2.3节所述方法构建。所提方法对应的逆问题(即利用OR1MP算法获得去背景荧光噪声后的可视域外表面观测所形成的逆问题)如(17)式所示。对比方法的逆问题则由含背景荧光噪声的可视域外表面观测形成,即(16)式。直接求解逆问题方程组(16)、(17)式的重建效果是不理想的,因为分子影像的逆问题都有非常严重的病态性。通常采用正则化方法对逆问题进行求解,求解的正则化数学模型有很多,如基于L₂-范数^[21]、L₁-范数^[22]、L₀-范数^[23]、L_p-范数(0<p<1)^[24]、TV-范数^[25]、(L₁-L₂)-范数^[26]等的正则化模型。目前,公认的在效率、精度和稳定性上都非常出色的优化模型是基于L₁-范数的数学模型,因此本节采用该模型对逆问题(16)、(17)式进行求解。对应于逆问题(16)、(17)式的L₁-范数优化模型分别为

\begin{matrix} \begin{matrix} \min_{x} \frac{1}{2} {(Ax - φ_{NO})}_{2}^{2} + τ {(x)}_{1}, (20) \\ \min_{x} \frac{1}{2} {(Ax - φ_{LR})}_{2}^{2} + τ {(x)}_{1} 。 (21) \end{matrix} \end{matrix}

式中:τ为正则参数。实验采用IVTCG(也称ITCG)算法^[27-28]对优化问题(20)、(21)式进行求解,这是一个在分子影像领域被广泛使用^[29-31]且程序公开的算法。为了叙述方便,将用于求解优化问题(20)式的IVTCG方法(即对比方法)记为IVTCG-NO,将用于求解优化问题(21)式的IVTCG方法(即所提方法)记为IVTCG-LR。

使用两个客观指标来评价上述两种方法的优劣。第一个指标是定位误差(LE),即重建目标的中心与真实目标中心的距离(单位为mm)。第二个指标是对比噪声比(CNR) $\begin{matrix} {R_{CN}}^{[32]} \end{matrix}$ ,其计算公式为

\begin{matrix} R_{CN} = \frac{μ_{noi} - μ_{nob}}{(w_{n o i} σ_{noi}^{2} + w_{nob} σ_{nob}^{2})^{\frac{1}{2}}}, (22) \end{matrix}

式中:μ_noi和μ_nob分别为感兴趣区域和背景区域节点对应能量值的均值; $\begin{matrix} σ_{noi}^{2} \end{matrix}$ 和 $\begin{matrix} σ_{nob}^{2} \end{matrix}$ 分别为感兴趣区域和背景区域节点对应能量值的方差;w_noi和w_nob分别为感兴趣区域和背景区域中节点数占总结点数的比例。计算CNR需要将截取的数字鼠躯干部分划分成两部分,第一部分为真实光源所在区域,称为感兴趣区域,第二部分为剩余的区域,称为背景区域。光源重建后,每一个有限元节点都会对应一个能量值。如果重建结果比较理想,即大部分能量大于0的节点位于感兴趣区域,且能量值分布比较均匀,只有少部分能量大于0的节点分布在背景区域,此时μ_noi很大,μ_nob较小,而由于感兴趣区域的节点能量值分布均匀,背景区域节点能量值大部分为0,故方差 $\begin{matrix} σ_{noi}^{2} \end{matrix}$ 和 $\begin{matrix} σ_{nob}^{2} \end{matrix}$ 都会比较小,因此获得的CNR会很大。反之,如果重建效果不理想,大部分能量大于0的节点位于背景区域,此时μ_noi很小,μ_nob较大,相应的方差 $\begin{matrix} σ_{nob}^{2} \end{matrix}$ 比较大,所以获得的CNR就会很小,甚至为负值。从上面的分析可以看到,指标CNR越大越好,LE越小越好。

IVTCG算法有两个参数,一个是正则参数τ,另一个是终止阈值ε,它们的取值范围分别为τ∈{10^-¹⁰,10^-¹¹,…,10^-¹⁷},ε∈{10^-⁹,10^-¹⁰,…,10^-¹⁵},参数组(τ,ε)共有56种可能的取值组合。在光源重建中,IVTCG-NO和IVTCG-LR分别在56个参数组合条件下运行,最优的LE对应的实验结果将采用表格的方式加以展示。对于矩阵低秩填充算法OR1MP来说,它也有一个参数,就是最大迭代次数K,实验均选择对应最优的相对误差值,见(19)式,对应的K为OR1MP算法的参数。

在实验中,真实目标都是圆柱体,其底面半径为1 mm,高为2 mm。真实目标位置的荧光产额设置为0.06。实验假设荧光探针在肝和肾部位富集,且分布均匀,由它们产生背景噪声,这两个器官处的背景荧光产额被设置为0.06/50、0.06/60、0.06/70、0.06/80、0.06/90、0.06/100共6个级别的不同噪声。不同噪声级别、不同光源个数下的重建结果将在后文分别说明。

这一研究的所有实验都在惠普Z420工作站上进行,该工作站的中央处理器(CPU)为至强E5-1680v2,其主频为3.0 GHz,内存为32 G。使用的程序在Windows 7环境下由MATLAB 2013b软件编写并运行。

3.2 单光源实验结果

下面通过重建单荧光目标检验所提方法的效果。在仿真实验中,需要使用两个有限元网格离散截取的数字鼠躯干部分。一个网格较细,用于模拟激发光在数字鼠体内的传播、激发光在真实光源和背景光源处转化形成发射光以及发射光在数字鼠体内传播的过程,目的是获得外表面观测数据。这个网格被称为前向网格。另一个网格相对较粗,被称为逆向网格,利用这个网格和前面获得的外表面观测数据,可以通过FMT逆问题重建光源。本实验将单光源设置在肝脏处,真实中心坐标为(17,14.3,15),前向网格含有25271个节点和141126个四面体,逆向网格含有2981个节点和15553个四面体。重建结果如表2所示,其中第1列是使用方法的名称,第2列是噪声级别(NL),第3列是两种方法各自最优的LE对应的参数组合,第4、第5和第6列分别为最优参数组合条件下的平均LE(±标准差)、平均CNR(±标准差)和平均算法运行时间(±标准差),每一列最优的重建结果用黑体表示。正常条件下OR1MP和IVTCG都是确定算法,不具有随机性,但是由于在程序中使用了MATLAB的svds.mat函数,该函数有一定的随机性,即对于相同的矩阵,其奇异值分解(SVD)结果会不同,这会导致同一噪声条件下多次运行IVTCG-LR得到的重建结果有所差异,因此表2中的数据都是5次随机实验平均后的结果,用于测试IVTCG-LR的稳健性。从表2中的数据可以看到,IVTCG-LR在LE和CNR两个指标上都有非常大的优势,两种方法的运行时间相差不多。值得注意的是,在两个最小的背景荧光噪声级别0.006/90和0.006/100下,IVTCG-LR的CNR突然变小了。考虑到FMT逆问题的病态性,这种数据上较大幅度的退化是正常的,即便退化比较严重,与比较方法相比,IVTCG-LR仍然是具有优势的。从稳健性测试结果也可以看出,由svds.mat函数导致的随机误差对IVTCG-LR基本没有影响。

表 2. 不同噪声级别下两种方法的重建结果

Table 2. Reconstruction using two methods at different noise levels

Method	NL	(τ,ε)	LE /mm	CNR	Time /s
IVTCG-LR	0.06/50	(10^-12,10^-12)	0.4679±9.2×10^-3	20.5455±1.8	71.0230±4.4
IVTCG-NO	0.06/50	(10^-14,10^-12)	0.7528±0	13.7432±0	74.3608±3.2
IVTCG-LR	0.06/60	(10^-12,10^-12)	0.4409±8.7×10^-3	25.2068±2.0	62.6312±4.3
IVTCG-NO	0.06/60	(10^-14,10^-12)	0.7420±0	14.2688±0	69.7136±1.3
IVTCG-LR	0.06/70	(10^-12,10^-12)	0.4154±1.1×10^-2	27.9832±3.3×10^-1	61.4392±7.0
IVTCG-NO	0.06/70	(10^-14,10^-12)	0.7313±0	14.8315±0	68.2468±0.6
IVTCG-LR	0.06/80	(10^-12,10^-12)	0.3928±3.9×10^-3	28.8322±4.0×10^-1	52.3152±3.5
IVTCG-NO	0.06/80	(10^-14,10^-12)	0.7134±0	16.1134±0	68.2628±1.9
IVTCG-LR	0.06/90	(10^-12,10^-13)	0.3810±4.4×10^-4	17.4285±1.2×10^-2	51.3362±4.8
IVTCG-NO	0.06/90	(10^-14,10^-12)	0.7031±0	16.1723±0	67.9030±1.5
IVTCG-LR	0.06/100	(10^-12,10^-13)	0.3591±3.9×10^-4	18.3511±5.9×10^-3	45.3048±1.9
IVTCG-NO	0.06/100	(10^-14,10^-12)	0.6915±0	16.4278±0	66.2242±0.5

查看所有表

图 4. 单光源实验中,0.06/50噪声级别下两种方法的2D和3D重建结果。(a) IVTCG-NO方法的2D重建结果;(b) IVTCG-LR方法的2D重建结果; (c) IVTCG-NO方法的3D重建结果; (d) IVTCG-LR方法的3D重建结果

Fig. 4. Reconstruction using two methods at noise level of 0.06/50 in single light source experiment. (a) 2D reconstruction using IVTCG-NO method; (b) 2D reconstruction using IVTCG-LR method; (c) 3D reconstruction using IVTCG-NO method; (d) 3D reconstruction using IVTCG-LR method

下载图片查看所有图片

下面在视觉上展示两种方法的重建结果。受篇幅所限,本组实验只给出0.06/50噪声级别下,5次实验的第1次实验的对比结果。图4给出了该次实验中两种方法的2D和3D重建结果。图中的第一行是2D重建结果图,该展示平面垂直于z轴,对应的z轴坐标为z=15 mm,与真实光源中心的z轴坐标相同。图中红色的圆圈为真实光源位置,不同的颜色对应不同的荧光产额密度。图4(c)、(d)为3D重建结果图,其中灰色区域对应肌肉,蓝色区域对应心脏,青色区域对应肺,绿色区域对应胃,褐色区域对应肾,黄色区域对应肝。图中的红色圆柱体是真实光源,它位于肝脏内,粉色区域是重建的荧光目标。

从图4的展示结果可以看出,两种方法的重建光源都可以很好地定位真实目标,但IVTCG-NO方法的重建目标区域体积非常大,在背景区域上的分布非常广,因此相应的CNR必然会相对较小,这与表2中的数据指标是一致的。

最后给出去噪前和去噪后外表面观测的平均相对误差和OR1MP算法的平均运行时间,两个误差分别利用(18)式和(19)式进行计算,结果如表3所示。可以看到,在利用OR1MP算法去除背景荧光噪声后,外表面观测的相对误差e出现了大幅降低。由于需要填充的矩阵较小,涉及到的外表面观测矩阵只有数百行18列,因此OR1MP算法的运行速度非常快,运行时间可以忽略不计。

表 3. OR1MP算法的运行时间和去噪前后外表面观测的相对误差

Table 3. Running time of OR1MP algorithm and relative errors of the boundary observation vectors before and after denoising operation

NL	0.06/50	0.06/60	0.06/70	0.06/80	0.06/90	0.06/100
Time /s	0.0052	0.0052	0.0058	0.0062	0.0056	0.0054
e_LR	1.0172	0.8858	0.7965	0.7331	0.6863	0.6508
e_NO	1.4944	1.2453	1.0674	0.9340	0.8302	0.7472

查看所有表

3.3 双光源实验结果

下面通过重建双荧光目标检验IVTCG-LR的效果。该实验将两个光源均设置在肝脏处,真实光源的中心坐标分别为(19,10.5,15)和(15,13.5,15),两个圆柱体光源的边缘相距3 mm。前向网格包含26246个节点和146735个四面体,逆向网格包含2981个节点和15553个四面体。对于每个实验,两种方法均运行5次,重建结果如表4和表5所示。

表4中每一列最优的重建结果用黑体表示。从表4的数据可以看到,IVTCG-LR方法在LE和CNR两个指标上都有非常大的优势。对于不同的噪声级别,两种方法的重建结果变化都不大,且重建结果都有随着噪声减小而变差的趋势。考虑到数据的变化确实非常小,本课题组认为这是正常的数据波动造成的。在运行时间方面,IVTCG-LR方法的运行时间约为IVTCG-NO方法的3倍,但其LE只为IVTCG-NO方法的1/3,在精度上远远优于IVTCG-NO方法。本课题组认为牺牲一定的运行时间来换取3倍的精度是完全值得的。

表5列出了IVTCG-LR和IVTCG-NO方法获得对应于第1和第2个重建光源的LE(LE1, LE2),最优的重建结果用黑体表示。从表5的数据中可以看出,对于每一个噪声级别,对于两个光源中的(任意的)单独一个光源,IVTCG-LR方法的定位精度优于IVTCG-NO方法。

下面在视觉上展示两种方法的重建结果。与3.2节的做法类似,本组实验也只给出0.06/50噪声级别下,5次实验中第1次实验的对比结果。图5给出了该次实验中两种方法的2D和3D重建结果。从2D展示图可以看出,IVTCG-LR方法的定位更加准确;从3D展示图可以看到,IVTCG-NO方法重建的光源在背景区域上分布得更加广泛,因此CNR相对较低,这与表4和表5的数值结果是一致的。

表 4. 不同噪声级别条件下两种方法的重建结果

Table 4. Reconstruction using two methods at different noise levels

Method	NL	(τ,ε)	LE /mm	CNR	Time /s
IVTCG-LR	0.06/50	(10^-14,10^-13)	0.8234±1.9×10^-3	18.8590±6.1×10^-2	148.0212±4.4
IVTCG-NO	0.06/50	(10^-16,10^-12)	2.7076±0	13.5230±0	43.6608±1.4
IVTCG-LR	0.06/60	(10^-14,10^-13)	0.8443±8.6×10^-4	18.6281±4.6×10^-2	146.4514±3.2
IVTCG-NO	0.06/60	(10^-15,10^-12)	2.7227±0	13.2980±0	48.4676±2.2
IVTCG-LR	0.06/70	(10^-14,10^-13)	0.8650±2.4×10^-3	18.2852±6.9×10^-2	146.5268±2.1
IVTCG-NO	0.06/70	(10^-15,10^-12)	2.7368±0	13.2485±0	47.4714±1.1
IVTCG-LR	0.06/80	(10^-14,10^-13)	0.8837±1.8×10^-3	18.1829±9.2×10^-2	151.7274±1.3
IVTCG-NO	0.06/80	(10^-16,10^-12)	2.7259±0	13.3588±0	51.4532±1.3
IVTCG-LR	0.06/90	(10^-14,10^-13)	0.9003±3.1×10^-3	17.8455±1.9×10^-1	147.2034±3.8
IVTCG-NO	0.06/90	(10^-16,10^-12)	2.7150±0	13.5971±0	50.9520±0.3
IVTCG-LR	0.06/100	(10^-15,10^-12)	0.9110±2.8×10^-3	18.6627±1.6×10^-1	70.1896±11
IVTCG-NO	0.06/100	(10^-17,10^-12)	2.7343±0	13.1110±0	48.5222±7.2

查看所有表

表 5. 光源1和光源2的LE

Table 5. LE of light source 1 and light source 2 mm

NL	0.06/50	0.06/60	0.06/70	0.06/80	0.06/90	0.06/100
IVTCG-LR(LE1)	0.4862±1.0×10^-3	0.5057±8.1×10^-4	0.5247±1.0×10^-3	0.5423±9.4×10^-4	0.5579±2.4×10^-3	0.5787±3.1×10^-3
IVTCG-LR(LE2)	0.3372±1.5×10^-3	0.3386±8.4×10^-4	0.3403±1.5×10^-3	0.3414±9.2×10^-4	0.3424±9.9×10^-4	0.3224±2.1×10^-3
IVTCG-NO(LE1)	1.8779±0	1.8957±0	1.9109±0	1.8988±0	1.8767±0	1.9088±0
IVTCG-NO(LE2)	0.8297±0	0.8270±0	0.8259±0	0.8271±0	0.8383±0	0.8255±0

查看所有表

图 5. 双光源实验中,0.06/50噪声级别下两种方法的重建结果。(a) IVTCG-NO方法的2D重建结果; (b) IVTCG-LR方法的2D重建结果; (c) IVTCG-NO方法的2D重建结果; (d) IVTCG-LR方法的3D重建结果

Fig. 5. Reconstruction using two methods at noise level of 0.06/50 in double light sources. (a) 2D reconstruction using IVTCG-NO method; (b) 2D reconstruction using IVTCG-LR method; (c) 3D reconstruction using IVTCG-NO method; (d) 3D reconstruction using IVTCG-LR method

下载图片查看所有图片

最后给出双荧光目标条件下,去噪前和去噪后外表面观测的平均相对误差和OR1MP算法的平均运行时间,两个误差分别利用(18)、(19)式计算。表6给出了相关结果。类似于3.2节的结论,在利用OR1MP算法去除背景荧光噪声之后,外表面观测的相对误差有了大幅降低,OR1MP算法的运行时间可以忽略不计。

表 6. OR1MP算法的运行时间和去噪前后外表面观测的相对误差

Table 6. Running time of OR1MP algorithm and relative errors of the boundary observation vectors before and after denoising operation

NL	0.06/50	0.06/60	0.06/70	0.06/80	0.06/90	0.06/100
Time /s	0.0046	0.0056	0.0050	0.0046	0.0054	0.0056
e_LR	1.3504	1.1595	1.0280	0.9332	0.8625	0.8082
e_NO	2.0465	1.7054	1.4618	1.2790	1.1369	1.0232

查看所有表

4 结论

提出了一种基于低秩矩阵填充技术的背景荧光噪声去除方法,该方法在较高背景荧光噪声条件下能很好地重建FMT荧光目标。本课题组深入研究了利用有限元方法重建FMT荧光目标的基础理论,通过一系列推导,得出了外表面观测矩阵是一个具有缺失元素的低秩矩阵这一结论,这使得利用低秩矩阵填充方法抑制或去除外表面背景荧光噪声的方法有了充分的理论基础。与以往的去除背景荧光的方法相比,这是一个解决该类问题的全新视角。数值实验表明,利用低秩矩阵填充方法去噪后,外表面观测矩阵包含的背景荧光噪声有了较大幅度降低,重建的荧光目标精度有了显著提高。

本课题组使用一个通用的低秩矩阵填充方法来重建外表面观测矩阵。从表3和表6的数据可以看出,尽管去噪后背景荧光噪声有了很大降低,但噪声仍然很大。因此,结合FMT问题的特点设计新的低秩矩阵填充算法,进一步降低背景荧光噪声,提高重建精度,是本课题组未来的研究工作。

参考文献

[1] Zhang J, Shi J, Guang H, et al. Iterative correction scheme based on discrete cosine transform and L1 regularization for fluorescence molecular tomography with background fluorescence[J]. IEEE Transactions on Biomedical Engineering, 2016, 63(6): 1107-1115.

[2] Gao M, Lewis G, Turner G M, et al. Effects of background fluorescence in fluorescence molecular tomography[J]. Applied Optics, 2005, 44(26): 5468-5474.

[3] Soubret A, Ntziachristos V. Fluorescence molecular tomography in the presence of background fluorescence[J]. Physics in Medicine & Biology, 2006, 51(16): 3983-4001.

[4] Ale A, Ermolayev V, Deliolanis N, et al. Fluorescence background subtraction technique for hybrid fluorescence molecular tomography/X-ray computed tomography imaging of a mouse model of early stage lung cancer[J]. Journal of Biomedical Optics, 2013, 18(5): 056006.

[5] Qin B J, Hu C, Huang S S. Target/background classification regularized nonnegative matrix factorization for fluorescence unmixing[J]. IEEE Transactions on Instrumentation and Measurement, 2016, 65(4): 874-889.

[6] Xu H, Rice B. In-vivo fluorescence imaging with a multivariate curve resolution spectral unmixing technique[J]. Journal of Biomedical Optics, 2009, 14(6): 064011.

[7] Plaza J. Hendrix E M T, Garcia I, et al. On endmember identification in hyperspectral images without pure pixels: a comparison of algorithms[J]. Journal of Mathematical Imaging and Vision, 2012, 42(2/3): 163-175.

[8] Pu H S, Zhang G L, He W, et al. Resolving fluorophores by unmixing multispectral fluorescence tomography with independent component analysis[J]. Physics in Medicine & Biology, 2014, 59(17): 5025-5042.

[9] Han D, Yang X, Liu K, et al. Efficient reconstruction method for L1 regularization in fluorescence molecular tomography[J]. Applied Optics, 2010, 49(36): 6930-6937.

[10] Chen X L, Yang D F, Zhang Q T, et al. L1/2 regularization based numerical method for effective reconstruction of bioluminescence tomography[J]. Journal of Applied Physics, 2014, 115(18): 184702.

[11] Ntziachristos V. Fluorescence molecular imaging[J]. Annual Review of Biomedical Engineering, 2006, 8(1): 1-33.

[12] Cong A, Wang G. A finite-element-based reconstruction method for 3D fluorescence tomography[J]. Optics Express, 2005, 13(24): 9847-9857.

[13] Wang D F, Liu Y, Chen Y P. A novel finite-element-based algorithm for fluorescence molecular tomography of heterogeneous media[J]. IEEE Transactions on Information Technology in Biomedicine, 2009, 13(5): 766-773.

[14] Cong W, Wang G, Kumar D, et al. Practical reconstruction method for bioluminescence tomography[J]. Optics Express, 2005, 13(18): 6756-6771.

[15] 易黄建. 基于正则化的荧光分子断层成像重建方法研究[D]. 西安: 西安电子科技大学, 2013: 15- 31.

Yi HJ. Regularization based reconstruction algorithms for fluorescence molecular tomography[D]. Xi’an: Xidian University, 2013: 15- 31.

[16] Liu P, Lewis J, Rhee T. Low-rank matrix completion to reconstruct incomplete rendering images[J]. IEEE Transactions on Visualization and Computer Graphics, 2017, 24(8): 2353-2365.

[17] Bhaskar S. Probabilistic low-rank matrix completion from quantized measurements[J]. Journal of Machine Learning Research, 2016, 17(1): 1-34.

[18] Wang Z, Lai M J, Lu Z S, et al. Orthogonal rank-one matrix pursuit for low rank matrix completion[J]. SIAM Journal on Scientific Computing, 2015, 37(1): A488-A514.

[19] Dogdas B, Stout D, Chatziioannou A, et al. Digimouse: a 3D whole body mouse atlas from CT and cryosection data[J]. Physics in Medicine & Biology, 2007, 52(3): 577-587.

[20] He X L, Wang X D, Yi H J, et al. Laplacian manifold regularization method for fluorescence molecular tomography[J]. Journal of Biomedical Optics, 2017, 22(4): 045009.

[21] Zhang X F, Badea C T, Johnson G A. Three-dimensional reconstruction in free-space whole-body fluorescence tomography of mice using optically reconstructed surface and atlas anatomy[J]. Journal of Biomedical Optics, 2009, 14(6): 064010.

[22] Shi J W, Liu F, Zhang G L, et al. Enhanced spatial resolution in fluorescence molecular tomography using restarted L1-regularized nonlinear conjugate gradient algorithm[J]. Journal of Biomedical Optics, 2014, 19(4): 046018.

[23] Ye J Z, Chi C W, Xue Z W, et al. Fast and robust reconstruction for fluorescence molecular tomography via a sparsity adaptive subspace pursuit method[J]. Biomedical Optics Express, 2014, 5(2): 387-406.

[24] Guo H B, Yu J J, He X W, et al. Improved sparse reconstruction for fluorescence molecular tomography with L1/2 regularization[J]. Biomedical Optics Express, 2015, 6(5): 1648-1664.

[25] Feng J C, Qin C H, Jia K B, et al. Total variation regularization for bioluminescence tomography with the split Bregman method[J]. Applied Optics, 2012, 51(19): 4501-4512.

[26] 张海波, 耿国华, 赵映程, 等. 基于非凸L1-2正则子的锥束X射线发光断层成像[J]. 光学学报, 2017, 37(6): 0617001.

Zhang H B, Geng G H, Zhao Y C, et al. Nonconvex L1-2 regularization for fast cone-beam X-ray luminescence computed tomography[J]. Acta Optica Sinica, 2017, 37(6): 0617001.

[27] He X W, Liang J M, Wang X R, et al. Sparse reconstruction for quantitative bioluminescence tomography based on the incomplete variables truncated conjugate gradient method[J]. Optics Express, 2010, 18(24): 24825-24841.

[28] Wang X D, Liu F, Jiao L C, et al. Incomplete variables truncated conjugate gradient method for signal reconstruction in compressed sensing[J]. Information Sciences, 2014, 288: 387-411.

[29] Zhang H, Geng G, Wang X, et al. Fast and robust reconstruction for fluorescence molecular tomography via L1-2 regularization[J]. BioMed Research International, 2016, 2016: 5065217.

[30] Liu M, Guo H, Liu H, et al. In vivo pentamodal tomographic imaging for small animals[J]. Biomedical Optics Express, 2017, 8(3): 1356-1371.

[31] 张旭, 易黄建, 侯榆青, 等. 基于局部保留投影的荧光分子断层成像快速重建[J]. 光学学报, 2016, 36(7): 0717001.

Zhang X, Yi H J, Hou Y Q, et al. Fast reconstruction in fluorescence molecular tomography based on locality preserving projections[J]. Acta Optica Sinica, 2016, 36(7): 0717001.

[32] Song X M, Pogue B W, Jiang S D, et al. Automated region detection based on the contrast-to-noise ratio in near-infrared tomography[J]. Applied Optics, 2004, 43(5): 1053-1062.

王晓东, 耿国华, 易黄建, 何雪磊, 贺小伟. 基于低秩矩阵填充的背景荧光噪声抑制方法[J]. 光学学报, 2018, 38(10): 1017003. Xiaodong Wang, Guohua Geng, Huangjian Yi, Xuelei He, Xiaowei He. Low-Rank-Matrix-Completion-Based Method for Suppressing Background Fluorescence[J]. Acta Optica Sinica, 2018, 38(10): 1017003.

基于低秩矩阵填充的背景荧光噪声抑制方法下载： 967次

1 引言

2 原理与方法

2.1 光辐射传输模型及其简化方程

2.2 外表面观测矩阵的低秩属性

图 1. 2D仿真模型

Fig. 1. 2D simulation model

2.3 提出方法

图 2. 可视域示意图

Fig. 2. Sketch of FOV

图 3. 算法1的流程图

Fig. 3. Flow chart of algorithm 1

3 仿真实验

3.1 实验设计

表 1. 数字鼠器官光学参数

Table 1. Optical parameters of the organs of digital mouse

3.2 单光源实验结果

表 2. 不同噪声级别下两种方法的重建结果

Table 2. Reconstruction using two methods at different noise levels

图 4. 单光源实验中,0.06/50噪声级别下两种方法的2D和3D重建结果。(a) IVTCG-NO方法的2D重建结果;(b) IVTCG-LR方法的2D重建结果; (c) IVTCG-NO方法的3D重建结果; (d) IVTCG-LR方法的3D重建结果

Fig. 4. Reconstruction using two methods at noise level of 0.06/50 in single light source experiment. (a) 2D reconstruction using IVTCG-NO method; (b) 2D reconstruction using IVTCG-LR method; (c) 3D reconstruction using IVTCG-NO method; (d) 3D reconstruction using IVTCG-LR method

表 3. OR1MP算法的运行时间和去噪前后外表面观测的相对误差

Table 3. Running time of OR1MP algorithm and relative errors of the boundary observation vectors before and after denoising operation

3.3 双光源实验结果

表 4. 不同噪声级别条件下两种方法的重建结果

Table 4. Reconstruction using two methods at different noise levels

表 5. 光源1和光源2的LE

Table 5. LE of light source 1 and light source 2 mm

图 5. 双光源实验中,0.06/50噪声级别下两种方法的重建结果。(a) IVTCG-NO方法的2D重建结果; (b) IVTCG-LR方法的2D重建结果; (c) IVTCG-NO方法的2D重建结果; (d) IVTCG-LR方法的3D重建结果

Fig. 5. Reconstruction using two methods at noise level of 0.06/50 in double light sources. (a) 2D reconstruction using IVTCG-NO method; (b) 2D reconstruction using IVTCG-LR method; (c) 3D reconstruction using IVTCG-NO method; (d) 3D reconstruction using IVTCG-LR method

表 6. OR1MP算法的运行时间和去噪前后外表面观测的相对误差

Table 6. Running time of OR1MP algorithm and relative errors of the boundary observation vectors before and after denoising operation

4 结论

Article Outline

关于本站 Cookie 的使用提示

全站搜索

基于低秩矩阵填充的背景荧光噪声抑制方法 下载： 967次

1 引言

2 原理与方法

2.1 光辐射传输模型及其简化方程

2.2 外表面观测矩阵的低秩属性

图 1. 2D仿真模型

Fig. 1. 2D simulation model

2.3 提出方法

图 2. 可视域示意图

Fig. 2. Sketch of FOV

图 3. 算法1的流程图

Fig. 3. Flow chart of algorithm 1

3 仿真实验

3.1 实验设计

表 1. 数字鼠器官光学参数

Table 1. Optical parameters of the organs of digital mouse

3.2 单光源实验结果

表 2. 不同噪声级别下两种方法的重建结果

Table 2. Reconstruction using two methods at different noise levels

图 4. 单光源实验中,0.06/50噪声级别下两种方法的2D和3D重建结果。(a) IVTCG-NO方法的2D重建结果;(b) IVTCG-LR方法的2D重建结果; (c) IVTCG-NO方法的3D重建结果; (d) IVTCG-LR方法的3D重建结果

Fig. 4. Reconstruction using two methods at noise level of 0.06/50 in single light source experiment. (a) 2D reconstruction using IVTCG-NO method; (b) 2D reconstruction using IVTCG-LR method; (c) 3D reconstruction using IVTCG-NO method; (d) 3D reconstruction using IVTCG-LR method

表 3. OR1MP算法的运行时间和去噪前后外表面观测的相对误差

Table 3. Running time of OR1MP algorithm and relative errors of the boundary observation vectors before and after denoising operation

3.3 双光源实验结果

表 4. 不同噪声级别条件下两种方法的重建结果

Table 4. Reconstruction using two methods at different noise levels

表 5. 光源1和光源2的LE

Table 5. LE of light source 1 and light source 2 mm

图 5. 双光源实验中,0.06/50噪声级别下两种方法的重建结果。(a) IVTCG-NO方法的2D重建结果; (b) IVTCG-LR方法的2D重建结果; (c) IVTCG-NO方法的2D重建结果; (d) IVTCG-LR方法的3D重建结果

Fig. 5. Reconstruction using two methods at noise level of 0.06/50 in double light sources. (a) 2D reconstruction using IVTCG-NO method; (b) 2D reconstruction using IVTCG-LR method; (c) 3D reconstruction using IVTCG-NO method; (d) 3D reconstruction using IVTCG-LR method

表 6. OR1MP算法的运行时间和去噪前后外表面观测的相对误差

Table 6. Running time of OR1MP algorithm and relative errors of the boundary observation vectors before and after denoising operation

4 结论

Article Outline

相关论文

相关资讯

关于本站 Cookie 的使用提示

全站搜索

基于低秩矩阵填充的背景荧光噪声抑制方法下载： 967次