专利名称：一种基于水平集检测目标体内光源分布的方法
技术领域：
本发明属于光学分子影像领域，涉及检测光源的定量信息和图像分析，尤其是一
种基于水平集的目标体内光源分布的检测方法。
背景技术：
光学分子影像是一种可以在细胞分子水平上实现生物组织内部生物物理过程进 行无创实时动态成像的一种技术。通过获取到的目标表面出射光强分布，检测在组织内部 的光源分布则能够对研究目标定量化研究，为光源提供三维的定位信息。
目前普遍采用的检测技术大多基于是基于梯度的传统方法，由于光学分子影像的 信号十分微弱，导致了采集的数据集噪声较大，加之该问题的固有病态性，造成在待检测过 程中对初始值的选择、可行域的选择、以及噪声的容忍度都有着很强的要求。这限制了光学 分子影像在实际中的应用。寻找一种稳定、鲁棒的检测方法一直是该领域的热点难点问题， 到目前为止仍然没有一种通用的技术很好的求解这个问题。 水平集是从界面传播等研究领域发展而来。水平集的灵活性使得能够应用于曲线 或者曲面动态演化和拓扑变化的情形。
发明内容
为了解决基于梯度的传统方法光学分子影像的信号微弱，导致了采集的数据集噪 声较大，使得检测过程中对初始值的选择、可行域的选择、以及噪声的容忍度要求高，限制 了光学分子影像应用的问题，本发明的目的是提供一种基于水平集对目标体光源分布进行 检测的方法。 为了实现上述目的，本发明基于水平集检测目标体内光源分布的方法，利用高性 能CCD相机在目标体表面捕捉到的出射光强分布，来检测目标体内部的光源，具体步骤如 下 步骤1 :利用高性能CCD相机在目标体表面捕捉到的出射光强分布； 步骤2 :使用水平集的分片恒定函数和分片常数值的联合迭代循环对目标体的光
源分布待检测区域的位置进行检测，由分片常数值控制得到光源分布每个子区域上的光强
值，得到目标体的光源分布的可视化定位和定量信息； 步骤3 :对目标体的光源分布的可视化定位和定量信息提供数据后处理。
本发明的有益效果是通过使用基于水平集的光源信息检测技术，由分片恒定函
数控制待检测目标体光源分布的位置，以及分片常数值控制在每个子区域上的光源值，通
过对分片恒定函数和分片常数值的联合迭代循环，最终检测得到目标体光源分布的定量和
定位信息。使用该方法，在检测过程中对噪声和初值不敏感，在不同初值、较强噪声以及较
大可行区域下均可得到可靠的目标体光源分布检测结果。还有，在子区域水平数未知的情
况下，仍然可以对目标体光源分布定量定位。此外，由于该水平集方法不要对水平集函数重
新初始化，提高了目标体光源分布进行检测的效率，拓展了光学分子影像的应用，通过获取到的目标体出射光强分布，检测在组织内部的光源分布则能够对研究目标定量化研究，为 光源提供三维的定位信息，由此本发明的方法具有很强的鲁棒性，可以有效地应用于成像 系统性能较低或者实验环境较差的情况，本发明可以推动和拓展在图像配准与分割，图形 学和仿真，计算机立体成像，隐式表面以及逆问题等领域的应用。
图1为本发明的整体过程图。 图2为本发明的检测方法流程图。 图3为仿真实验所使用的非匀质圆柱仿体。 图4为不同初值下使用本发明检测光源的位置与强度演化曲线以及与使用 Powell方法的对比。 图5为使用本发明技术检测光源过程中的收敛曲线。 图6为不同高斯噪声下使用本发明检测光源的位置与强度演化曲线以及与使用 Powell方法的对比。 图7为使用本发明在不同噪声下检测光源强度相对于信噪比的线性回归拟合曲 线。
具体实施例方式
下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。应指出的是，
所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。 本发明是一种基于水平集的对目标体光源分布信息进行检测的方法，通过分片恒
定函数控制目标体光源分布的位置，以及分片常数值控制在每个子区域上的光源值，通过
迭代循环，最终检测得到光源分布的定量和定位信息。在光学分子影像中，例如通过对小动
物转基因，使得小动物在病灶处产生荧光蛋白酶，进而通过静脉注射荧光底物，从而在该蛋
白酶的催化作用下发射出可见或近红外光，光在生物组织中历经多次散射。最终有一部分
到达小动物皮肤出射到自由空间中被检测器接收。说发明就是使用检测得到的出射光强、
小动物的几何结构信息以及各个组织的光学特性参数等信息，来检测在小动物体内的光源
信息。本发明能在不同初值，较强的噪声下在使用比较大的感兴趣区域的情况下准确检测
出光源的信息。 下面结合附图l，详细描述本发明的对目标体内光源分布进行检测的方法。
首先采用水平集作为光学分子影像中检测光源的技术手段，利用高性能CCD相机 在目标体表面捕捉到的出射光强分布，采用水平集检测目标体，形成以下水平集目标体函 数 L(小，c， A) =F((Kc)+ / 。AK((j))dr+ / fiK2 (小)dr/(2 ii )
其中，采用增益拉格朗日策略强化技术将增益算子/ fl AK(小)dr和惩罚算子 / 。K2((t)dr/(2iO结合，用于保证迭代循环的稳定性；小表示待检测目标体内光源分布 的位置和分布的分片恒定函数，c表示目标体强度值的分片常数值，A表示拉格朗日因子，
Q为待检测区域，P表示惩罚因子，水平约束K(小)表示为K(f^fl^ —0 ;F((Kc)
5表示为尸(#) = |构(-,c)一0>"
1+ .
A 2
|M>，C)||，，是检测光源信息的经典目标函数；A
为加权矩阵，e为正则化参数，M表示描述生物自发光成像的系统矩阵，q表示光源强度分 布，①m代表在目标体表面捕捉到的出射光强分布。 然后使用水平集方法对目标体内的光源进行检测，也即对上述目标函数对分片恒 定函数和相应的分片常数值进行迭代，检测得到光源的信息。 最后经过数据后处理对数据进行分析，最终将得到光源的定位和定量信息可视 化。 作为一种具体的实现方案，请参阅图2，将检测区域Q演化成几个不同的子区域 Qi， i = 1，2，... ，n，检测光源的迭代循环步骤如下 步骤21 :设定迭代步数k二O，初始化用以表示待检测光源分布的位置和分布的分 片恒定函数小°、用以表示光源强度值的分片常数值(3°，以及拉格朗日因子^和惩罚因子 步骤22 :计算第k+l/2次被检测光源的分片恒定函数4k+1/2 :
/+1/2 —A0f/^K/,cA，》);其中，ok，ck，入k分别为第k次迭代的分辨恒定函数，分 片常数值和拉格朗日因子，At表示虚拟时间间隔，F表示F((Kc); 步骤23:如果H^ k+1/2—) k||/||) kj|<(^PF((tk+1/2) 〉F(一)为真，跳到步骤 24，否则，一+1= 一+1/2，跳到步骤25;这里，|^ k+1/2—) k ||< 5表示前后两侧迭代
所得的分片恒定函数小取整后相对误差小于某一个值S时，而F(小k+"2) >F(4k)表示本 次迭代较上次迭代的目标函数值要大时；^表示取与小最接近的整数，I I I I表示2范数， S表示一个小于1的正数； 步骤24 :利用牛顿法计算水平集的演化多项式方程 小W-小k+"2+At入K' (ctk+1)+T2K(^k+1)K' ((K+1)=0 求得第k+l次目标体的光源分布的水平集分片恒定函数4k+1 ; 步骤25 :如果k+l > Y为真，跳到步骤26，否则，跳到步骤27。其中，Y为更新分
片常数值的最小迭代步数； 步骤26 :使用共轭梯度法更新k+l次分片常数值Ck+1 :
ck+1 = ck+ P kdk，〃"1 = —VF(/+1, c"1, /L" + ， 其中，分片常数值(^+1的搜索步长Pk使用最速下降法线性搜索得到
A =argminJP(, —/ ^尸(/ ;^),/^)，(115为1^1次迭代的共轭方向，^通过弗 莱彻_李威斯方法得到； 步骤27 :更新拉格朗日因子Ak+1 :入w = A k+K (小k+1) / P ;采用A k+1 =
入k+K(C/i！更新k+l次迭代中拉格朗日因子入k。
步骤28 :更新目标体的光源信息函数《(/+V+1) = ^>，+V,(^+1)，
6
其中，Vi为基于分片恒定函数(K在每个子区域的特征函数；所述子区域的特征 函数在整个待检测区域的不同子区域有不同的特征函数值
(J) = i in Qi， i = 1，2， ， n. 同时采用不同分片常数值Ci， i = 1，2， ... ， n代表在不同子区域内的光源强度 分布函数q(r) =Ci，rG Qi，其中r代表Q的位置。所述光源强度分布函数的表达式
为WO)-Sc,^W)，其中，水平集特征函数1^按下列公式进行^= J! W-力/"i，
n w ，其中cii表示第i个子区域uTi上的特征常数。 步骤29 :根据迭代步数判断k+1 = N，若为假，跳到步骤22继续迭代；若为真，检测 停止，其中，N为检测光源信息的最大迭代步数。 采用算子分解技术对于分片恒定函数的求解，根据水平集随虚拟时间演化的函数 A +/30 = 0 ，在稳态境况下小t = 0， = 0 ，在k+1次迭代上将其分解为以下两个方
程 (/+1/2—+a^/^(,y,A" = 0 (cj)k+1_(j)k+1/2)/At+AK' (ctk+1)+K(^k+1)K' (cj)k+1)/ii=0
其中，A t表示水平集演化时间间隔。 在本发明中，使用一种水平集分片恒定函数控制待检测光源的位置，以及分片常 数值控制在每个子区域上的光源值，通过对两者的迭代循环，最终检测得到光源分布的定 位和定量信息。与传统的基于梯度的方法相比较，该发明的方法对于给定不同的初值均可 以检测得到准确的结果，而且对噪声有更大的容忍性。此外，在可获得的先验信息较少时， 特别是对于感兴趣区域较大时，较之原有的方法，本发明能够检测得到可靠的光源信息。
请参阅图3，使用计算机，在Tecplot软件中显示仿真实验中所使用的非匀质圆柱 仿体。该仿体由类似于小动物器官组织的各部分组成，分别模拟生物组织的心脏、脾脏、肺 部(两个)、骨骼和肌肉的非匀质特性。各个部分都提供了类似于小动物体内各器官的先验 光学特性参数。在(-5， 5， 15)和(-5， -5， 15)坐标处分别放置了一个半径为1.8mm的球形 光源，它的光强为0.238纳瓦/毫米3。通过检测在椭圆柱表面上得的出射光强(单位为纳 瓦/毫米2)，加之前文所述的其它信息，来检测光源的信息。 请参阅图4为不同初值下使用本发明检测光源的位置与强度演化曲线以及与使 用Powell方法的对比。图4(a)-(c)为在初值分别为(a)0. 8，纳瓦/毫米3， (b) 1. 0，0. 0纳 瓦/毫米3和(c)1.5，0.0纳瓦/毫米3时，检测光源的结果和真实光源位置的图示。以上 图示由检测光源的强度函数的半幅全宽的等值面生成。图4(d)为使用Powell方法在与本 发明的方法一样感兴趣区域条件下得到的检测结果。在以上四个图中左上角的角度为以仿 体正面为基准，转过的角度(逆时针为正)。图4(a)-(c)与图4(d)比较，前三个实验在不 同初值下均检测出了准确的光源位置，而(d)只是得到散乱的碎片，据此很难判断出光源 的实际位置。 图4(e)-(g)为三组初值下的对应检测光源强度的演化曲线。检测的强度用实线 表示，实际的强度用虚线表示。在三组初值下，均在一定迭代次数的演化后，逐渐接近并且 收敛与真实强度。以上实验均在1.0%的高斯噪声条件下进行的。
请参阅图5为使用本发明检测光源过程中的收敛曲线。该曲线是在初值为0. 8，
0. 0纳瓦/毫米3时检测过程中目标函数F(q)随着迭代步数的变化。在整个过程中，能量 值单调下降，并最终趋于平缓，能量值不再下降。图中纵坐标为对原有目标函数值取10的 对数。对于其它情况下(不同噪声，不同初值)，亦有类似的收敛曲线。 请参阅图6为不同高斯噪声下使用本发明检测光源的位置与强度演化曲线以及 与使用Powell方法的对比。图6(a)-(d)为在高斯噪声分别为(a)0.25%， (b)6.25%， (c) 10% ， (d) 25%时，在初值为1. 0， 0. 0纳瓦/毫米3下，检测光源的结果和真实光源位置的 图示。以上图示由检测光源的强度函数的半幅全宽的等值面生成。图6(e)为使用Powell 方法在与本发明的方法一样感兴趣区域条件下并且噪声水平为25%时得到的结果。在以上 四个图中左上角的角度为以仿体正面为基准，转过的角度(逆时针为正)。图6(a)-(d)与 图6(e)比较，前四个实验在不同初值下均检测出了准确的光源位置，而(e)只是得到散乱 的碎片，据此很难判断出光源的实际位置。 图6(f)-(i)为四个噪声下的对应检测光源强度的演化曲线。检测的强度用实线 表示，实际的强度用虚线表示。在四组噪声水平下，均在一定迭代次数的演化后，逐渐接近 并且收敛与真实强度，其误差均在20%以内。因此，该发明可以有效地适用于成像系统性能 较低或者实验环境较差的情况。 请参阅图7为使用本发明在不同噪声下检测光源强度相对于信噪比的线性回归 拟合曲线。其中，五次实验的信噪比(噪声水平)分别为26.02分贝(0. 25%)，20.00分贝 (1.0% )，12.04分贝(6.25% )，10.00分贝(腦)，6.02分贝(25%)。光强的相对误差 与信噪比的线性回归表达式为y = 17. 993-0. 275x，其方差为0. 433.以上实验均在初值为
1. O，O. 0纳瓦/毫米3条件下进行。 以上所述，仅为本发明中的具体实施方式
，但本发明的保护范围并不局限于此，任 何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在 本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。
权利要求
一种基于水平集检测目标体内光源分布的方法，其特征在于，包括步骤如下步骤1利用高性能CCD相机在目标体表面捕捉到的出射光强分布；步骤2使用水平集的分片恒定函数和分片常数值的联合迭代循环对目标体的光源分布待检测区域的位置进行检测，由分片常数值控制得到光源分布每个子区域上的光强值，得到目标体的光源分布的可视化定位和定量信息；步骤3对目标体的光源分布的可视化定位和定量信息提供数据后处理。
2. 根据权利要求1所述检测目标体内光源分布的方法，其特征在于，还包括建立水平 集目标体函数L(小，c， A) =F((Kc)+ / 。AK(cj))dr+ / fiK2 (小)dr/(2 y )采用增益拉格朗日策略强化技术将增益算子/ 。 AK(小)dr和惩罚算子/ 。K2((t)dr/ (2iO结合，用于保证迭代循环的稳定性；小表示待检测目标体的光源分布的位置和分布 的分片恒定函数，c表示光源强度值的分片常数值，A表示拉格朗日因子，Q为待检测区域，p表示惩罚因子，水平约束K(小)表示为《W)-flW-01，2......n， F(小，A+|/%(ML '是检测光源信息的经典目标函数；c)表示为F^,c)^^lM《^,c)-^A为加权矩阵，13为正则化参数，M表示描述生物自发光成像的系统矩阵，q表示待检测的 光源强度分布，。m代表在目标体表面捕捉到的出射光强分布。
3.根据权利要求1所述检测目标体内光源分布的方法，其特征在于，所述迭代循环包 括步骤如下步骤21 :设定迭代步数k = O，初始化用以表示目标体的光源分布的位置和分布的分 片恒定函数小°、用以表示光源强度值的分片常数值(3°，以及拉格朗日因子^和惩罚因子步骤22 :计算第k+1/2次目标体的光源分布的分片恒定函数4k+1/2 :/""-^-△^/^(^Z,^)^*, (J3k，Ck，入k分别为第k次迭代的分辨恒定函数，分 片常数值和拉格朗日因子，At表示虚拟时间间隔，F表示F((Kc);步骤23:如果ll一 k+1/2—^〈柳F(C >F(<)为真，跳到步骤24，否则，一+1 = 一+1/2，跳到步骤25 ;这里，ll- k+1/2—^ k ||/||- k ||<^表示前后两侧迭代所得的分片恒定函数小取整后相对误差小于某一个值S时，而F(小k+"2) >F(4k)表示本次迭代较上次迭代的目标函数值要大时^表示取与小最接近的整数，表示2范数，S表示一个小于1的正数；步骤24 :利用牛顿法计算水平集的演化多项式方程 cj)k+1_(j)k+1/2+AtAK' (ctk+1)+T2K(^k+1)K' (c^k+1)=0 求得第k+l次目标体的光源分布的水平集分片恒定函数4k+1 ;步骤25:如果k+l〉 y为真，跳到步骤26，否则，跳到步骤27，其中，y为更新分片常 数值的最小迭代步数；步骤26 :使用共轭梯度法更新k+l次分片常数值Ck+1 : ck+i = ck+ p dk，<formula>formula see original document page 3</formula>其中，分片常数值Ck+1的搜索步长P k使用最速下降法线性搜索得到A = arg min , 一 — p VF(/+1 ， ), f) ， dk为k+1次迭代的共轭方向，4 k通过弗 莱彻_李威斯方法得到；步骤27:更新拉格朗日因子Ak+1:Ak+1= P+K(C/i！; 步骤28 :更新目标体的光源分布信息函数<formula>formula see original document page 3</formula>其中，^为基于分片恒定函数小，在每个子区域的特征函数；步骤29 :根据迭代步数判断k+1 = N，若为假，跳到步骤22继续迭代；若为真，检测停止，其中，N为检测目标体的光源分布信息的最大迭代步数。
4. 根据权利要求3所述检测目标体内光源分布的方法，其特征在于采用算子分解技术对于分片恒定函数的求解，根据水平集随虚拟时间演化的函数《+ = o ，在稳态境况下4t = 0，^:/^^^0，在k+l次迭代上将其分解为以下两个方程<formula>formula see original document page 3</formula> 其中，A t表示水平集演化时间间隔。
5. 根据权利要求3所述的检测方法，其特征在于采用Ak+1=入k+K(C/i！更新k+l次迭代中拉格朗日因子入k。
6. 根据权利要求3所述检测目标体内光源分布的方法，其特征在于所述子区域的特 征函数在整个待检测区域的不同子区域有不同的特征函数值(J) = i in Q i， i = 1， 2， . . . ， n.同时采用不同分片常数值Ci， i = 1，2， . . . ， n代表在不同子区域内的光源强度分布函 数q(r) = Ci，r G Qi，其中r代表Q的位置。
7. 根据权利要求6所述检测目标体内光源分布的方法，其特征在于所述光源强度分布函数的表达式为^>) = 2]Ci^(^)，其中，水平集特征函数L按下列公式进行'■=1K= ll W —(z'—"，其中cii表示第i个子区域^上的特征常数。
全文摘要
本发明涉及一种基于水平集检测目标体内光源分布的方法，利用高性能CCD相机在目标体表面捕捉到的出射光强分布，使用波尔兹曼传输方程的一阶球谐近似的扩散方程的前向模型，采用水平集作为检测光源的方法，由其分片恒定函数控制待检测光源分布的位置和分布，以及分片常数值控制在每个子区域上的光源值，同时融合算子分解技术和增益拉格朗日技术来加强该检测技术的效率，以对光源进行定位定量化，以及对光源的分布区域进行检测。该方法在检测过程中对噪声和初值不敏感，而且，在水平集的水平数未知时仍可检测出可靠的结果，具有很强的鲁棒性。该技术可以有效地应用于成像系统性能较低或者实验环境较差的情况。
文档编号A61B5/00GK101744610SQ20091009163
公开日2010年6月23日 申请日期2009年8月26日 优先权日2009年8月26日
发明者刘凯, 杨鑫, 田捷, 秦城虎 申请人:中国科学院自动化研究所