CN101342075B - 基于单视图的多光谱自发荧光断层成像重建方法 - Google Patents

基于单视图的多光谱自发荧光断层成像重建方法 Download PDF

Info

Publication number
CN101342075B
CN101342075B CN2008101168184A CN200810116818A CN101342075B CN 101342075 B CN101342075 B CN 101342075B CN 2008101168184 A CN2008101168184 A CN 2008101168184A CN 200810116818 A CN200810116818 A CN 200810116818A CN 101342075 B CN101342075 B CN 101342075B
Authority
CN
China
Prior art keywords
tau
light source
phi
meas
mesh generation
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.)
Expired - Fee Related
Application number
CN2008101168184A
Other languages
English (en)
Other versions
CN101342075A (zh
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.)
Beijing University of Technology
Original Assignee
Beijing University of Technology
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 Beijing University of Technology filed Critical Beijing University of Technology
Priority to CN2008101168184A priority Critical patent/CN101342075B/zh
Publication of CN101342075A publication Critical patent/CN101342075A/zh
Application granted granted Critical
Publication of CN101342075B publication Critical patent/CN101342075B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Investigating Or Analysing Materials By Optical Means (AREA)

Abstract

基于单视图的多光谱自发荧光断层成像重建方法属于光学分子影像成像领域。该方法基于扩散方程模型,并考虑了小动物体的非均匀特性,同时也考虑自发荧光光源的光谱及实际应用的特点。为达到此目的,本发明提出的基于单视图的多光谱自发荧光断层成像重建方法的实现步骤主要包括:1)数据获取;2)有限元离散化处理;3)优化可行光源区域选择;4)光源的重建。本发明克服了自发荧光断层成像重建方法重建光源不准确,重建速度比较慢,不利于实时处理以及在多光谱数据采集时可能存在误差的缺陷。

Description

基于单视图的多光谱自发荧光断层成像重建方法
技术领域
本发明涉及到光学分子影像成像模态-自发荧光断层成像(BLT)技术,尤其涉及到一种基于单视图的多光谱自发荧光断层成像重建方法。
背景技术
在过去的几年里,分子影像由于能在体揭示分子和细胞的信息而受到了越来越多的关注。因此它已成为一种疾病诊断、药物疗效评价等的重要工具。作为一种小动物分子影像成像模态,光学成像技术特别是自发荧光断层成像因为其高性能、低价格和无创伤等特性已经得到了迅速发展和广泛的应用。自发荧光断层成像是最近刚发展起来的,是一种用来在活体小动物体内重建自发荧光光源分布的光学成像技术。当荧光素与荧光素酶在有氧气和三磷酸腺苷(ATP)的条件下,就会发出荧光。而因为荧光素酶含有萤火虫素酶(firefly),磕头虫素酶(click beetle),所以发出的荧光有不同的光谱,波长一般在400nm-750nm。发出的荧光穿透生物体而到达生物体表面,然后用高灵敏度的液氮制冷电荷耦合器件(CCD)获得。CCD获得的动物表面的数据形成了BLT重建的基础。为采集整个生物体的表面数据,通常将生物体每隔90度旋转一次,用CCD探测器来获取生物体的前后左右四个视图。当光源在生物体的位置比较深时,采集一个视图的数据需要的时间大约是5-10分钟。而当采集时间超过一个小时之后荧光强度就会衰减,所以在多谱的情况下,如果在每个单谱段都采集四个视图的数据,那么采集时间远远超过1个小时,可能就会使得采集的数据不够准确。另外,在实际情况下存在类似像X射线(x-ray),计算机断层成像(CT)一样的物理限制,采集数据时就受到角度限制,因此采集到的数据仅是动物体表面的一部分。多谱采集也会使得系统矩阵的维数增加,导致计算量过大。如何降低测量时间和减小系统矩阵的维数是目前的一个难点。另一方面,在进行光源重建时,因为BLT自身具有的不适定特点导致重建结果不准确。
BLT的目标是能在体实时连续的观测。因此,发展快速的重建方法就成为亟待解决的问题。目前的重建方法常用的是有限元方法。理论上有限元网格越细,得到的结果越好,但重建的时间就越长。此外,BLT是病态问题,如何降低病态性仍需要进一步的探索。
发明内容
本发明的目的在于克服了自发荧光断层成像重建方法重建光源不准确,重建速度比较慢,不利于实时处理以及在多光谱数据采集时可能存在误差的缺陷,提出了一种基于单视图的多光谱自发荧光断层成像重建方法,该方法基于扩散方程模型,并考虑了小动物体的非均匀特性,同时也考虑自发荧光光源的光谱及实际应用的特点。
为达到此目的,本发明提出的基于单视图的多光谱自发荧光断层成像重建方法的实现步骤主要包括:1)数据获取;2)有限元离散化处理;3)优化可行光源区域选择;4)光源的重建。
1)数据获取是指利用CCD探测器在被测物体表面对逃出的光子进行探测,获得边界处的光流密度;
有限元离散化处理,是指利用有限元方法将扩散方程离散为矩阵方程,并把重建区域网格剖分为{T1,...,Tk,...}的网格序列;
优化可行光源区域的选择,主要是在对扩散方程离散化处理的基础上,利用迭代方法快速选择光源大致存在的区域,即优化可行光源区域;
光源的重建,主要是在优化可行光源区域选择的基础上,确立未知光源密度变量与被测物体表面测量值之间的线性关系,并利用Tikhonov正则化方法建立目标函数,并最终求解光源密度分布。
本发明具体包括以下步骤:
(1)数据获取:
在多光谱下,用滤波片将荧光波长λ离散为m个波段τ1,...τm,其中τl=[λl-1,λl],l=1,2,...m-1,τm=[λm-1,λm],在多光谱情况下m=5即可提供足够的先验信息;利用CCD探测器在被测物体的表面对逃出的光子进行探测,获得每个波段τl的光流密度Q(x,τl),其中CCD探测到的表面为γ,x表示在被测物体中的位置;根据下面公式计算被测物体表面的光子流密度Φmeas(x,τl):
Q ( x , τ l ) = - D ( x , τ l ) ( v ( x ) · ▿ Φ meas ( x , τ l ) ) = Φ meas ( x , τ l ) 2 A ( x ; n , n ′ ) ( x ∈ γ )
D(x,τl)=1/(3(μa(x,τl)+(1-g)μs(x,τl)))是生物组织的扩散系数,其中μa(x,τl)是生物组织的吸收系数,μs(x,τl)是生物组织散射系数,g是生物组织各向异性参数;v(x)是被测物体的表面
Figure G2008101168184D00031
的单位法向量;A(x;n,n′)≈(1+R(x))/(1-R(x)),n′为外界媒质的折射系数;当外界媒质是空气时,R(x)可以近似为R(x)≈-1.4399n-2+0.7099n-1+0.6681+0.0636n,其中n为空气的折射率;
(2)有限元离散化处理:
光子在生物组织中传输的数学模型用下面的扩散方程表示:
- ▿ · ( D ( x , τ l ) ▿ Φ ( x , τ l ) ) + μ a ( x , τ l ) Φ ( x , τ l ) = S ( x , τ l ) ( x ∈ Ω )
Φ ( x , τ l ) + 2 A ( x ; n , n ′ ) D ( x , τ l ) ( v ( x ) · ▿ Φ ( x , τ l ) ) = 0 ( x ∈ ∂ Ω )
其中Ω是被测物体;
Figure G2008101168184D00034
是被测物体的表面;Φ(x,τl)是光子流密度;S(x,τl)是光源密度;根据有限元理论,得到下面的弱解形式:
∫ Ω ( D ( x , τ l ) ( ▿ Φ ( x , τ l ) ) · ( ▿ Ψ ( x , τ l ) ) + μ a ( x , τ l ) Φ ( x , τ l ) Ψ ( x , τ l ) ) dx
+ ∫ ∂ Ω 1 2 A n ( x ) Φ ( x , τ l ) Ψ ( x , τ l ) dx = ∫ Ω S ( x , τ l ) Ψ ( x , τ l ) dx ( ∀ Ψ ( x , τ l ) ∈ H 1 ( Ω ) )
H1(Ω)是Sobelev空间,Ψ(x,τl)是任意给定的试探函数。在自适应有限元分析的框架下,基于自适应的网格细分,把重建区域网格剖分表示为{T1,...,Tk,...}的网格序列,其中网格剖分Tk包括
Figure G2008101168184D00037
个离散单元。利用有限元方法,对弱解形式进行离散,得到下面矩阵方程:
(Kkl)+Ckl)+Bkl))Φkl)=Fkl)Skl)
矩阵元素的具体形式为:
k ij k ( τ l ) = ∫ Ω D ( x , τ l ) ( ▿ Ψ i k ( x , τ l ) ) · ( ▿ Ψ j k ( x , τ l ) ) dx c ij k ( τ l ) = ∫ Ω μ a ( x , τ l ) Ψ i k ( x , τ l ) Ψ j k ( x , τ l ) dx b ij k ( τ l ) = ∫ ∂ Ω Ψ i k ( x , τ l ) Ψ j k ( x , τ l ) / ( 2 A ( x ; n , n ′ ) ) dx f ij k ( τ l ) = ∫ Ω Ψ i k ( x , τ l ) Ψ j k ( x , τ l ) dx
令Kkl)+Ckl)+Bkl)=Mkl),Mkl)是稀疏正定矩阵,所以得到:
Φ k ( τ l ) = M k - 1 ( τ l ) F k ( τ l ) S k ( τ l ) - - - ( 1 )
通过删除[Mk -1l)Fkl)]中对应非测量点的行得到下面的方程:
Φkl)=Akl)Skl)                  (2)
因此对于m个波段即得到m个单一波段的矩阵方程;
(3)优化可行光源区域选择
I)在网格剖分T1上,将上述含有m个单一波段的矩阵方程组合成含有m个波段的一个方程:
根据光源谱信息的性质,每个谱段τl上所占整个谱段的能量表示为S(τl)=ω(τl)·S,其中ω(τl)≥0,
Figure G2008101168184D00041
ω(τl)可通过实验测定;S表示光源总的能量。在网格剖分T1上,考虑光源谱信息的性质,通过方程(2)得到:
AS1=Φ                        (3)
其中 Φ = Φ 1 ( τ 1 ) Φ 1 ( τ 2 ) . . . Φ 1 ( τ m ) , A = ω ( τ 1 ) A 1 ( τ 1 ) ω ( τ 2 ) A 1 ( τ 2 ) . . . ω ( τ m ) A 1 ( τ m )
Φ1l)和A1l)都由方程(2)计算得到,S1是网格剖分T1上的未知光源密度变量;在(3)式两边同乘AT,得到:
ATAS1=ATΦ                    (4)
AT是A的转置矩阵,AT·A是一个的对称矩阵,所以(4)式是一个标准线形方程。
II)在网格剖分T1上,确定优化光源可行区域:
具体方法为:在网格剖分T1上,假定整个被测物体Ω为光源可行区域,给定未知光源密度的初始值S1 0,根据公式(4),网格剖分T1上的未知光源密度变量S1用下面的关系式进行迭代求解:
S 1 n + 1 = S 1 n + α n β n
其中步长
Figure G2008101168184D00046
搜索方向
Figure G2008101168184D00047
n是迭代次数;
在迭代的过程中,用Φ1 measl)代替Φ1l),这是因为噪声在迭代的过程中产生的影响很小;Φ1 measl)是网格剖分T1上的边界测量值,由Φmeas(x,τl)通过最近邻域插值得到;当||βn||≤δ或n>Nmax,迭代停止,此时便得到了未知光源密度S1的优化解S*,优化解S*的大致区域为Ω*,这个区域Ω*我们称之为优化的光源可行区域;在每次迭代的过程中,如果光源的密度小于零,则将其置为零;δ取值范围为10-5-10-2,Nmax为最大迭代次数,取值不超过1000。
(4)光源重建
在网格序列{T1,...,Tk,...}上,根据选择好的优化光源可行区域Ω*,重新将含有m个单一波段的矩阵方程组合成含有m个波段的一个方程:
考虑未知光源密度Skl)与边界测量值Φk measl)之间的线性关系,由方程(1): Φ k ( τ l ) = M k - 1 ( τ l ) F k ( τ l ) S k ( τ l ) , 得到:
Φ k meas ( τ l ) = G k ( τ l ) S k ( τ l )
其中Gkl)通过移出[Mkl)-1Fkl)]中对应非测量点的行而得到;考虑到光源谱信息和优化可行光源区域Ω*,得到:
Φ k meas = G k W k S k
其中 G k = w ( &tau; 1 ) G k ( &tau; 1 ) w ( &tau; 2 ) G k ( &tau; 2 ) . . . w ( &tau; m ) G k ( &tau; m ) , &Phi; k = &Phi; k meas ( &tau; 1 ) &Phi; k meas ( &tau; 2 ) . . . &Phi; k meas ( &tau; m ) , W k = Diag ( w k ( 11 ) , w k ( 22 ) , . . . , w k ( ii ) , . . . ) , w k ( ii ) = 1 s k ( i ) &GreaterEqual; &gamma; k S k max 0 s k ( i ) < &gamma; k s k max
Gk表示多谱下的边界测量值与光源位置量之间的关系,Wk为k层网格上的对角矩阵,wk(ii)是对角矩阵Wk的第i个对角元素;Φk measl)是波段τl在网格剖分Tk上的边界测量值,由Φmeas(x,τl)通过最近邻域插值得到;当k=1时,sk(i)是通过迭代方法得到的优化解S*,sk max是优化解S*的最大值;当k≥2时,sk(i)和sk max是网格剖分Tk上的未知光源密度的初始值和初始值的最大值,通过插值网格剖分Tk-1上的重建结果获得,即
Figure G2008101168184D00055
Sk-1是在网格剖分Tk-1上重建得到的结果,Ik-1 k是插值因子,通过子单元继承父单元的重建结果值的方式来实现;γk是比例因子,是随k变化的分段常数;最终通过保留GkWk中对应可行光源区域的列,得到光源可行区域的未知光源密度变量Sk p与表面测量值Φk meas在网格剖分Tk上的关系: A k S k p = &Phi; k meas
II)利用Tikhonov正则化方法建立目标函数并优化从而得到重建结果:
上一步虽然得到了未知光源密度变量Sk p与表面测量值Φk meas之间的关系,由于Ak矩阵的病态,不能通过直接求解的方法得到重建结果,因此利用Tikhonov正则的方法,并考虑未知光源密度变量不能为负值的特性,得到在网格剖分Tk上的目标函数:
min 0 &le; s k p &le; s sup k &Theta; k ( S k p ) = { | | A k S k p - &Phi; k meas | | &Lambda; + &lambda; k &CenterDot; ( S k p - S k init ) T ( S k p - S k init ) }
考虑到多光谱以及非匀质模型会极大的增加计算量,因此采用谱梯度大规模优化算法对该目标函数进行优化;Ssup k是光源密度上界,经验取值为不大于1000,||V||Λ=VTΛV,λk正则化系数,经验取值10-9-10-6,Sk init是网格剖分Tk上的未知光源密度初始值,取值范围为10-7-10-4;在判断优化过程是否中止时,采用当前优化梯度
Figure G2008101168184D00058
以及优化迭代次数Nk te作为评判准则,即当
Figure G2008101168184D00059
Figure G2008101168184D000510
时,优化停止,得到重建结果,即光源可行区域的光源密度分布;上标te表示当前的迭代次数;εg是梯度阈值,其经验取值不超过10-6;Nk max是最大优化迭代次数,经验取值不超过104
III)当优化完成后,判断重建是否停止:
利用优化求得的光源密度求解任意一个谱段τl上的光密度分布Φk cl),并计算
Figure G2008101168184D00061
Figure G2008101168184D00062
或k≥Kmax,重建停止;Φk cl)是通过将优化求得的光源密度代入扩散方程而求得的边界处光子流密度,Φk measl)是波段τl在网格剖分Tk上边界测量点的测量值,由Φmeas(x,τl)通过最近邻域插值得到;εΦ经验取值范围10-9-10-7;Kmax是网格最大剖分次数;如果不能满足
Figure G2008101168184D00063
或k≥Kmax,根据重建的结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域;通过以上方法对可行和禁止光源区域内的网格单元的选择,对区域剖分网格进行自适应的细分;并从网格剖分Tk转到Tk+1,然后转到第(4)(I)步,重复前面的步骤,直至重建停止;
本发明基于单视图测量数据,利用多光谱提供的先验信息,通过迭代方法进行可行光源区域的选择,重建光源的不准确问题得到了很好的解决,重建的速度和质量由于光源可行区域的使用而被进一步提高,此方法对自发荧光断层成像的发展具有重要的应用价值。
附图说明:
图1.非均匀仿体
图中:(a)非均匀仿体,1代表肌肉、2骨骼、3心脏、4肺、5肝和6单光源;(b)重建算法用到的初始网格;(c)包含肌肉、骨骼、肺的横截面,光源在右肺中;箭头所示的方向为CCD探测的方向,7前视图方向,8左视图方向,9后视图方向,10右视图方向。
图2.程序主流程图。
图3.确立光源可行区域的子流程图。
图4.确定光源可行区域的算法示意图。
图5.重建结果示意图。
表1.各组织在不同波长下的光学特性参数。
具体实施方式:
下面结合附图详细描述本实施例。
(1)数据获取:
为了验证本方法,我们设计了非匀质模型来近似模仿小鼠腹部内的组织,如图1(a)所示,生物组织分别为:肌肉、骨骼、心脏、肺、肝,各个组织的光学特性参数如表1所示。在此实验中,我们将重建自发光源的光源密度分布作为重建目标,在重建实验中,这个光源大小为半径1毫米,光源密度为0.238nano-Watts/mm3,位置为(-3,5,15)。它的波谱范围从500纳米到750纳米,为了进行多光谱实验,我们依据当前的探测水平,把这一波谱范围分成了5个谱段进行分别探测,分别为τ1=[500,550)nm,τ2=[550,600)nm,τ3=[600,650)nm,τ4=[650,700)nm,τ5=[700,750)nm。边界测量数据由蒙特卡罗方法产生,并按照图1(c)所示的方法取前视图的边界数据,并求得每个波段的光子流密度Φmeas(x,τl);
(2)有限元离散化处理:
光子在生物组织中传输的数学模型用下面的扩散方程表示:
- &dtri; &CenterDot; ( D ( x , &tau; l ) &dtri; &Phi; ( x , &tau; l ) ) + &mu; a ( x , &tau; l ) &Phi; ( x , &tau; l ) = S ( x , &tau; l ) ( x &Element; &Omega; )
&Phi; ( x , &tau; l ) + 2 A ( x ; n , n &prime; ) D ( x , &tau; l ) ( v ( x ) &CenterDot; &dtri; &Phi; ( x , &tau; l ) ) = 0 ( x &Element; &PartialD; &Omega; )
其中Ω是被测物体;
Figure G2008101168184D00073
是被测物体的表面;Φ(x,τl)是光子流密度;S(x,τl)是光源密度;根据有限元理论,得到下面的弱解形式:
&Integral; &Omega; ( D ( x , &tau; l ) ( &dtri; &Phi; ( x , &tau; l ) ) &CenterDot; ( &dtri; &Psi; ( x , &tau; l ) ) + &mu; a ( x , &tau; l ) &Phi; ( x , &tau; l ) &Psi; ( x , &tau; l ) ) dx
+ &Integral; &PartialD; &Omega; 1 2 A n ( x ) &Phi; ( x , &tau; l ) &Psi; ( x , &tau; l ) dx = &Integral; &Omega; S ( x , &tau; l ) &Psi; ( x , &tau; l ) dx ( &ForAll; &Psi; ( x , &tau; l ) &Element; H 1 ( &Omega; ) )
H1(Ω)是Sobelev空间,Ψ(x,τl)是任意的试探函数。在自适应有限元分析的框架下,基于自适应的网格细分,把重建区域网格剖分表示为{T1,...,Tk,...}的网格序列,其中网格剖分Tk包括
Figure G2008101168184D00076
个离散单元,网格剖分T1包含6878个离散单元,如图1(b)所示。利用有限元方法,对弱解形式进行离散,在每个波段τl(l=1,2,3,4,5)上得到下面矩阵方程:
(Kkl)+Ckl)+Bkl))Φkl)=Fkl)Skl)
矩阵元素的具体形式为:
k ij k ( &tau; l ) = &Integral; &Omega; D ( x , &tau; l ) ( &dtri; &Psi; i k ( x , &tau; l ) ) &CenterDot; ( &dtri; &Psi; j k ( x , &tau; l ) ) dx c ij k ( &tau; l ) = &Integral; &Omega; &mu; a ( x , &tau; l ) &Psi; i k ( x , &tau; l ) &Psi; j k ( x , &tau; l ) dx b ij k ( &tau; l ) = &Integral; &PartialD; &Omega; &Psi; i k ( x , &tau; l ) &Psi; j k ( x , &tau; l ) / ( 2 A ( x ; n , n &prime; ) ) dx f ij k ( &tau; l ) = &Integral; &Omega; &Psi; i k ( x , &tau; l ) &Psi; j k ( x , &tau; l ) dx
令Kkl)+Ckl)+Bkl)=Mkl),Mkl)是稀疏正定矩阵,所以得到:
&Phi; k ( &tau; l ) = M k - 1 ( &tau; l ) F k ( &tau; l ) S k ( &tau; l ) (1)通过删除
[Mk -1l)Fkl)]中对应非测量点的行得到下面的方程:
Φkl)=Akl)Skl)                     (2)
对于5个波段得到5个单一波段的矩阵方程;
(3)优化可行光源区域选择
I)在网格剖分T1上,将上述含有5个单一波段的矩阵方程组合成含有5个波段的一个方程:
根据光源谱信息的性质,每个谱段τl上所占整个谱段的能量表示为
S(τl)=ω(τl)·S,ω(τl)≥0,
Figure G2008101168184D00082
其中S表示光源总的光源能量,值为1;ω(τl)通过实验测定分别为:ω(τ1)=0.222,ω(τ2)=0.167,ω(τ3)=0.222,ω(τ4)=0.167,ω(τ5)=0.222。
在网格剖分T1上,考虑光源谱信息的性质,通过方程(2)得到:
AS1=Φ                                     (3)
其中 &Phi; = &Phi; 1 ( &tau; 1 ) &Phi; 1 ( &tau; 2 ) &Phi; 1 ( &tau; 3 ) &Phi; 1 ( &tau; 4 ) &Phi; 1 ( &tau; 5 ) , A = 0.222 A 1 ( &tau; 1 ) 0.167 A 1 ( &tau; 2 ) 0.222 A 1 ( &tau; 3 ) 0.167 A 1 ( &tau; 4 ) 0.222 A 1 ( &tau; 5 )
Φ1l)和A1l)都由方程(2)计算得到,S1是网格剖分T1上的未知光源密度变量。在(3)式两边同乘AT,得到:
ATAS1=ATΦ                                  (4)
AT是A的转置矩阵,AT·A是一个34390×34390的对称矩阵,所以(4)式是一个标准线形方程。
III)在网格剖分T1上,确定优化光源可行区域:
具体方法为:在网格剖分T1上,假定整个被测物体Ω为光源可行区域,给定未知光源密度的初始值
Figure G2008101168184D00085
根据公式(4),网格剖分T1上的未知光源密度变量S1根据下面的关系式进行迭代求解:
S 1 n + 1 = S 1 n + &alpha; n &beta; n
其中步长
Figure G2008101168184D00087
搜索方向
Figure G2008101168184D00088
n是迭代次数。
在迭代的过程中,用Φ1 measl)代替Φ1l),Φ1 measl)是τl在网格剖分T1上的边界测量值,由Φmeas(x,τl)通过最近邻域插值得到;参数取值分别为:δ=5×10-3,Nmax=500。经过5次迭代后,||βn||=4.9e-5,||βn||≤δ,所以迭代停止,此时便得到了未知光源密度S1的优化解S*,优化解S*的大致区域为Ω*,这个区域Ω*我们称之为优化的光源可行区域。在每次迭代的过程中,如果光源的密度小于零,则将其置为零。
(4)光源重建
(I)在被测物体Ω的网格剖分{T1,...,Tk,...}的网格序列上,根据选择好的优化光源可行区域Ω*,重新将含有5个波段的矩阵方程组合成含有5个波段的一个方程:
考虑未知光源密度Skl)与边界测量值Φk measl)之间的线性关系,由方程(1): &Phi; k ( &tau; l ) = M k - 1 ( &tau; l ) F k ( &tau; l ) S k ( &tau; l ) , 得到:
&Phi; k meas ( &tau; l ) = G k ( &tau; l ) S k ( &tau; l )
其中Gkl)通过移出[Mkl)-1Fkl)]中对应非测量点的行而得到。考虑到光源谱信息和优化可行光源区域Ω*,得到:
&Phi; k meas = G k W k S k
其中 G k = w ( &tau; 1 ) G k ( &tau; 1 ) w ( &tau; 2 ) G k ( &tau; 2 ) w ( &tau; 3 ) G k ( &tau; 3 ) w ( &tau; 4 ) G k ( &tau; 4 ) w ( &tau; 5 ) G k ( &tau; 5 ) , &Phi; k = &Phi; k meas ( &tau; 1 ) &Phi; k meas ( &tau; 2 ) &Phi; k meas ( &tau; 3 ) &Phi; k meas ( &tau; 4 ) &Phi; k meas ( &tau; 5 ) , W k = Diag ( w k ( 11 ) , w k ( 22 ) , . . . , w k ( ii ) , . . . ) , w k ( ii ) = 1 s k ( i ) &GreaterEqual; &gamma; k S k max 0 s k ( i ) < &gamma; k s k max
Gk表示多谱下的边界测量值与光源位置量之间的关系,Wk为k层网格上的对角矩阵,Φk measl)是波段τl在网格剖分Tk上的边界测量值。当k=1时,sk(i)是通过迭代方法得到的优化解S*,sk max是优化解S*的最大值;当k≥2时,sk(i)和sk max是网格剖分Tk上的初始值和初始值的最大值,通过插值网格剖分Tk-1上的重建结果得到,即
Figure G2008101168184D00097
Sk-1是在网格剖分Tk-1上重建得到的结果,Ik-1 k是插值因子,通过子单元继承父单元的重建结果值的方式来实现;γk是比例因子,是随k变化的分段常数;k=1时,参数γk=0;k>1时,γk初始为0.1,并随着k增加,有γk+1=10γk。最终通过保留GkWk中对应可行光源区域的列,得到光源可行区域的未知光源密度变量Sk p与表面测量值Φk meas在网格剖分Tk上的关系:
A k S k p = &Phi; k meas
(II)利用Tikhonov正则化方法建立目标函数并优化从而得到重建结果:
上一步虽然得到了未知光源密度变量Sk p与表面测量值Φk meas之间的关系,由于Ak矩阵的病态,不能通过直接求解的方法得到重建结果,利用Tikhonov正则的方法,并考虑未知光源密度变量不能为负值的特性,得到在网格剖分Tk上的目标函数:
min 0 &le; s k p &le; s sup k &Theta; k ( S k p ) = { | | A k S k p - &Phi; k meas | | &Lambda; + &lambda; k &CenterDot; ( S k p - S k init ) T ( S k p - S k init ) }
考虑到多光谱以及非匀质模型会极大的增加计算量,因此采用谱梯度大规模优化算法对该目标函数进行优化。Ssup k是光源密度上界,取值为100,||V||Λ=VTΛV,λk正则化系数,取值1.0×10-8,Sk init是在网格剖分Tk上未知光源密度变量的初始值,取值1.0×10-5。在判断优化过程是否中止时,采用了当前优化梯度
Figure G2008101168184D00103
以及优化迭代次数Nk i作为评判准则,即当
Figure G2008101168184D00104
时,优化停止,得到重建结果,即光源可行区域的光源密度分布,其中εg是梯度阈值,取值1.0×10-8,Nk max是最大优化迭代次数,取值1.0×104
(III)当优化完成后,判断重建是否停止:
利用优化得到的光源密度分布求解任意一个谱段τl上的光子流密度Φk cl),并计算
Figure G2008101168184D00106
Figure G2008101168184D00107
Φ为正数)或k≥Kmax(Kmax是网格最大剖分次数),重建停止。Φk cl)是通过将优化求得的光源密度代入扩散方程而求得的边界处光子流密度,Φk measl)是波段τl在第k层网格上边界测量点的测量值;如果不能满足
Figure G2008101168184D00108
或k≥Kmax,利用重建的结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域。通过以上方法对可行和禁止光源区域内的网格单元的选择,并对区域剖分网格进行自适应的细分,从网格剖分Tk转到Tk+1,然后转到第(4)(I)步,重复前面的步骤,直至重建停止。
停止阈值εΦ和网格最大细分次数Kmax分别是1×10-8和3。当在网格剖分T1上重建完成后,计算在τ1=[500-550)nm的光密度分布Φk cl),求得
Figure G2008101168184D00109
此时
Figure G2008101168184D001010
(k=1)<(Kmax=3)所以重建未停止;利用网格剖分T1上的重建结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域。通过以上方法对可行和禁止光源区域内的网格单元的选择,并对区域剖分网格进行自适应的细分,并从网格剖分T1转到网格剖分T2,然后转到第(4)(I)步,重复前面的步骤;优化完成后计算在τ1=[500-550)nm的光密度分布Φk cl),求得
Figure G2008101168184D00111
因此重建停止,重建的结果如图5所示,重建的位置为(-1.53,5.3,14.98),此时重建的最大光源密度为0.224nano-Watts/mm3
表1
Figure G2008101168184D00112

Claims (1)

1.一种基于单视图的多光谱自发荧光断层成像重建方法,其特征在于,包括以下步骤:
1)数据获取:
在多光谱下,用滤波片将荧光波长λ离散为m个波段τ1,...,τm,其中τl=[λl-1,λl],l=1,2,...,m-1,τm=[λm-1,λm];利用CCD探测器在被测物体的表面
Figure F2008101168184C00011
对逃出的光子进行探测,获得每个波段τl的光流密度Q(x,τl),其中CCD探测到的表面为γ,x表示在被测物体中的位置;根据下面公式计算被测物体表面的光子流密度Φmeas(x,τl):
D(x,τl)=1/(3(μa(x,τl)+(1-g)μs(x,τl)))是生物组织的扩散系数,其中μa(x,τl)是生物组织的吸收系数,μs(x,τl)是生物组织散射系数,g是生物组织各向异性参数;v(x)是被测物体的表面
Figure F2008101168184C00013
的单位法向量;A(x;n,n′)≈(1+R(x))/(1-R(x)),n′为外界媒质的折射系数;当外界媒质是空气时,R(x)近似为R(x)≈-1.4399n-2+0.7099n-1+0.6681+0.0636n,其中n为空气的折射率;
2)有限元离散化处理:
光子在生物组织中传输的数学模型用下面的扩散方程表示:
- &dtri; &CenterDot; ( D ( x , &tau; l ) &dtri; &Phi; ( x , &tau; l ) ) + &mu; a ( x , &tau; l ) &Phi; ( x , &tau; l ) = S ( x , &tau; l ) ( x &Element; &Omega; )
&Phi; ( x , &tau; l ) + 2 A ( x ; n , n &prime; ) D ( x , &tau; l ) ( v ( x ) &CenterDot; &dtri; &Phi; ( x , &tau; l ) ) = 0 ( x &Element; &PartialD; &Omega; )
其中Ω是被测物体;
Figure F2008101168184C00016
是被测物体的表面;Φ(x,τl)是光子流密度;S(x,τl)是光源密度;根据有限元理论,得到下面的弱解形式:
&Integral; &Omega; ( D ( x , &tau; l ) ( &dtri; &Phi; ( x , &tau; l ) ) &CenterDot; ( &dtri; &Psi; ( x , &tau; l ) ) + &mu; a ( x , &tau; l ) &Phi; ( x , &tau; l ) &Psi; ( x , &tau; l ) ) dx
+ &Integral; &PartialD; &Omega; 1 2 A n ( x ) &Phi; ( x , &tau; l ) &Psi; ( x , &tau; l ) dx = &Integral; &Omega; S ( x , &tau; l ) &Psi; ( x , &tau; l ) dx ( &ForAll; &Psi; ( x , &tau; l ) &Element; H 1 ( &Omega; ) )
H1(Ω)是Sobelev空间,Ψ(x,τl)是任意给定的试探函数;在自适应有限元分析的框架下,基于自适应的网格细分,把重建区域网格剖分表示为{T1,...,Tk,...}的网格序列,其中网格剖分Tk包括
Figure F2008101168184C00019
个离散单元;利用有限元方法,对弱解形式进行离散,得到下面矩阵方程:
(Kkl)+Ckl)+Bkl))Φkl)=Fkl)Skl)
矩阵元素的具体形式为:
k ij k ( &tau; l ) = &Integral; &Omega; D ( x , &tau; l ) ( &dtri; &Psi; i k ( x , &tau; l ) ) &CenterDot; ( &dtri; &Psi; j k ( x , &tau; j ) ) dx c ij k ( &tau; l ) = &Integral; &Omega; &mu; a ( x , &tau; l ) &Psi; i k ( x , &tau; l ) &Psi; j k ( x , &tau; l ) dx b ij k ( &tau; l ) = &Integral; &PartialD; &Omega; &Psi; i k ( x , &tau; l ) &Psi; j k ( x , &tau; l ) / ( 2 A ( x , n , n &prime; ) ) dx f ij k ( &tau; l ) = &Integral; &Omega; &Psi; i k ( x , &tau; l ) &Psi; j k ( x , &tau; l ) dx
令Kkl)+Ckl)+Bkl)=Mkl),Mkl)是稀疏正定矩阵,所以得到:
&Phi; k ( &tau; l ) = M k - 1 ( &tau; l ) F k ( &tau; l ) S k ( &tau; l ) - - - ( 1 ) 通过删除
[Mk -1l)Fkl)]中对应非测量点的行得到下面的方程:
Φkl)=Akl)Skl)           (2)
因此对于m个波段即得到m个单一波段的矩阵方程;
3)优化可行光源区域选择
I)在网格剖分T1上,将上述含有m个单一波段的矩阵方程组合成含有m个波段的一个方程:
根据光源谱信息的性质,每个谱段τl上所占整个谱段的能量表示为S(τl)=ω(τl)·S,其中ω(τl)≥0,ω(τl)通过实验测定;S表示光源总的能量;在网格剖分T1上,考虑光源谱信息的性质,通过方程(2)得到:
AS1=Φ             (3)
其中 &Phi; = &Phi; 1 ( &tau; 1 ) &Phi; 1 ( &tau; 2 ) &CenterDot; &CenterDot; &CenterDot; &Phi; 1 ( &tau; m ) , A = &omega; ( &tau; 1 ) A 1 ( &tau; 1 ) &omega; ( &tau; 2 ) A 2 ( &tau; 2 ) &CenterDot; &CenterDot; &CenterDot; &omega; ( &tau; m ) A 1 ( &tau; m )
Φ1l)和A1l)都由方程(2)计算得到,S1是网格剖分τ1上的未知光源密度变量;在(3)式两边同乘AT,得到:
ATAS1=ATΦ            (4)
AT是A的转置矩阵,AT·A是一个
Figure F2008101168184C00026
的对称矩阵,所以(4)式是一个标准线形方程;
II)在网格剖分T1上,确定优化光源可行区域:
具体方法为:在网格剖分T1上,假定整个被测物体Ω为光源可行区域,给定未知光源密度的初始值S1 0,根据公式(4),网格剖分T1上的未知光源密度变量S1用下面的关系式进行迭代求解:
S 1 n + 1 = S 1 n + &alpha; n &beta; n
其中步长搜索方向
Figure F2008101168184C00032
n是迭代次数;
在迭代的过程中,用Φ1 measl)代替Φ1l),Φ1 measl)是网格剖分T1上的边界测量值,由Φmeas(x,τl)通过最近邻域插值得到;当||βn||≤δ或n>Nmax,迭代停止,此时便得到了未知光源密度S1的优化解S*,优化解S*的大致区域为Ω*,这个区域Ω*我们称之为优化的光源可行区域;在每次迭代的过程中,如果光源的密度小于零,则将其置为零;δ取值范围为10-5-10-2,Nmax为最大迭代次数,取值不超过1000;
4)光源重建
I)在网格序列{T1,...,Tk,...}上,根据选择好的优化光源可行区域Ω*,重新将含有m个单一波段的矩阵方程组合成含有m个波段的一个方程:
考虑未知光源密度Skl)与边界测量值Φk measl)之间的线性关系,由方程(1):
Figure F2008101168184C00033
得到:
&Phi; k meas ( &tau; l ) = G k ( &tau; l ) S k ( &tau; l )
其中Gk(τl)通过移出[Mkl)-1Fkl)]中对应非测量点的行而得到;考虑到光源谱信息和优化可行光源区域Ω*,得到:
&Phi; k meas = G k W k S k
其中 G k = w ( &tau; 1 ) G k ( &tau; 1 ) w ( &tau; 2 ) G k ( &tau; 2 ) &CenterDot; &CenterDot; &CenterDot; w ( &tau; m ) G k ( &tau; m ) , &Phi; k meas = &Phi; k meas ( &tau; 1 ) &Phi; k meas ( &tau; 2 ) &CenterDot; &CenterDot; &CenterDot; &Phi; k meas ( &tau; m ) , Wk=Diag(wk(11),wk(22),…wk(ii),…),
w k ( ii ) = 1 s k ( i ) &GreaterEqual; &gamma; k s k max 0 s k ( i ) < &gamma; k s k max
Gk表示多谱下的边界测量值与光源位置量之间的关系,Wk为k层网格上的对角矩阵,wk(ii)是对角矩阵Wk的第i个对角元素;Φk measl)是波段τl在网格剖分Tk上的边界测量值,由Φmeas(x,τl)通过最近邻域插值得到;当k=1时,sk(i)是通过迭代方法得到的优化解S*,sk max是优化解S*的最大值;当k≥2时,sk(i)和sk max是网格剖分Tk上的未知光源密度的初始值和初始值的最大值,通过插值网格剖分Tk-1上的重建结果获得,即
Figure F2008101168184C00039
Sk-1是在网格剖分Tk-1上重建得到的结果,Ik-1 k是插值因子,通过子单元继承父单元的重建结果值的方式来实现;γk是比例因子,是随k变化的分段常数;最终通过保留GkWk中对应可行光源区域的列,得到光源可行区域的未知光源密度变量Sk p与表面测量值Φk meas在网格剖分Tk上的关系:
A k S k p = &Phi; k meas
II)利用Tikhonov正则化方法建立目标函数并优化从而得到重建结果:
上一步虽然得到了未知光源密度变量Sk p与表面测量值Φk meas之间的关系,由于Ak矩阵的病态,不能通过直接求解的方法得到重建结果,因此利用Tikhonov正则的方法,并考虑未知光源密度变量不能为负值的特性,得到在网格剖分Tk上的目标函数:
min 0 &le; s k p &le; s sup k &Theta; k ( S k p ) = { | | A k S k p - &Phi; k meas | | &Lambda; + &lambda; k &CenterDot; ( S k p - S k init ) T ( S k p - S k init ) }
考虑到多光谱以及非匀质模型会极大的增加计算量,因此采用谱梯度大规模优化算法对该目标函数进行优化;Ssup k是光源密度上界,经验取值为不大于1000,||V||Λ=VTΛV,λk正则化系数,经验取值10-9-10-6,Sk init是网格剖分Tk上的未知光源密度初始值,取值范围为10-7-10-4;在判断优化过程是否中止时,采用当前优化梯度以及优化迭代次数Nk te作为评判准则,即当
Figure F2008101168184C00044
时,优化停止,得到重建结果,即光源可行区域的光源密度分布;上标te表示当前的迭代次数;εg是梯度阈值,其经验取值不超过10-6;Nk max是最大优化迭代次数,经验取值不超过104
III)当优化完成后,判断重建是否停止:
利用优化得到的光源密度求解任意一个谱段τl上的光密度分布Φk cl),并计算||Φk cl)-Φk measl)||,当或k≥Kmax,重建停止;Фk cl)是通过将优化求得的光源密度代入扩散方程而求得的边界处光子流密度,Φk measl)是波段τl在网格剖分Tk上边界测量点的测量值,由Φmeas(x,τl)通过最近邻域插值得到;εФ经验取值范围10-7-10-9;Kmax是网格最大剖分次数;如果不能满足
Figure F2008101168184C00047
或k≥Kmax,根据重建的结果分别用直接最大值选择方法和等级值亏修正基的误差估计方法选取可行和禁止光源区域;通过以上方法对可行和禁止光源区域内的网格单元的选择,对区域剖分网格进行自适应的细分;并从网格剖分Tk转到Tk+1,然后转到第4)(I)步,直至重建停止。
CN2008101168184A 2008-07-18 2008-07-18 基于单视图的多光谱自发荧光断层成像重建方法 Expired - Fee Related CN101342075B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN2008101168184A CN101342075B (zh) 2008-07-18 2008-07-18 基于单视图的多光谱自发荧光断层成像重建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN2008101168184A CN101342075B (zh) 2008-07-18 2008-07-18 基于单视图的多光谱自发荧光断层成像重建方法

Publications (2)

Publication Number Publication Date
CN101342075A CN101342075A (zh) 2009-01-14
CN101342075B true CN101342075B (zh) 2010-06-02

Family

ID=40244275

Family Applications (1)

Application Number Title Priority Date Filing Date
CN2008101168184A Expired - Fee Related CN101342075B (zh) 2008-07-18 2008-07-18 基于单视图的多光谱自发荧光断层成像重建方法

Country Status (1)

Country Link
CN (1) CN101342075B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012071682A1 (zh) * 2010-11-30 2012-06-07 中国科学院自动化研究所 基于特异性的多模态三维光学断层成像系统和方法

Families Citing this family (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101947103B (zh) * 2010-09-20 2012-03-28 西安电子科技大学 全光学生物发光断层成像方法
CN102034266B (zh) * 2010-11-30 2013-03-06 中国科学院自动化研究所 激发荧光断层成像的快速稀疏重建方法和设备
CN102393969B (zh) * 2011-06-02 2013-04-03 西安电子科技大学 基于生物组织特异性的光学三维成像方法
CN102871646B (zh) * 2012-08-16 2014-09-10 清华大学 一种大数据量荧光分子断层成像重建方法
CN102940482B (zh) * 2012-11-22 2015-03-25 中国科学院自动化研究所 自适应的荧光断层成像重建方法
CN103207946A (zh) * 2013-03-08 2013-07-17 西安交通大学 基于截断奇异值和全变分的闪光照相客体正则化重建方法
CN104921706B (zh) 2015-07-08 2018-02-09 北京工业大学 基于多任务贝叶斯压缩感知方法的生物发光断层成像重建算法
CN105326475B (zh) * 2015-09-16 2018-01-19 西北大学 一种基于多光源分辨的生物发光断层成像重建方法
CN105455780A (zh) * 2015-11-17 2016-04-06 西北大学 基于双网格的有限投影的荧光分子断层成像重建方法
CN105581779B (zh) * 2015-12-13 2018-07-31 北京工业大学 一种直接融合结构成像的生物发光断层成像重建的方法
CN105559750B (zh) * 2015-12-13 2018-06-01 北京工业大学 组织结构引导的复合正则化生物发光断层成像重建方法
CN107358653A (zh) * 2017-08-15 2017-11-17 北京数字精准医疗科技有限公司 成像重建方法及装置
CN107392977B (zh) * 2017-08-22 2021-04-13 西北大学 单视图切伦科夫发光断层成像重建方法
CN108814550A (zh) * 2018-04-16 2018-11-16 北京工业大学 一种基于神经网络的近红外光谱断层成像重建方法
CN110047085B (zh) * 2019-04-17 2023-03-31 泰山学院 一种针对肺ct图像阈值分割结果的肺膜粘连结节区域精确修复方法
CN112089434B (zh) * 2020-10-16 2024-05-03 陕西师范大学 一种多光谱生物发光断层成像方法和系统

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1953701A (zh) * 2004-03-11 2007-04-25 通用医院有限公司 使用荧光蛋白进行断层摄影成像的方法和系统
CN101221128A (zh) * 2007-04-18 2008-07-16 中国科学院自动化研究所 一种基于自适应有限元的多光谱重建方法
CN101292863A (zh) * 2007-04-25 2008-10-29 中国科学院自动化研究所 一种基于单谱或混合段测量的自适应有限元光源重建方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1953701A (zh) * 2004-03-11 2007-04-25 通用医院有限公司 使用荧光蛋白进行断层摄影成像的方法和系统
CN101221128A (zh) * 2007-04-18 2008-07-16 中国科学院自动化研究所 一种基于自适应有限元的多光谱重建方法
CN101292863A (zh) * 2007-04-25 2008-10-29 中国科学院自动化研究所 一种基于单谱或混合段测量的自适应有限元光源重建方法

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2012071682A1 (zh) * 2010-11-30 2012-06-07 中国科学院自动化研究所 基于特异性的多模态三维光学断层成像系统和方法

Also Published As

Publication number Publication date
CN101342075A (zh) 2009-01-14

Similar Documents

Publication Publication Date Title
CN101342075B (zh) 基于单视图的多光谱自发荧光断层成像重建方法
CN105894562B (zh) 一种光学投影断层成像中的吸收和散射系数重建方法
JP4271040B2 (ja) 実時間の光学的トモグラフィーの正規化された差方法の変更
CN103239255B (zh) 一种锥束x射线发光断层成像方法
Davis et al. Image-guided diffuse optical fluorescence tomography implemented with Laplacian-type regularization
CN101221128B (zh) 一种基于自适应有限元的多光谱重建方法
CN101856220B (zh) 定量光学分子断层成像装置和重建方法
CN109584323A (zh) 超声反射信息约束的腹部病变电阻抗图像重建方法
CN105559750B (zh) 组织结构引导的复合正则化生物发光断层成像重建方法
Cong et al. Multispectral bioluminescence tomography: methodology and simulation
CN111915733B (zh) 基于LeNet网络的三维锥束X射线发光断层成像方法
CN103300829B (zh) 一种基于迭代重加权的生物自发荧光断层成像方法
US20120049088A1 (en) Systems, methods and computer-accessible media for hyperspectral excitation-resolved fluorescence tomography
JP2003528291A (ja) 相対的検出器値を使用した散乱媒体の画像化
CN104921706B (zh) 基于多任务贝叶斯压缩感知方法的生物发光断层成像重建算法
CN102988026A (zh) 基于乘子法的自发荧光断层成像重建方法
CN107220961A (zh) 一种基于半阈值追踪算法的荧光分子断层成像重建方法
CN106097441A (zh) 基于l1范数与tv范数的复合正则化生物发光断层成像重建方法
US10775308B2 (en) Apparatus and methods for determining optical tissue properties
CN103393410A (zh) 一种基于交替迭代运算的荧光分子断层成像重建方法
CN107330968A (zh) 基于五层脑结构模型的ot导引脑功能半三维dot方法
Slavine et al. Semi-automated image processing for preclinical bioluminescent imaging
CN107374588B (zh) 一种基于同步聚类的多光源荧光分子断层成像重建方法
Saiko et al. Visibility of capillaries in turbid tissues: an analytical approach
Guggenheim et al. Development of a multi-view, multi-spectral bioluminescence tomography small animal imaging system

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
C17 Cessation of patent right
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20100602

Termination date: 20100718