CN1427690A - 混浊媒体中的光学计算断层摄影术 - Google Patents

Abstract

光子迁移法被用来使埋在诸如组织之类的混浊媒体中的吸收物体成像。对于改善的分辨率，早到达的光子被检测以便以光学计算断层摄影术(CT)为基础提供图像重建数据。CT法是广义的，以便考虑光子路径的分布。点展宽函数(PSF)是按照用于迁移方程的格林函数表达的。这个PSF提供在广义级数展开法中使用的权重函数。用散射和吸收性质测量混浊媒体包括在双投影中收集到的同轴传输扫描。与多重散射相关联的图像模糊被消除，高分辨率图像能够被获得。

Description

混浊媒体中的光学计算断层摄影术

相关的专利申请

这份申请享受于2000年5月5日申请的美国专利申请第60/201,938号的利益。上述申请的全部教导在此通过引证全部并入本文。政府的支持

本发明通过专用拨款第P41-RR02594号全部或部份地得到来自国立健康研究所的资助。政府在本发明中有某种权利。

本发明的现有技术

X射线计算断层摄影术(CT)在人体内部结构成像方面已经是非常成功的。它已提供一种供临床使用的准确的、微分辨的和实时的医学成像工具。然而，使用X射线有一些缺点。反差对于某些类型的肿瘤(例如早期乳癌)是低的。X射线乳房照相术的误诊率是高的，而且X射线是诱导有机体突变的物质。近几年，基于光量子的成像技术已显著地引人关注。在被称为治疗窗的700和900nm之间的光谱区域中，光子不引起诱导有机体突变的效应，而且由于在这些波长下光线的吸收弱，所以它们能够深深地穿透到组织之中。对光学反差的敏感性是高的，而且能够获得光谱信息。再者，反差能够通过注射瞄准肿瘤细胞的外源性染剂得到提高。因此，由光量子提供的信息能够补足X射线CT的信息，也许还能提供用来检测身体内部的肿瘤和其它反常的替代的诊断工具。

本发明的概述

光学成像的主要障碍是人体组织非常混浊。注入组织的光子在它们被检测之前遭受多重散射事件。为了提取空间信息，不得不摆脱散射的影响。为了处理这个问题，已经开发了各种不同的光学成像方法。这些方法各有各的优点和缺点。利用弹道光子的成像技术对于深入组织内部1毫米的成像是非常好的，但是在深层组织成像的状况下不工作。为了去混浊度卷积和提高空间分辨率，一些图像重建方案已被设计。有限元素法和有限差分法已分别在时域和频域中得到发展。滤波逆投影已被应用于CW成像。首先在X射线CT中建议的用于个别voxel对测量的贡献的权重因数也已经使用Monte Carlo模拟计算出来并且被用在反演模型中。本发明涉及级数展开法对光学状况的使用。使用提前检测法，组织的三维空间图像能够通过考虑混浊度的影响而被重建。

问题可以通过比较光量子和X-光线在人体组织中的传播而被理解。如同众所周知的那样，穿越身体的X-射线沿着直线传播(图1a)。反之，由于散射，光量子的传播具有三维空间的传播。光量子路径的分布可以被想象成连接光源和检测器的管子(图1b)。管子的横截面的宽度按照到达光子被收集的时间变化。

光子路径分布的一些特征可以被注意：

首先，在飞行时间限制方面这种管子减少到一条直线，而且问题被减少到传统的X射线CT的问题。然而，这些所谓的弹道光子所生产的信号极其微弱，而且就厚的组织而言本质上是检测不到的；其次，对于早期检测时间[在飞行时间之后数百微微秒(ps)]，管子仍然非常狭窄，而且光子路径被很好地定义。被传输的信号是重要的，而且空间信息仍然得到很好的保护；第三，在后期检测时间[几纳秒(ns)]，管子能够完全充满感兴趣的器官，而空间分辨率被大大降低。这种状况被称为扩散极限。

通过用早期时间光子路径分布的狭窄管子代替X-光线的直线路径，光学CT程序可以被使用，即在X射线CT中使用的那个程序的修改。

在光学状况下，早期到达光子类似于X射线光子。与高度扩散的光子相比它们遭受较少的散射事件，因此保持大量的空间信息。就早期到达的光子而言信号水平仍然可能比较高。使用早期时间检测进行的测量与用连续波(CW)和频域技术获得的那些相比具有比较高的分辨率。清晰的图像可以使用光子路径密度的概念予以重建。众所周知，扩散近似解没有很好地描述早期到达的光子。直到t-t₀≈600ps时间窗用点展宽函数预测扩散近似值不适合正确的图像重建。考虑因果关系的点展宽函数(PSF)产生好得多的结果。用于早期时间的点展宽函数的正确形式对于图像重建起重要的作用。

本发明的图像重建方法是以在散射占主导地位而且光源和检测器之间的光子路径分布必须予以考虑的光学状况下使用级数展开方法为基础的。为了完成这个，PSF被用来产生权重函数矩阵。先前，已经使用扩散近似、微观的Beer-Lambert定律和Monte Carlo模拟对权重函数进行过讨论和计算。然而，PSF的使用提供深入选择权重函数的指导。根据PSF物理解释是比较清楚的，而且关于光子迁移的不同理论能够被检验，因为它们预测不同的PSF。此外，如同下面在方程(12)中所展示的那样，实际上有二组权重函数，一组用于散射反差，另一组用于吸收反差。例如，乳房组织中的肿瘤呈现吸收和散射两种反差。光子迁移曲线的早期部分对散射反差比对吸收反差更敏感。

光学CT的分辨率受一些因素限制，例如，散射效应和重建程序中不确定的性质。此外，投影和测量的总数可以增加，扇形光束的几何形状可以被用来改善数据收集效率。纤维光学系统可以在诊察中用于传输和/或收集来自患者的组织的光。使用重建程序用计算机实现的精密的时域光子迁移仪器能够提供高质量的乳房、脑和体内其它地方的光学CT图像。

附图简要说明

本发明的上述和其它的目的、特征和优势从下面在附图中图解说明的本发明的优选实施方案的更具体的描述将是显而易见的，其中相同的参考符号在不同的视图中始终表示同一零部件。这些附图不必按比例绘制，而是把重点放在图解说明本发明的原则上。

图1A和1B是把代数的重建技术用于光学CT的示意图，其中正在研究的样品被分成N×N×N元音，而吸收分布是用每个voxel之内的平均吸收表示的。

图2A和2B图解说明每个voxel被赋予一个包括用于X射线的权重因素，其中轨道上的voxel有100％贡献，而脱离轨道者有0％贡献，而且就光量子而言，权重因素是由光子路径分布分别确定的。

图3A和3B是依照本发明用来完成光学计算断层摄影的系统示意图。

图4图解说明散射媒体的尺寸和扫描几何学的实例。

图5A-5D是强度的等值线图：1-物体(a)-(b)和2-物体(c)-(d)；时间窗：(a)-(b)Δt＝534ps，宽度50ps；(c)-(d)Δt＝607ps，宽度50ps。

图6A-6D是点展宽函数，其中光源和检测器分别在(0，-3.2，0)厘米和(0，3.2，0)厘米，而且四个组合是给定的：(A)用扩散近似计算的z＝0平面的侧视图；(B)用因果关系修正计算的Z＝0平面的侧视图；(C)用扩散近似计算的y＝0平面的俯视图；(D)因果关系修正计算的y＝0平面的俯视图。

图7A-7D是有一个植入物体的重建图像，包括(A)用直接的X射线算法重建；(B)用扩散近似重建；(C)用因果关系修正重建；(D)精确的结构。

图8A-8D是两个植入物体的图像，包括(A)用直接的X射线算法重建；(B)用扩散近似重建；(C)与因果关系修正重建；(D)精确的结构。用直接的X射线算法重建不分辨两个植入物体。另一方面，用扩散近似重建去卷积做得过分而且导致比较小的图像。比较小的物体从3D视图中是看不到的。

图9A-9D展示图8C中的重建图像的四个切片的等值线图，其中这些切片是(a)Z＝10mm；(b)Z＝-6mm；(c)Z＝2mm；和(d)Z＝10mm。

图10图解说明本发明的优选实施方案的处理顺序。

本发明的详细说明

穿过人体的X射线强度的衰减的表达被很好地建立起来。假设正在研究的人体组织对于单色的X射线具有衰减分布μ_a ^x(r’)。然后，越过身体传输的X射线强度是： $I\left(r\right)=I_0\exp [-{\∫}_{r_0}^r\mathrm{dl}{\mathrm{\μ}}_a^x\left(r^{\′}\right)]-----\left(1\right)$

其中积分沿着连接在r₀处的光源和在r处的检测器的直线。方程(1)可以被重新写成： $-[\ln l\left(r\right)-\mathrm{In}{I}_0]=\underset{r_0}{\overset{r}{\∫}}\mathrm{dl}{\mathrm{\μ}}_a^x\left(r^1\right)----\left(2\right)$

上述的线积分往往被称为Radon变换。

不同于X-光线，近红外线在人体组织中遭受强烈的散射而且不遵循直线路径。光量子在它们被检测或被组织吸收之前遭受多重散射事件。在均匀一致的散射媒体中光子路径的分布已经采用各种不同的方法进行研究。业已确定，光子路径的分布在早期检测时间是狭窄的，但是随着时间增加而展开。

一般地说，肿瘤的出现产生光学的不均匀性，这种不均匀性是作为组织的散射和吸收性质的局部变化出现的。对于早期阶段的肿瘤，这些变化可以被看作是小的扰动。

近红外线光子在人体组织中传播是用迁移方程很好地描述的。对于散射分布为μ_s(r)而吸收分布为μ_a(r)的媒体，辐射强度L(r，，t)满足： $\frac{1}{c}\frac{\∂L(r,\hat{s},t)}{\∂t}+\hat{s}\·\▿L(r,\hat{s},t)=-[{\mathrm{\μ}}_a\left(r\right)+{\mathrm{\μ}}_s\left(r\right)]L(r,\hat{s},t)-----\left(3\right)$ $+{\mathrm{\μ}}_s\left(r\right)\underset{4}{\∫}\underset{\mathrm{\π}}{\∫}L(r,\hat{s},t)f_r(\hat{s}\·{\hat{s}}^{\′})d{\mathrm{\Ω}}^{\′}+Q(r,\hat{s},t)$

而通常称之为相位函数的标准化的微分散射截面f_Γ(·＇)满足： $\underset{4}{\∫}\underset{\mathrm{\π}}{\∫}{f}_r(\hat{s}\·{\hat{s}}^{\′})d{\mathrm{\Ω}}^{\′}=1------\left(4\right)$

在方程(3)中，c是在媒体中的光速，而Q(r，，t)是光源项。对于扰动系统，吸收、散射和相位函数的分布有小的形式变化：其中 $\mathrm{\δ}{\widetilde{\mathrm{\μ}}}_s(r,\hat{s}\·{\hat{s}}^{\′})={\mathrm{\δ\μ}}_s\left(r\right)f(\hat{s}\·{\hat{s}}^{\′})+{\mathrm{\μ}}_s[f_r(\hat{s}\·{\hat{s}}^{\′})-f(\hat{s}\·{\hat{s}}^{\′})]$

在方程(5)中，μ_a和μ_s分别是平均的吸收系数和散射系数；δμ_a(r)(＜＜μ_a)和δM_s(r)(＜＜μ_s)  是由肿瘤引起的扰动。一般地说，来自全局相位函数f(·＇)的相位函数f_r(·＇)的局部变化对于某些角度可能不小，但是方程(4)要求这样的变化在所有的立体角上积分之后彼此抵消。因此，把 (r，·＇)处理成小扰动。

在方程(5)给定的变化下，方程(3)能够能用扰动理论求解。对于点光源的情况， $Q(r,\hat{s},t)=\frac{1}{4\mathrm{\π}}\mathrm{\δ}(r-r_0)\mathrm{\δ}(t-t_0),------\left(6\right)$

方程(3)是用于有展开方程(7)的格林函数G(r，t|r₀，t₀)的方程式， $G_{\hat{s}}(r,t|r_0,t_0)=G_{\hat{s}}^{\left(0\right)}(r,t^0|r_0,t_0)+G_{\hat{s}}^{\left(1\right)}(r,t^0|r_0,t_0)+...----\left(7\right)$

其中G_s ⁽⁰⁾(r，t|r₀，t₀)是没有扰动时正常的格林函数，而G_s ⁽¹⁾(r，t|r₀，t₀)是由扰动造成的一阶修正。把方程(5)-(7)代入方程(3)而且仅仅保留到一阶项，我们有： $\frac{1}{c}\frac{\∂G_{\hat{s}}^{\left(0\right)}(r,t|r_0{t}_0)}{\∂t}+\hat{s}\·\▿G_{\hat{s}}^{\left(0\right)}(r,t|r_0{t}_0)+({\mathrm{\μ}}_a+{\mathrm{\μ}}_s)G_s^{\left(0\right)}(r,t|r_0{t}_0)-----\left(8\right)$ $-{\mathrm{\μ}}_s\underset{4}{\∫}\underset{\mathrm{\π}}{\∫}{G}_{\hat{s}}^{\left(0\right)}(r,t|r_0,t_0)f(\hat{s}\·{\hat{s}}^{\′})d{\mathrm{\Ω}}^{\′}=\frac{1}{4\mathrm{\π}}\mathrm{\δ}(r-r_0)\mathrm{\δ}(t-t_0)$ 和 $\frac{1}{c}\frac{\∂G_{\hat{s}}^{\left(1\right)}(r,t|r_0{t}_0)}{\∂t}+\hat{s}\·\▿G_{\hat{s}}^{\left(1\right)}(r,t|r_0{t}_0)+({\mathrm{\μ}}_a+{\mathrm{\μ}}_s)G_s^{\left(1\right)}(r,t|r_0{t}_0)-----\left(9\right)$ $-{\mathrm{\μ}}_s\underset{4}{\∫}\underset{\mathrm{\π}}{\∫}{G}_{\hat{s}}^{\left(1\right)}(r,t|r_0,t_0)f(\hat{s}\·{\hat{s}}^{\′})d{\mathrm{\Ω}}^{\′}$ $=-[\mathrm{\δ\μ}\left(a\right)\left(r\right)+\mathrm{\δ\μ}\left(r\right)]G_{\hat{s}}^{\left(0\right)}(r,t|r_0,t_0)+\underset{4}{\∫}\underset{\mathrm{\π}}{\∫}{G}_{\hat{s}}^{\left(0\right)}(r,t|r_0,t_0)\mathrm{\δ}{\widetilde{\mathrm{\μ}}}_s(r,\hat{s}\·{\hat{s}}^{\′})d{\mathrm{\Ω}}^{\′}$

方程(9)可以通过方程(8)的伴随方程求解，该伴随方程是： $\frac{1}{c}\frac{\∂G_{\hat{s}}^{\left(0\right)}(r,t|r_0{t}_0)}{\∂t_0}-\hat{s}\·{\▿}_0{G}_{\hat{s}}^{\left(0\right)}(r,t|r_0{t}_0)+({\mathrm{\μ}}_a+{\mathrm{\μ}}_s)G_{\hat{s}}^{\left(0\right)}(r,t|r_0{t}_0)-----\left(10\right)$ $-{\mathrm{\μ}}_s\underset{4}{\∫}\underset{\mathrm{\π}}{\∫}{G}_{\stackrel{^\′}{s}}^{\left(0\right)}(r,t|r_0,t_0)f(\hat{s}\·{\hat{s}}^{\′})d{\mathrm{\Ω}}^{\′}=\frac{1}{4\mathrm{\π}}\mathrm{\δ}(r-r_0)\mathrm{\δ}(t-t_0).$

显然，方程(9)的解能满足： $\frac{1}{4\mathrm{\π}}\underset{4}{\∫}\underset{\mathrm{\π}}{\∫}{\mathrm{d\ΩG}}_{\hat{s}}^{\left(0\right)}(r,t|r_0,t_0)={\∫}_{t_0}^{{t_1}^{+}}-{\mathrm{dt}}^{\′}\∫d^3{r}^{\′}[{\mathrm{\δ\μ}}_a\left(r^{\′}\right)+{\mathrm{\δ\μ}}_s\left(r^{\′}\right)]\underset{4}{\∫}\underset{\mathrm{\π}}{\∫}d{G}_{\hat{s}}^{\left(0\right)}(r,t|r_0,t_0)G_{\hat{s}}^{\left(0\right)}(r^{\′},t^{\′}|r_0,t_0)---\left(11\right)$ $+{\∫}_{t_0}^{t_1^{+}}{\mathrm{dt}}^{\′}\∫d^3{r}^{\′}\underset{4}{\∫}\underset{\mathrm{\π}}{\∫}d{\mathrm{\Ω}}^{\′}\underset{4}{\∫}\underset{\mathrm{\π}}{\∫}{\mathrm{d\ΩG}}_{\hat{s}}^{\left(0\right)}(r,t|r_0,t_0)G_{\hat{s}}^{\left(0\right)}(r^{\′},t^{\′}|r_0,t_0)\mathrm{\δ}{\widetilde{\mathrm{\μ}}}_s(r,\hat{s},{\hat{s}}^{\′})$

假设迁移方程的解是唯一的，方程(11)可以被减化到： $\frac{1}{4\mathrm{\π}}{G}_{\hat{s}}^{\left(0\right)}(r,t|r_0,t_0)=-{\∫}_{t_0}^{t^{+}}{\mathrm{dt}}^{\′}\∫d^3{r}^{\′}{\mathrm{\δ\μ}}_a\left(r^{\′}\right)G_{\hat{s}}^{\left(0\right)}(r,t|r^{\′},t^{\′})G_{\hat{s}}^{\left(0\right)}(r^{\′},t^{\′}|r_0,t_0)---\left(12\right)$ $-{\∫}_{t_0}^{t^{+}}{\mathrm{dt}}^{\′}\∫d^3{r}^{\′}{\mathrm{\δ\μ}}_s\left(r^{\′}\right)G_{\hat{s}}^{\left(0\right)}(r,t|r^{\′},t^{\′})G_{\hat{s}}^{\left(0\right)}(r^{\′},t^{\′}|r_0,t_0)$ $+{\∫}_{t_0}^{t^{+}}{\mathrm{dt}}^{\′}\∫d^3{r}^{\′}\underset{4}{\∫}\underset{\mathrm{\π}}{\∫}d{\mathrm{\Ω}}^{\′}{G}_{\hat{s}}^{\left(0\right)}(r,t|r^{\′},t^{\′})G_{\hat{s}}^{\left(0\right)}(r^{\′},t^{\′}|r_0,t_0)\mathrm{\δ}{\widetilde{\mathrm{\μ}}}_s(r,\hat{s},{\hat{s}}^{\′})$

方程(12)右侧的第一项是归因于吸收变化的修正，第二和第三项包含归因于散射变化的修正；第三项还包括归因于相位函数变化的附加修正。在式(12)中最后一项是面积分并且与边界条件有关。对于普遍采用的边界条件，例如在我们的系统中下面描述的零-边界条件，面积分在高阶作贡献，因此可以舍弃。为了使公式简化，把点展宽函数定义为： $\mathrm{PSF}(r,t;r^{\′}{;r}_0,t_0)=4{\mathrm{\π}\∫}_{-\∞}^{t^{+}}{\mathrm{dt}}^{\′}{G}_{\hat{s}}^{\left(0\right)}(r,t|r^{\′},t^{\′})G_{\hat{s}}^{\left(0\right)}(r^{\′},t^{\′}|r_0,t_0).---\left(13\right)$

所以，方程(7)和方程(12)可以被重新写成： $-[G_{\hat{s}}(r,t|r_0,t)-G_{\hat{s}}^{\left(0\right)}(r,t|r_0,t_0)]={\∫d}^3{r}^{\′}{\mathrm{\δ\μ}}_a\left(r^{\′}\right)\mathrm{PSF}(r,t;r^{\′};r_0,t_0).-----\left(14\right)$

请注意在方程(2)和方程(14)之间显著的相似性。实际上，PSF(r，t；r’；r₀，t₀)是光子在时间t₀在源点r₀被注入、在时间t之前途经场测量点r＇而且在时间t在点r被检测器收集到的概率。它代表在时间t在光源和检测器之间的光子路径分布。在弹道范围内，横向(即，垂直于r-r₀的方向)的光子路径分布收缩到δ函数： $\mathrm{PSF}(r,t;r^{\′};r_0,0)\stackrel{t\→|r-r_0|\mathrm{fc}}{\→}{\mathrm{\δ}}^{\⊥}(r-r_0)-----\left(15\right)$

因此，方程(14)对r，h，s的体积分被简化成沿着r-r_o的线积分，而且方程(2)和方程(14)本质上是同一的。

重要的是注意方程(14)的推演不使用扩散近似的假设。方程(13)中的格林函数G_ ⁽⁰⁾(r，t|r＇，t＇)可以使用各种不同的模型计算出来。三个例子是路径积分解，随机游动解和传统的扩散近似解。关于定时选通成像法的进一步的细节在美国专利第5,919,140号中能够找到，在此通过引证将其全部内容并入。事实上，这些重建表明传统的扩散解不适合用早期到达的光子成像。它预示太宽的点展宽函数并且提供太小的图像。因果关系令人满意的解是必需的。

为了比较，考虑扩散近似：

其中μ_tr是定义为μ_tr＝μ_a+μ＇_s的散射输运特性，而μ＇_s＝(1-g)μ_s，其中g是向前散射角的平均余弦。从方程(11)可以看出对格林函数的一阶修正是： $G^{\left(1\right)}(r,t|r_0,t_0)=-{\∫d}^3{r}^{\′}{\mathrm{\δ\μ}}_a\left(r^{\′}\right){\∫}_{t_0}^{t^{-+}}{\mathrm{dt}}^{\′}[G^{\left(0\right)}(r,t|r^{\′},t^{\′})G^{\left(0\right)}(r^{\′},t^{\′}|r_0,t_0)------\left(17\right)$ ${+\frac{1}{3{\mathrm{\μ}}_{\mathrm{tr}}^2}\▿G}^{\left(0\right)}(r,t|r^{\′},t^{\′})\·\▿\′G^{\left(0\right)}(r^{\′},t^{\′}|r_0,t_0)]$ $-{+\frac{1}{3{\mathrm{\μ}}_{\mathrm{tr}}^2}{\∫d}^3{r}^{\′}{\mathrm{\δ\μ}}_s^{\′}\left(r^{\′}\right){\∫}_{t_0}^{t^{+}}{\mathrm{dt}}^{\′}\▿G}^{\left(0\right)}(r,t|r^{\′},t^{\′})\·\▿\′G^{\left(0\right)}(r^{\′},t^{\′}|r_0,t_0).$

上述的方程与用扰动理论直接依据散射方程计算的修正是一致的。

本发明包括用于光学CT的级数展开方法和代数重建技术(ART)。待成像的体积被分割成N×N×N立方体的voxel。每个voxel被赋予一个代表吸收分布的本地平均值的数值。成像问题在X射线的情况下是二维的。然后，方程(2)中的线性衰减被简化成这些voxel数值沿着光源-检测器线的总和。光线总数可能是完全地表示成： $y_j=\underset{i=1}{\overset{N^2}{\mathrm{\Σ}}}{b}_{\mathrm{ij}}{w}_{\mathrm{ij}}{\mathrm{\μ}}_i------\left(18\right)$

在方程(18)中，y_j是第j次测量的数据；w是与第j次测量的倾斜角相关的几何学因子，而且它等于在voxel i之内第j条射线的线段长度；μ_i是在voxel i之内吸收的局部平均值。在方程(18)右侧的总和是在检测平面上N²个voxel的总和。在X射线的情况下，因子b_ij被给定为：

如图2A所示，通过引入权重因素b_ij，方程(18)中的平面求和实际上是线求和。那些落在线上的voxel对X射线函数做出100％的贡献，而那些落在线外的给出0％的贡献。

我们对光学CT的拓展是通过权重因子b_ij。因为散射，光子沿着不同的路径行进具有不同的概率。所以，每voxel对光源-检测器测量的贡献是用横越它的光子路径密度加权的。对于光学情况，b的数值是光子路径概率的数值(图2B)；然后，可以直接应用于方程(18)。回忆一下，

取消光子路径是三维的管子，在三维上求和是必要的。广义方程是： $y_j=\underset{i=1}{\overset{N^3}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}{w}_{\mathrm{ij}}{\⟨{\mathrm{\δ\μ}}_a\⟩}_i.-----\left(20\right)$

这完全地是在固定的时间t方程(14)的离散形式，如果我们设定：

其中Δv是每个voxel的体积元素。方程(20)可以用如下紧凑的矩阵形式重写：

                                        (22)

                  y＝Rx+n，

其中y是测量矢量(维数M)；x是图像矢量(维数N³)；R是包含几何和光子路径信息的射影矩阵(维数M×N³)；而n表示误差矢量。

在图像重建方面，把X射线CT的加法ART直接应用于光学CT。由于成像问题可能不是超定的就是欠定的，因此直接求解方程y＝Rx是不适当的，因为它可能全然没有任何解，或者它可能有许多解，但是没有一个是适当的。

成像矢量的估计通常是根据向前模型的预测结果与实验数据之间的误差以及关于成像区域的先验的信息使用优化判据完成的。一个例子是作为Bayesian估计的特殊情况的所谓的规则化的最小二乘准则。重建试图寻找的图像是通过使下面的函数减至最小完成的：

                                           (23)

       r²‖y-Rx‖²+‖x-δμ₀‖²

在方程(23)中，r是在文献中通常称之为信噪比的参数，δμ₀对于成像矢量(N³×1)是预先了解的并且也被用作图像的初始估计，而‖...‖²表示矢量的模的平方。保持第二项小将保证照片离预先了解的不太远，而保持第一项小则保证照片与测量结果一致。

一种使方程(20)最小化的叠代程序将分别引入N³维和M维的矢量x^(k)和u^(k)的两个序列。最初，设定x⁽⁰⁾等于δμ₀，而将u⁽⁰⁾设定为零矢量。通过引证在此并入的叠代步骤【G.T.Herman，Reconstruction From Projection：The Fundamentalsof Computerized Tomography(依据投影重建：计算机化的断层摄影术的基本原理)(Academic Press，New York，New York，1980)】是用下式给定的：其中： $c^{\left(k\right)}=\mathrm{\λ}\frac{r(y_{i_k}-\underset{j=1}{\overset{N^3}{\mathrm{\Σ}}}{R}_{i_k,j}{x}_j^{\left(k\right)})-u_{i_k}^{\left(k\right)}}{1+r^2{\left|\right|R_{i_k}\left|\right|}^2}-----\left(25\right)$

其中i_k＝k(mod N³)+1，R_i是矩阵R的第i行的转置矩阵，e是在第i行中为1、否则为0的M维的列矢量，而λ是称之为驰豫参量的实数。业已证明只要0＜λ＜2，序列{x^(k)}在极限情况下收敛到方程(23)的最小值(minimizer)。

为了证明这项光学CT技术，测量是在装满了像组织一样的散射媒体的边长6.35cm的立方体玻璃容器中进行的。选择立方体的几何形状是因为它的边界规则便于简化理论模型。散射媒体是用1.072微米直径的苯乙烯珠子(PolyScience公司)以0.27％的浓度制备的储备溶液。媒体的散射和吸收性质都是通过调整透射的漫射光实时地确定的。在图3A和3B中给出系统100的示意图。光源102是按毫微微-模式操作并且用Coherent Innova 400型多线氩离子激光器104抽运的CoherentMira 900型Ti/蓝宝石锁模激光器。波长是800nm，而重复频率是76MHz。脉宽是大约150毫微微秒。检测系统106是带M 5675型同步扫描单元的Hamamatsu条纹照相机C5680。激光束光的一小部分被石英板反射到快速光敏二极管108(Hamamatsu C 1808-02)，后者产生用于条纹照相机的触发信号110。立方体玻璃容器被安装在用来在光源、检测器和待扫描的物体之间激励相对运动的平移平台。奔腾-2计算机120作为中央控制和数据获得单元。为了在扫描期间监视激光功率漂移，DT2801-A卡(Data Translation，Inc.)被安装在计算机上并且为了记录来自光源的控制盒的激光器供电电压而编程，另外它还为驱动平移平台的步进电机服务。用于全行扫描的数据获得的整体自动化是通过给条纹照相机软件的用户界面编程实现的。条纹照相机是按模拟模式操作的。它把光信号的瞬时演变转换成垂直的条纹图像。透射的光是用四条粘在一起的纤维束130收集的。每条光纤束都具有500μm的芯径而且是由一万条单模石英纤维组成的。四条纤维的近端被捆在一起以增加收集面积。检测系统的总时间分辨率被设定为30ps。

容器被放在平移平台上的样品架之内。多个物体在混浊媒体内部被镶嵌在固定的位置。样品架是这样设计的，以致容器的顶面、底面和左右边界全都是黑色的。激光脉冲从前面传输到媒体之中，而透射的信号是对面收集的。入射的激光束和收集纤维按同轴的几何形状对齐。表面扫描是在XZ和YZ方向上实施的，所以吸收分布的3D图像可以被重建。扫描面积是沿着每个方向4cm×4cm见方。扫描在水平和垂直两个方向上都是每步2毫米，从而导致2个投影，总共882次扫描测量。用于每个点的数据获得时间是8秒。一两个不透明的物体被镶嵌在媒体之中。在每种情况下都首先对XZ平面进行扫描，然后容器旋转90°，继续对YZ平面扫描。

点展宽函数的计算需要了解散射媒体的平均散射和吸收性质(μ’_s和μ_a)。为了确定μ’_s和μ_a，随着吸收体从散射媒体中移开和在X-表面的中心对齐，2.2ns时间窗的同轴传输信号被收集起来。这条随时间变化的曲线是使用零边界条件的散射近似解拟合的。最好的拟合是通过μ’_s＝7.38cm^-1和μ_a＝0.03cm^-1给出的。这些数值类似于人类乳房组织的那些数值。

在这些测量中，数据是为两种配置收集的。在第一种情况下，8毫米直径的黑色球体被放在容器的中心(“一个物体的配置”)。在第二种情况下，8毫米直径的黑色球体放在距中心(0.56，56，-0.56)厘米的被置，而6毫米直径的第二个黑色球体被放在距中心(-0.56，-0.56，56)厘米的位置(“两个物体的配置”)。对于每种情况，扫描是按2毫米的步幅在表面上进行的，而透射信号的882条时间演变曲线是实测的，每个扫描点一条。因为光源-检测器距离对于所有的882次测量保持一样的，时间刻度是一样的，而强度可以在同一时间点进行比较。当光源-检测器线穿越吸收体的时候信号水平下降，反之亦然。投影强度图通过绘制所有的882条时间演变曲线在同一时间点的强度等值线图就能够轻易地产生。为了达到比较高的计数还降低噪音，这样的强度可以是在可与检测系统的时间分辨率比较的时间窗上的计数总和。时间点的选择是以空间分辨率和信噪比水平的折衷为基础的。早期时间部分有比较好的空间分辨率，但是低信号水平遭受相对比较高的噪音并且将扭曲图像。后期时间部分有比较高的信号水平，但是空间分辨率下降。用于求和的时间窗的选择如下：对于一个物体的配置，时间窗是534-600ps；对于两个物体的配置，时间窗是607-657ps。在这里计时零是飞行时间。投影强度图是在图5中提交的。如同预期的那样，穿越吸收体的投影线具有比较弱的强度。吸收体在投影强度图上是作为阴影出现的。嵌入一个吸收体的投影(图5A和5B)不同于嵌入两个吸收体的那些(图5C和5D)。图5a-5D被称为零阶图像。它们已经展现了诸如嵌入物体的中央位置之类的特征，尽管这些图像由于散射是杂乱的。

在处理图5a-5D的数据之前，由于极小的气泡和表面的不对称性在边缘造成的人为现象是人工去除的，并且格子被产生。在这项分析中，选择2毫米的扫描步幅作为每个格子的大小，从而导致21³＝9261个voxel。当9261个voxel数值是依据882个测量结果重建的时候，逆元是欠定的。

然后，成像重建程序是在两个步骤中应用的。首先，X射线CT的直线PSF被应用于在初始估计值δμ＝0(见方程(23))和“信噪比”r＝60的情况下使用ART(方程(24)和(25))的数据。其次，用直线PSF获得的合成图像随后被作为初始估计值使用并且被应用于使用特定的光子迁移PSF的数据。方程(23)中的“信噪比”再一次被设定为r＝60，而方程(25)中的驰豫参量被设定为λ＝1。重建是用从扩散近似解和因果关系修正解(见下文)两者计算出来的点展宽函数计算的。叠代步骤的次数被设定为2×10⁵。在Pentium-II 450MHz PC上，每次重建花费大约40秒。叠代步骤的次数增加到2×10⁶没有观察到更多的改善。

第二个步骤包括依据迁移方程的解计算PSF。扩散近似解在在文献中是被广泛使用的。然而，众所周知，这个解违反因果关系而且在早期时间状况中崩溃。具体表现因果关系的解对于以随机游动和路径积分理论为基础的模型已被解决。具体表现因果关系并且对于早期到达的光子有效的格林函数已被使用。这个格林函数是通过用t₀+t_ir(其中t_ir是未被散射的光子从光源传播到检测器的飞行时间)代替把脉冲注入媒体的时间t₀依据扩散近似格林函数G(r，t|r₀，t₀)构成的。具体地说，这个程序考虑到散射在光脉冲横越媒体之后才开始这一事实。

图6A-6D展示用扩散近似解和因果关系修正解计算的点展宽函数。请注意，图6A-6D中所有的等值线图都对应于供两个物体情况使用的实验数据的时间窗。使用扩散近似计算的点展宽函数比使用因果关系修正程序计算的宽。

对于X射线(直线PSF)、扩散近似、因果关系修正和实际配置的图像重建结果是在图7A-7D和图8A-8D的等值线图表中展现的。图9A-9D用二维等值线图给出图8C的特征切片。这些切片是在z＝-10毫米、z＝-6毫米、z＝2毫米和z＝10毫米的平面上挑选的。如同在图7A-7D和8A-8D中清楚地看到的那样，用直线路径(X射线CT)重建的图像没有高分辨率，而用扩散近似重建的图像与球体的实际尺寸相比太小，对于两个物体的情况尤为如此，对于这种情况较小的物体是看不见的。由于实验数据中的噪声，对于一个物体的配置重建的图像呈现一些扭曲。显然，对于两个物体的配置，用因果关系修正重建的图像给出镶嵌物体的正确的尺寸。如图9所示，voxel数值离得越远物体就越近乎于零，这自然是起因于反演程序。在图8C中，8毫米物体的大小和位置和6毫米物体的大小都是正确的，而6毫米物体的位置偏离其实际位置2毫米。这可能表明当镶嵌两个物体的时候外来干扰存在于逆元之中。

用传统的扩散近似解计算的点展宽函数就我们的数据的早期时间窗而言太宽。在这种情况下重建的图像与物体的实际尺寸相比太小。尤其是对于两个物体的配置，比较小的吸收体在重建的图像中基本上是看不见的。另一方面，用因果关系修正进行的同样的计算提供改善的分辨率。

图10图解说明依照本发明的优选实施方案的处理序列200。在把光分布函数和/或任何参考数据储存在电子存储器中202和为了处理针对特定类型的剖析或组织结构(例如，癌性损伤)收集的图像数据为计算机编程204之后，用内窥镜或探针扫描患者或活检样品206。为了产生用于显示和进一步处理210的图像处理收集到的数据208。

虽然这项发明已参照其优选实施方案被具体地展示和描述，但是熟悉这项技术的人应该理解不脱离权利要求书所囊括的本发明的范围可以在形式和细节方面做出各种各样的改变。

Claims (40)

1.一种使物体成像的方法，该方法包括：

提供光分布函数，该函数有散射分量和吸收分量；

把光引向待成像的物体；

检测物体发出的光；

用检测到的光形成物体的电子表达；以及

使用光分布函数处理该电子表达以便形成物体的图像。

2.根据权利要求1的方法进一步包括把波长在700nm到900nm范围内的光引向在物体。

3.根据权利要求1的方法，其中物体包括这样的组织，以致该方法进一步包括形成组织的图像。

4.根据权利要求1的方法进一步包括提供光源、检测器和接到检测器上的数据处理器。

5.根据权利要求1的方法进一步包括提供包括级数展开的光分布函数。

6.根据权利要求1的方法进一步包括提供用来检测光的收集时间，该收集时间少于1000ps。

7.根据权利要求4的方法，其中提供光源的步骤包括提供激光器。

8.根据权利要求4的方法，其中提供检测器的步骤包括提供条纹照相机。

9.根据权利要求1的方法，其中光分布函数包括点展宽函数。

10.根据权利要求9的方法进一步包括提供众多权重函数。

11.根据权利要求1的方法进一步包括确定组织中癌性损伤的大小。

12.一种用来使物体成像的系统，该系统包括：

有光发布函数的数据处理器，该函数有散射分量和吸收分量；

把光线传输到待成像的物体上的光传输系统；

检测物体发出的光线的光传感器，该传感器被这样连接到数据处理器上，以致物体的电子表达是用检测到的光线形成的，该电子表达是用光分布函数进行处理的，以便形成物体的图像。

13.根据权利要求12的系统进一步包括把波长在700nm到900nm范围内的光线引导到物体上。

14.根据权利要求12的系统，其中物体包括组织，而且进一步包括接在数据处理器上显示组织图像的显示器。

15.根据权利要求12的系统进一步包括与检测器排成直线的光源。

16.根据权利要求12的系统进一步包括其中包括级数展开的光分布函数。

17.根据权利要求12的系统进一步包括控制用来检测光线的收集时间的控制器，该收集时间少于1000ps。

18.根据权利要求15的系统，其中光源包括激光器。

19.根据权利要求12的系统，其中传感器包括条纹照相机。

20.根据权利要求12的系统进一步包括在待成像的物体和传感器之间提供相对运动的扫描器。

21.根据权利要求20的系统进一步包括控制扫描器、门检测器、光源和数据处理的控制器。

22.根据权利要求12的系统进一步包括众多光分布函数。

23.根据权利要求12的系统进一步包括纤维光学光耦合器。

24.根据权利要求12的系统进一步包括用于插入身体把光线输送给组织的探头。

25.根据权利要求12的系统进一步包括众多权重函数。

26.一种使患者成像的方法，该方法包括：

提供光分布函数，该函数有散射分量和吸收分量；

提供患者体内组织的电子表达；以及

使用光分布函数处理该电子表达以便形成物体的图像。

27.根据权利要求26的方法，其中从患者处收集到的光线具有在700nm到对900nm范围内的波长。

28.根据权利要求26的方法进一步包括形成组织内癌性损伤的图像。

29.根据权利要求26的方法进一步包括提供用光分布函数编程的数据处理器。

30.根据权利要求26的方法，其中光分布函数包括级数展开。

31.根据权利要求26的方法进一步包括提供收集时间，在此期间从患者处收集光线，该收集时间少于1000ps。

32.根据权利要求26的方法进一步包括诸如条纹照相机之类的检测器。

33.根据权利要求26的方法进一步包括提供有级数展开分量的光分布函数。

34.根据权利要求26的方法，其中光分布函数包括点展宽函数。

35.根据权利要求34的方法进一步包括提供众多权重函数。

36.根据权利要求26的方法进一步包括确定组织内癌性损伤的大小。

37.根据权利要求26的方法进一步包括用光纤装置收集光。

38.根据权利要求26的方法进一步包括在待成像的身体之内定义具有众多voxel的图像体积，每个voxel都有权重因子。

39.根据权利要求26的方法，其中光分布函数包括迁移方程近似。

40.根据权利要求26的方法，其中光分布函数定义众多有横截面积的光路径有跨部份的区域，该面积小于其扩散近似。

Images