CN1678900A - 用于内部光源的3d成像的方法和设备 - Google Patents

Abstract

本发明提供了一种用于获得在诸如哺乳动物之类的样本内的一个或多个光源的三维(3D)表示的系统和方法。哺乳动物组织是混浊介质，意味着当光子通过组织传播时都被吸收和散射了。如果与吸收相比较散射更大，例如红色到近红外光穿过组织，那么就用扩散理论来描述光在样本内的传输。使用成像数据和由计算机实现的光子扩散模型，本发明的实施例生成在样本内光源的3D表示，例如这种光源的3D位置、大小和亮度。

Description

用于内部光源的3D成像的方法和设备

技术领域

本发明总体上涉及利用光进行的成像。特别是，本发明涉及用于获得在混浊样本(turbid sample)内光源或光源分布的三维表示的系统和方法，这在生物医学成像和研究应用中尤为有用。

背景技术

生物发光成像是一种非侵入技术，用于在医学研究、病理学和药物发现和开发领域内对动物研究对象执行活体内诊断性研究。生物发光典型地由已经被诸如荧光素酶之类的发光指示器感染的细胞生成，并且可以用作区分特定组织类型(例如瘤)、监视生理机能、跟踪研究对象服用的治疗化合物的分布或疾病的进展。已经证实了大量的应用，包括肿瘤学、传染病和转基因动物领域。用发光指示器来标记细胞的活体内成像是相关的技术，其最近在绿色荧光蛋白(GFP)和诸如Cy5.5之类的近红外(NIR)染料中也已经得到证实。

由生物发光细胞所发射的光子在研究对象的组织中被强烈地散射，以致实际上传播是扩散的。虽然通过组织扩散的许多光子被吸收了，但是小部分会到达研究对象表面并且可以被检测到。一般说来，在哺乳动物组织中的吸收在光谱的蓝绿色部分(＜600nm)高，而在光谱的红色和NIR部分(600-900nm)低。萤光虫的荧光素酶具有相当宽的发射光谱范围从500-700nm，因此至少部分发射处于低吸收区域。因为用于在组织中散射的平均自由路径不长，大约～0.5mm，所以来自深源的光子在到达表面之前常常被多次散射。生物发光成像系统实际上记录这些发自研究对象表面的光子的空间分布。

然而，最重要的定量信息并不直接与表面发射相关，而是代之以与在研究对象内的生物发光源有关。重要的参数是源强度(与光发射细胞的数目相关)、位置和几何结构。迄今为止公布的大部分生物发光成像工作涉及单视图2D成像系统的使用。图像分析通常涉及定量在研究对象表面上的光发射感兴趣区域(ROI)。虽然此分析方法简单并且提供了好的光发射的相对测量，但是它没有考虑源深度和所产生的穿过组织的衰减。

因此，有兴趣开发改进的成像系统和重构算法，其会根据在样本(例如动物)表面上测量的图像来提供在所述样本内光子发射的三维分布。

发明内容

本发明提供了用于获得在诸如哺乳动物之类的样本内一个或多个光源的三维(3D)表示的系统和方法。哺乳动物组织是混浊介质，这意味着当光子传播通过组织时即被吸收也被散射。在与吸收相比散射更大的情况下-例如红色到近红外光穿过组织，那么就用扩散理论来描述光在样本内的传输。本发明的实施例，使用成像数据和由计算机实现的光子扩散模型，生成在样本内光源的3D表示，例如这种光源的3D位置、大小和亮度。

在一个方面，本发明涉及一种用于获得在混浊样本内光源分布的三维表示的方法。根据发自所述样本表面的光来提供表面光图像数据，所述光来源于在样本内光源的分布。然后，根据这样提供的表面光图像数据来获得在样本表面内的光数据。在一个具体实施例中，把表面光图像数据转换为刚好在所述样本表面下的光子密度，然后，此光子密度用于获得在所述样本表面内的光数据，由此产生在样本内光源分布的三维表示。

在又一方面，本发明涉及一种用于获得在样本内光源分布的三维表示的成像系统。所述成像系统包括具有一组壁的成像室，所述一组壁围住内腔。所述成像室还包括照相机架和可移动载物台(stage)设备，其中所述照相机架被配置成用于安置照相机。所述可移动载物台设备包括传送机构和载物台，所述载物台被配置为支撑在所述内腔之内的样本。所述载物台与所述传送机构连接，用于把所述样本移动到在所述内腔中多个位置之一。所述成像室还包括用于传送发自样本表面的光的光传输机构。所述成像室还包括处理器，该处理器被设计或配置成根据发自样本表面的光来提供表面光图像数据，所述光来源于在所述样本内光源的分布。还把所述处理器配置为获得在所述样本表面内的光数据，以便根据所述表面光图像数据来获得所述光源分布的三维表示。

在又一方面，本发明涉及一种计算机程序产品。所述计算机程序产品包括计算机可读介质和经由所述计算机可读介质提供的程序指令。所述程序指令包括用于获得位于样本内光源的三维表示的重构指令。所述重构指令能够执行一个或多个上述方法操作。

结合下列附图，下面在本发明的具体实施方式中将要更详细地描述本发明的这些及其它特征。

附图说明

在附图的图中作为举例而不是作为限制说明了本发明，并且在附图中相同的参考数字指的是同样的元，其中：

图1是依照本发明一个实施例的适于生成位于样本内的光源的3D表示的成像系统的透视图。

图2A是依照本发明一个实施例的具有用于使样本的多视图便于进行的内部组件的图1的室的剖视图。

图2B是依照本发明一个实施例的用于使样本的多视图便于进行的图2A的内部组件的透视图。

图2C示出了所述载物台相对于图2A和2B的光传输装置的四个不同的示例位置：0度、90度、180度和270度。

图2D示出了用于获得研究对象表面形貌的结构光投射器的配置。

图3A示出了依照本发明一个实施例在混浊介质中光成像的简化说明。

图3B说明了示出在表面边界处平面近似的示意图。

图3C说明了圆柱体的体积元(体元)和面积元之间的关系。

图4A说明了用于在2D照相机数据和示例性面积元的表面数据之间转换的示例性关系。

图4B-4H说明了对应于图5C的处理流程的结构光成像的图示。图5A说明了依照本发明一个实施例的用于获得在样本内一个或多个光源的3D表示的处理流程。

图5B说明了依照本发明一个实施例的用于使用图1A的成像系统以获得成像数据的处理流程。

图5C说明了依照本发明一个实施例的用于使用图1A的成像系统以获得表面形貌数据的处理流程。

图5D说明了依照本发明一个特定实施例的用于把表面光发射数据转换为在所述表面内的光的处理流程。

图5E和5F表示了用于在仿真老鼠(例如具有嵌入式光纤的塑料老鼠)内点状光源的重构结果。

图6说明了依照本发明一个实施例的示例性的计算机系统。

具体实施方式

在本发明下面的详细说明中，阐述了大量具体实施例以便提供对本发明彻底的理解。然而，正如对那些本领域内技术人员来说显而易见的一样，在没有这些具体细节的情况下或通过使用可替换的元素或过程可以实施本发明。在其它例子中，没有详细描述众所周知的过程、组件和设计，以免不必要地使本发明的特征变得不明显。

1. 成像系统

本发明的3D光重构技术可以在任何适当的系统内加以实现。其中适于实现本发明技术的成像系统的几个实施例还在美国专利申请号09/905,668中进行了描述，所述专利申请由Nilson等人在2001年7月13日提交，题目为Multi-view imaging apparatus(多视图成像设备)。

图1说明了依照本发明一个实施例的适于生成在样本内光源的3D表示的成像系统10。所述成像系统10可以用于各种任务。例如，这些任务可以包括捕获照相的、发光的以及结构光图像。这些任务还可以包括执行各种处理和成像任务，例如获得、分析和处理3D光源表示。

所述系统10一般使用感光照相机对低亮度光源进行成像并加以处理。来自光源的光指的是波长范围在400-950nm的可见到近红外(NIR)光谱部分中任何地方的光子和电磁能量。当然，在系统10中成像并处理的一些亮度不容易被人类视觉所察觉。例如，发自样本的低亮度光的表面辐射度可以在大约10³到大约10¹⁰光子/秒/平方厘米/球面度之间，其中此范围的下端一般完全在人类检测等级之下。

在所说明的实施例中，所述成像系统10包括适于接收发光样本的成像室，其中将检测诸如基于荧光素酶的荧光之类的低亮度光。高灵敏度照相机20与成像室12耦合，所述照相机20例如是增强或电荷耦合器件(CCD)照相机。所述照相机20能够在成像室12内捕获样本的荧光的、照相的(即基于反射的图像)和结构的光图像。所述照相机20可以可选地由诸如致冷装置22之类的适当的源来冷却，所述致冷装置22经由导管24使低温流体穿过所述照相机循环。

图像处理单元26可选地分别通过电缆30和32在照相机20和计算机28之间连接。可以是任何适当类型的计算机28典型地包括包含硬件的主机36，所述硬件包括处理器、诸如随机存取存储器(RAM)和只读存储器(ROM)之类的存储器组件和磁盘驱动器组件(例如硬盘驱动器、CD、软盘驱动器等)。所述计算机28还包括显示器38和诸如键盘40和鼠标42之类的输入装置。所述计算机2 8经由电缆34与在成像室12中的各种组件通信。作为选择，计算机可以与成像设备集成在单个单元中。

为了向系统10中的组件提供控制和通信，计算机28可以用任何用来控制并监视在系统10中的任何装置的适当的硬件和软件的组合来实现。由计算机28控制的组件可以包括照相机20、负责照相机20聚焦的电动机、照相机镜头、f制光圈、与包括在室12内的可移动载物台耦合的一个或多个电动机等。计算机28还与显示器38通信，用来向用户呈现成像信息并且还担当控制所述成像系统10的接口。

计算机28还包括用于如这里所述的图像处理和数据操作的适当的处理硬件和软件。例如，可以用软件来配置计算机28以使用来自样本表面的光发射测量来构建光源的3D表示。另外，可以把所述计算机28配置成使用2D结构光图像来生成3D表面形貌，所述2D结构光图像取自在内腔中载物台的一个或多个位置。下面进一步描述这些3D光构造技术。作为选择，所述计算机28可以与所述成像系统10的其它组件结合以形成独立的装置。

在一个实施方式中，下面进一步讨论的3D光重构技术需要在整个样本表面上的生物发光表面辐射度和表面形貌测量。因此，所述成像系统还包括用于捕获来自各种视图和位置的样本表面的图像。优选的是，所述成像系统满足下列目标。首先，希望所述照相机保持固定，因为用于冷却照相机的各种电缆和电源线太过笨重以致于在成像期间不能移动。其次，所述样本优选的是在成像期间保持水平。最后，优选地把所述成像系统配置为获得样本的3D表面形貌来作为下面进一步描述的重构技术的输入。可以使用任何适当的硬件和软件的组合来实现这些优选的目标。

图2A和2B说明了依照本发明一个实施例的用于使多视图成像便于进行的室12的各种内部组件。图2A是室12和其内部组件的剖视图，而图2B是依照本发明一个实施例的图2A的内部组件的透视图。在所说明的实施例中，样本位于连同镜子一起旋转的水平透明平台上，所述镜子绕水平光轴旋转以便从各个视角捕获来自所述样本的光，并且把所捕获的光指向照相机以便成像。

如同所示，把可移动载物台设备以载物台204和传送机构202的形式布置在室12(图2A)的内腔中，所述载物台204用于支撑发光样本106，所述传送机构202用于移动所述载物台。光传输装置120被包括，以在其各种位置上收集来自所述样本的光。在图2A和2B所说明的实施例中，以光束252a的形式发自样本的光一般从光传输装置120反射出来为光束252b。所述光传输装置把此反射光252b导过孔径122和透镜262(图2B)，以便由照相机20(图1)的检测器264(图2B)捕获图像。

在一个实施方式中，光传输装置120包括把来自在载物台204上样本106的光反射穿过孔径122的斜镜121。外壁123基本上是圆柱形的并且包括孔径122，所述孔径122使光能够从载物台204上的样本106经由镜子121传到成像透镜262(图2B)。光传输装置120的外壁123还防止在室12的内腔中的余光到达透镜262，所述余光不直接与载物台204的当前视角相关联。这通过把镜子121配置得足够长以便跨过载物台204的长度来部分地执行。由于把所述载物台定位在关于所述光传输装置120的固定轴的各种位置上，所以外壁123和镜子121协作来主要地从载物台204的角度方向收集光，然后，把所述光反射向透镜262。

在所说明的实施例中，光传输装置120可旋转地与室12耦合，并且在方向258上绕着固定轴旋转。所述可移动载物台设备具有两个自由度，以在相对于所述光传输装置120的多个位置定位所述载物台204(和样本106)以便捕获来自所述样本不同侧面的光，从而形成所述样本的多视图图像。在没有相对于载物台204重新定位样本姿势的情况下，所述可移动载物台设备允许照相机从不同的视图、角度和位置在所述成像室内捕获所述样本或其部分的图像。图2C示出了所述载物台相对于光传输装置(120)的四个不同的示例位置：0度、90度、180度和270度。

所述可移动载物台设备可以由任何适当的数目、类型和布置的组件来组成，以便实现所述样本的多个定位和观察。在上面引用的Nilson等人的专利申请中描述了多视图成像系统的几个实施例。优选的是，所述可移动载物台设备还包括用于检测所述样本位置的传感器机构，以用于下面进一步描述的后续数据分析过程。

在一个实施例中，在镜子部件120的镜筒(未示出)内提供了光源以照亮在所述成像室12中的样本或标本。所述光源可以连续地照明或闪光以收集样本的照相图像，并且当捕获荧光图像时关闭。

优选的是把结构光发光器集成到用于重构所述样本的表面形貌的系统中。在通常的实施方式中，所述成像系统优选地包括用于从某一角度把方格线投射到所述样本上的机构，所述角度例如正常来说是从20到30度。当方格线遇到有限高度的物体时，则相对于所述平台移动或相移所述方格线。使用用于结构光分析的常规公式，可以根据所述线的相移来确定在所述样本全部表面上的表面形貌数据。所述线距可以基于样本表面纹理和样本大小来变化，但是在每毫米0.5到2线范围内的线距可以适于老鼠样本。较紧密的线距提供了高分辨率，但是在诸如毛皮之类的粗糙表面上所述线更难以被追踪。依照下面所描述的技术使用表面形貌数据来构建3D表示，例如为3D重构设置边界条件。

在具体实施例中，把简单的柯勒投影方案用作结构光源。在这种情况下，可以由扩散LED源来照亮所述划线，并且然后，以大约10x的放大倍数把所述划线投射到动物载物台上。在图2D中示出了该系统合并到系统10中的例子。所述投射器模块170跨在所述旋转镜部件120的后面，以致在所有视角上始终把线投射到所述样本106上。水平地投射所述光照图案并且在较大旋转镜的基部通过小镜子173反射以照亮样本106。

相对于图2A-2D所描述的这种成像系统可以从CA.Alameda的Xenogen公司中获得。尽管图2A-2D说明了用于获得来自样本106相对于照相机20的多个位置的数据的一个系统，但是本发明并不局限于这种系统。

2.  3D源重构

在诸如哺乳动物研究对象之类的许多混浊样本中，由在所述样本内的光源传播并且发自所述样本表面的光在所述样本内部(例如哺乳动物组织)被散射，以致实际上通过所述样本的光传播是扩散的(图3A)。当所述光通过样本扩散时，一些光被吸收了，而一小部分光到达所述样本表面。

通常以单元辐射度规定来自样本表面的光线发射，定义为光子/秒/平方厘米/球面度。校准这里所描述的成像系统以便以单元辐射度来报告表面亮度。使用用于在组织-空气接口处的光子传播模型，所述表面辐射度可以与刚好在样本表面内的光子密度有关。然后，使用扩散模型刚好在所述表面内的光子密度能与在所述样本内发光指示器的分布有关。因此，本发明使由成像系统测量的混浊样本的表面辐射度与在所述样本内发光的分布相关。更具体地说，本发明使用重构技术来生成内部光源的3D表示，所述重构技术利用了发自样本表面的光数据。所述重构技术使用由下列组成的输入数据集：a)样本的表面形貌表示和b)在至少一部分表面上的一组光辐射度的测量(例如表面图像)。为了使处理便于进行，可以把所述表面表示分成表面元，并把所述样本内部分成构成体积元网格的体积元或体元。由所述体积元网格的元来描述在样本内光源的分布。

图3A示出了依照本发明一个实施例的光成像的简化说明。图3A显示了由成像系统10成像的哺乳动物样本302，所述成像系统10使用用于光数据捕获的CCD照相机20。哺乳动物组织和许多其它混浊介质的特性使从在所述样本内的光源305到所述样本表面309传播的光出现光子散射307，导致在表面309上扩散成像。使用光子扩散模型和成像数据，本发明根据一个或多个表面309图像生成光源305的3D位置、大小和亮度。为了方便用户，例如可以经由在图1中的显示器38把3D表示表达为图示。

本发明的一个有益应用是使用表面光发射的一个或多个图像来重构在诸如老鼠之类的活体动物研究对象内的一个或多个生物发光或荧光源、或源分布。这在例如肿瘤学、传染病研究、基因表达研究和毒理学方面得到应用。本发明适用于具有复杂表面的样本，例如动物研究对象(例如老鼠)。这里使用的术语，复杂表面是指不能只用单个多边形描述的任何表面。这里描述的重构技术不对源分布有任何限制，例如在所述样本中的光源数目或所述源的大小和形状，并且不对所述表面的几何结构、大小或形状有任何限制。

在样本表面内的光数据通常指的是在所述样本内部的实际光的数学表示或近似。这可以包括一组点或体积元，每个以3D位置和源强度为特征。在一个实施例中，本发明把样本内部分成体积元，其中每个体积元被认为在其中心包含点光源。然后，这些体积元的立体网格定义了许多点源。这些点源近似于在所述样本内实际光源的分布。例如立方体积元的立体网格是合适的。

图3C说明了被分为体积元的立体网格的说明性的圆柱体，示出了其中的一个立体网格。每个体积元包含点光源314S_i，其在表面上产生光子密度316ρ_j。对于立方体积元表示，所述体积元大小可以依照所述样本大小和所要求解的精确性来变化。一般说来，增加立体网格密度会改善内部光数据表示的精确性。使所述立体网格的密度随样本的体积变化是可能的。特别是，希望增加在光源附近的所述立体网格的密度以在这个区域内提供更精确的表示，而在那些没有发生感兴趣动作的区域内(在对应于没有或最小光产生的区域中)降低所述立体网格的密度。

所述表面表示指的是样本的表面形貌，并且由实际表面的数学描述或近似来将其具体化。所述表面表示不必包括整个样本，而是可以包括与特定成像方案有关的样本的一部分。例如对于老鼠，所述表面表示可以不必包括整个老鼠，例如尾部的末梢部分和每只脚的末梢部分。因此，所述表面表示意指广泛地涉及所述样本的任何表面部分而不必是整个样本。典型地，所述表面表示包括在所述样本上一个或多个表面元或感兴趣区域，其生成与内部光源相关的表面光发射数据。为了方便用户，常常以诸如来源于结构光构造技术的3D描绘之类的图示描绘来显示所述表面表示。

包含所述光源的样本包括某种混浊的内部。在这种情况下的混浊内部是指不允许光无阻碍传输的体积。所述混浊内部可以包括一种或多种介质、结构、固体、液体、气体等。在一个实施例中，为样本建立模型为均匀的，以致在所述样本中的每个典型体积元具有相同的光传输属性的特性。在另一实施例中，把样本表示为不均匀的，以致在所述样品中各种典型的体积元具有不同的光传输属性的特性。在诸如老鼠之类的哺乳动物样本中，例如其内部可以包括组织、骨头、器官等的混合物，其中每一个都可以具有依照非均匀模型的独立的光传输属性的特性。

本发明给出了几种用于根据所述表面光发射来重构在样本内三维光源分布的方法，正如下面将要具体描述的一样。简短地说，首先获得一组结构光和荧光图像，例如利用上述的系统来获得。使用结构光算法重构所述动物的表面形貌并且产生表面网格。令整数j来计数所述表面元。使用一个或多个荧光图像、所述成像系统的方向和表面网格的方向，可以确定刚好在第j个元的表面下的光子密度ρ_j。刚好在所述表面下的光子密度与发自表面的光强有关(公式1)，正如下面将要描述的一样。可以在所述表面内的体积中构造一组体积元。令整数i来计数体积元。在每个体积元S_i中的源强度与刚好在所述表面下的光子密度通过格林函数核G_ij线性相关，其产生了一组线性方程(公式4)，正如下面将要描述的一样。

图5A说明了依照本发明一个实施例的用于获得在样本内一个或多个光源的3D表示的处理流程500。依照本发明的过程可以包括多达几个的附加步骤，这里不描述或说明它们以便不致使本发明变得不明显。

处理流程500使用处理发自样本表面的光数据的3D重构技术。所述重构技术使用其提供关于所述样本几何结构和来自样本表面光发射的空间分布的信息的数据。这可以包括下列数据，所述数据包括a)所述样本的表面表示和b)经由至少一部分表面的光发射的一组测量(例如图像)。用于此数据的一个格式包括(1)定义所述样本表面的表面表示，(2)来自样本表面的一组光发射图像和(3)定义在所述样本和每个图像之间的空间关系的一组参数。

处理流程500从获得用于3D重构技术的输入数据(502和504)开始。例如，处理流程500可以从获得表面表示(502)开始。为了使处理便于进行，可以把表面表示划分为近似于所述样本实际表面的表面元或适当的表面网格(502)。例如，可以由一组相连的多边形定义表面网格，其中每个多边形是表面元。表面元的数目可以依照样本大小和所要求解的精确性来变化。表面元的密度也可以在表面网格上从一点到另一点变化。一般说来，通过增加表面元的数目来改进表面表示的精确性。

处理流程还可以获得来自所述样本的表面光发射数据(504)。对于成像系统10，所述表面光发射数据包含在由照相机20拍摄的图像中。所述图像包括经由一个或多个表面部分的一组光发射的测量。在一个实施例中，根据样本相对于照相机的多视图来获得所述样本的多个图像。在这种情况下，每个图像提供了来自所述样本的表面光发射的不同的二维视图。可以拍摄多个图像以提供来自多个角度的附加数据。

然后，可以把所述图像数据映射回到所述样本表面上。因为照相机生成2D数据，所以依照在样本表面和照相机镜头之间的几何结构来处理所述图像数据以得到表面上的光发射强度(或辐射度)值。

图4A是用于把2D照相机数据转换为用于示例性表面元402的表面数据的说明性的关系。图4A示出了在表面元402、图像404和所述成像系统的入瞳406之间的关系。发自表面元402的光穿过入瞳406并且记录在图像404中。发射相对于表面法线的角度是θ₂。所述入瞳406对着小立体角dΩ。成像系统10可以以角度θ₂(相对于表面元402的法线来测量)把发自在样本上表面元402的光收集到由入瞳对着的立体角dΩ中。然后，可以使用此信息来把由照相机获得的图像数据转换成对应于表面几何形状的表面发射亮度。

然后，处理流程500把表面光发射数据转换为在所述表面内的光数据以获得三维表示(506)。所述表面发射亮度与刚好在样本表面内的光子密度有关。因此本发明把每个表面元的光发射亮度值转换为刚好在所述表面内的光子密度。再次参照图4A，在表面元处的发射亮度值I(θ₂)与在所述表面元下面的光子密度ρ有关。该关系的确切形式取决于用于描述穿过所述表面边界的光子传输的模型。根据分流边界条件[Haskell等人]，给出了此关系的一个实施例：

$I\left({\mathrm{\θ}}_2\right)=\frac{c}{4\mathrm{\π}{n}^2}T\left(\mathrm{\θ}\right)\cos {\mathrm{\θ}}_2\mathrm{d\Ω}[1+\frac{3}{2}\frac{1-R_{\mathrm{eff}}}{1+R_{\mathrm{eff}}}\cos \mathrm{\θ}]\mathrm{\ρ}---\left(1\right)$

这里c是光速，n是样本介质的折射率，T是光通过表面元离开样本的传输系数，以及θ是内部发射角，其通过Snell定律与外部发射角θ₂有关：

nsinθ＝sinθ₂                      (2)

所述参数R_eff是由下列公式计算的平均内部反射系数：

$R_{\mathrm{eff}}=\frac{R_{\mathrm{\φ}}+R_j}{2-R_{\mathrm{\φ}}+R_j}$

$R_{\mathrm{\φ}}=\underset{0}{\overset{\frac{\mathrm{\π}}{2}}{\∫}}2\sin \mathrm{\θ}\cos \mathrm{\θR}\left(\mathrm{\θ}\right)\mathrm{d\θ}$

$R_j=\underset{0}{\overset{\frac{\mathrm{\π}}{2}}{\∫}}3\sin \mathrm{\θ}{\cos }^2\mathrm{\θR}\left(\mathrm{\θ}\right)\mathrm{d\θ}$

因此，所述内部反射率R_eff取决于在表面元下面的介质的折射率。例如在组织中，R_eff典型地在0.3-0.5范围内。

因此可以使用公式(1)和(2)来把在每个表面元处测量的表面发射数据转换为在所述表面下面的光子密度值。作为这里所使用的项，下标j计数一组表面元。然后，ρ_j是在第j个表面元处计算的光子密度值。在图5D中提供了把表面光发射数据转换为在所述表面内的光数据以获得三维表示的进一步描述。

图5B说明了依照本发明一个实施例的用于使用图1的成像系统10来获得成像数据的处理流程520(从处理流程500的504)。处理流程520从把样本放置在成像室12内的载物台204上开始(521)，所述样本例如是用于光发射成像的标本或试样。使用计算机28，用户输入用于载物台204的所要求的位置。作为选择，基于自动化数据采集方法，用于载物台204的所要求位置是已知的。传送机构202依照由计算机28提供的控制信号把载物台204移动到所要求的位置(522)。光传输装置120还可以依照由计算机28提供的控制信号重新定位。然后，为捕获样本的照相图像来准备成像室12和相关的图像组件(523)。准备可以包括在计算机28上装入成像和采集软件(例如由CA Alameda的Xenogen公司提供的“LivingImage”)以及初始化照相机20。进一步的准备可以包括关闭门18，激活软件中的照相捕获选项，把照相机20聚焦到所述样本或动物的特定深度，并打开室12中的光。准备还可以包括调照相机20的镜头的焦距，有选择地放置合适的照相机20的镜头滤光镜，设置照相机20的f制光圈等。

然后，捕获照相图像(524)。在可选择的实施例中，在所述样本的照相成像期间使用“实况模式”以实时地观察所述样本。所述实况模式包括频繁地拍摄足够的照相图像序列以模拟实况视频。当完成照相捕获时，把照相图像数据转送到图像处理单元26和/或在计算机系统28中的处理器。可以使用这些来处理并存储所述照相图像数据以及处理所述用于显示在计算机监视器38上的数据。

在一个实施例中，本发明在图像捕获期间使用结构光。因此，随着载物台204在所要求的位置上，可以捕获结构光图像(526)。结构光图像捕获可以使用如上所述的结构光投影系统来完成，并且其还可以包括准备所述结构光投影系统(525)和用于所述结构光图像捕获中的成像系统10的任何其它组件。当完成时，把所述结构光图像数据转送到图像处理单元26和/或在计算机系统28中的处理器，其可以用于构建所述样本的3D表示。

随后，随着载物台204在所要求的位置上或在如先前进行照相或结构光图像的相同位置上，成像设备10为荧光图像捕获作准备(527)。这种准备可以包括使用所述计算机28来选择荧光的曝光时间和等级，并且关掉在室12内的光。当准备好时，所述CCD照相机20然后经由一个设定时段(多达几分钟)来捕获(514)所述荧光图像。把所述荧光图像数据转送到图像处理单元和/或在计算机28中的处理器。

此时，用户可以处理并存储所述荧光图像数据以及处理它以便将其显示在计算机显示器38上。所述操作还可以包括用照相图像来覆盖所述荧光图像以及把这两个图像一起显示作为2D“覆盖”图像，典型地以伪彩色示出所述荧光数据以便示出亮度。如所提及，所述光子发射数据可以表示在照相机20上的特定像素，所述照相机20在图像捕获期间的持续时间上检测光子。然后，此覆盖图像可以用于用户分析的基础；并且可以按要求分析并处理。特别的是，分析可以包括在一部分荧光表示内的像素上的照明量的求和。注意，尽管该论述将集中在用于所述覆盖图像的单个荧光表示，但是处理流程520可以包括从载物台204的相同位置取出多个荧光表示，例如同时或稍后(530)。

然后，可以把载物台204移到第二位置(531)。当所述载物台位于第二位置时，如上所述可以捕获所述样本的一个或多个照相的、结构光和/或荧光图像。可以通过捕获来自所述样本的可替换位置和视图的所述样本的图像来进一步继续图像采集。对于3D重构，可以从许多位置捕获所述样本的照相的、结构光和/或荧光图像。例如，图像捕获可以发生在成像室内样本的从2到200位置的任何地方。图像的优选数目是8，每45度隔开。在一个实施例中，经由计算机28使此过程自动化并控制该过程。

在已经从各种角度捕获适当数目的结构光图像并存储之后，计算机28然后可以处理来自每个角度的结构光数据以产生结构光表示(532)。每个结构光图像提供了所述样本大约半个面的表面形貌。

在另一实施例中，处理流程520和成像设备10使用图像序列来重构所述样本的3D表面形貌。通过从几个视角来获取图像，例如大约每45度，通过把从每个视图获得的部分表面重构”缝合”在一起可以重构所述样本的整个3D表面。然后，可以从不同视角获取图像序列并将其用于重构所述样本的3D表面形貌。还可以把3D表面形貌和图像数据用于重构在所述样本内光源的3D位置、亮度和大小。一旦由处理器28接收所述图像，则对所述数据应用适当的重构算法以获得3D表面形貌。正如本领域的技术人员理解的一样，存在许多用于根据结构光图像来重构表面的算法。例如，可以根据2D傅立叶变换来确定在所述图像上所有点处每条线的相移。在M.Takeda，H.Ina和S.Kobayshi，JOSA72，156-160(1982)的文章题目为“Fourier-transform method offringe-pattern analysis for computer-based topography andinterferometry(用于基于计算机的地形学和干涉量度学的边缘模式分析的傅立叶变换方法)”中详细描述了这种过程。然后，通过“展开”所述相位图来计算实际表面高度。在由D.C.Ghiglia和M.D.Pritt(John Whiley和Sons，纽约，纽约，1998)所写的教科书题目为“Two-Dimensional Phase Unwrapping，Theory，Algorithms，andSoftware(二维相位展开、理论、算法和软件)”中详细描述了此过程。可以一同将所述样本的结构光照相表示和所述样本的荧光表示相组合以形成结构光叠加或3D覆盖图像，典型地以伪彩色示出荧光数据以在视觉上表示亮度的特征。

图5C(从处理流程500的502开始)说明了依照本发明另一实施例的用于使用图1的成像系统10以获得表面形貌数据的处理流程510。图4B-4H说明了对应于处理流程510的结构光成像的图示。

处理流程510从对参考结构光进行成像以生成没有所述样本的图案开始(511和图4B)。这可以通过在成像样本之前把结构光施加到所述样本停放的载物台或表面上来执行。在所述样本的图像捕获期间，把所述载物台移到一般的位置，正如那些在没有所述样本的情况下用于图像捕获的位置。

随后当所述样本处于成像室时，用结构光来使所述样本成像(512和图4C)。结构光使用一连串的光线，所述光线以和表面法线成一定角度(例如大约30度)向下投射到样本上。当所述线越过所述样本时发生弯曲，并且该线的弯曲可用于确定由结构光投射器(例如上述的结构光投射器170)照亮的所有位置上的表面高度。如图2D所示，结构光投射器170被固定在光传输装置120上并与其一起旋转。在这种情况下，结构光投射器170由柯勒照明系统组成，其中幻灯片由光源照亮继而把所述幻灯片的图像投射到所述样本或背景上。投射角度要足够大以便得到在线内的充分的“弯曲”来实现空间分辨率，但是足够小不会出现大的阴影。

然后处理流程510以没有结构光来成像所述样本来进行(513和图4D)。这可以对应于如上所述的照相图像捕获(524)。可以根据2D傅立叶变换来确定在背景和样本上所有点处每条线的相移。

然后把背景数据转换为包裹相位(514和图4E)。这里，在计算包裹相位之前对所述背景数据进行傅立叶变换和滤波。类似地，通过傅立叶变换和滤波所述样本数据、以及计算所述样本数据的包裹相位来把所述样本数据转换为包裹相位(515和图4F)。

然后计算所述样本的表面形貌(516和图4G)。在这种情况下，这通过“展开”所述相位图来执行。为此任务几个展开算法对本领域技术人员来说是可以采用的。例如，可以根据使用傅立叶轮廓测定法技术来确定在所述图像上所有点处每条线的相移。使用这些方法，对所述边缘数据(图4D)进行2D快速傅立叶变换(FFT)来确定在所述图像中各处的线的相移(图4F)。因为对于典型的物体所述相位将移动2□的多倍，所以如图4F所示所述相位显示2□跳跃。“展开”这些相位跳跃以便确定实际的表面。

然后可以根据不同的视图和位置重复(517)上述过程(511-516)。根据多视图来成像样本提供了附加信息，其有助于这里所描述的技术提供更精确的3D表面再现。然后，把根据在所述3D成像系统中的每个视图所获得的多个图像或部分表面一起对准以形成完整的3D表面(518和图4H)。可以通过使用非线性最小二乘拟合技术来完成对准，以使在连接的两个表面上的网格元之间的距离最小化。典型地，所述表面应该具有相当接近于最后对准位置的开始方向。换句话说，用这种方法只可以提供对表面位置的精调。例如，另一对准技术是在图像中提供某种绝对参考线或基准，其给出了任何部分表面相对于所述载物台的绝对位置。如果每个表面的绝对定位足够精确，那么可以跳过上述的非线性拟合法。

图5D说明了依照本发明具体实施例的、用于把表面光发射数据转换为在物体的体积内光源分布的重构的处理流程540。依照本发明的过程可以包括多达几个附加步骤，这里不描述或说明它们以便不致使本发明变得不明显。

处理流程540假定拥有了所述样本表面表示和拥有了来自所述样本的光发射数据。处理流程540从把表面表示划分为表面元开始(541)。表面元的数目将依照样本表面区域和所要求解的精确性来变化。表面元的数目应该足够大以便捕获穿过所述表面的光子密度的细节和变化。例如，在大约100到大约10,000之间的表面元可以适于老鼠。

此外，处理流程540把所述样本内部体积分为体积元(542)。在一个实施例中，每个体积元被认为在其中心包含点光源。然后，这些体积元的立体网格定义了许多点源，这些电源近似于实际光源分布。有时，增加在光源附近的所述立体网格的密度以在此空间内提供增加的信息，而在那些没有发生感兴趣动作(没有光产生或传输)的空间内降低所述立体网格的密度。此外，如下面相对于循环556将要描述的，所述体积元大小可以在获得解期间依照各种自适应建网技术来变化。例如，初始体积元大小可以从大约0.1cm³到大约1cm³变化，并且接近于所述源的体积元的最终体积元大小可以减少从大约1*10^-3cm³到大约1*10^-2cm³。在特定的例子中，初始体积元大小可以大约是1cm³，而接近于所述源的体积元的最终体积元大小可以降低为大约8*10^-3cm³。

然后处理流程540在表面元和体积元之间建立关系(544)。在一个实施例中，重构方法使用在表面上的源发射强度和光子密度之间的线性关系。在具体实施例中通过格林函数来描述该线性关系。所述格林函数包含所有关于在所述样本内光子传输的信息，其包括在体积中不均匀性和边界上内部反射的效果。所述格林函数描述了在所述样本内从所述源分布中每点到所述样本每个表面元内的光子传输。

当把在所述样本内的介质假定为均匀的或为其建立模型是均匀之时，格林函数的一个有益的形式是简化的近似，其中把所述样本的表面局部地对待为方向与连接体积元中心和表面元的线垂直的平面界面。使用分流边界条件在表面上的光子密度是在半无限平板中点源的解析解。因为所述解只是在所述体积元和所述表面之间距离的函数，所以在最小的计算代价下可以对所有的体积元和表面顶点对计算简化的格林函数。

利用格林函数G_ij描述的在每个体积元中的源强度和在每个表面元的光子密度之间的线性关系，可以通过来自所有体积元的成分之和来估计在第j个表面元处的光子密度：

${\mathrm{\ρ}}_i\≅\underset{i}{\mathrm{\Σ}}{G}_{\mathrm{ij}}{S}_i---\left(4\right)$

其中下标i计数所述体积元并且S_i是在第i个体积元内的点源的强度值(光子/秒)。

在诸如组织之类的混浊介质中的光传输是由散射来占统治地位的，并且实际上其基本上是扩散的。扩散传输的条件是散射系数μ_s大于吸收系数μ_a，以致在光子密度方面的变化在散射事件之间很小。可以由扩散方程来表示由在均匀介质中源功率密度U_i生成的光子密度：

$D{\▿}^2\mathrm{\ρ}-{\mathrm{\μ}}_a\mathrm{c\ρ}=-U_i\left(\underset{\‾}{x}\right)---\left(5\right)$

其中扩散系数D是

$D=\frac{c}{3({\mathrm{\μ}}_A+{\mathrm{\μ}}^{\′}s)}---\left(6\right)$

在这种情况下，格林函数是在受到由所述样本表面施加的边界条件下方程(5)的解。对于为样本建立模型是均匀的情况而言，格林函数的有益的近似解使用在所述表面边界上的平面近似。

图3B说明了示出此平面近似的示意图。通过第j个表面元的中心，垂直于连接在x_i处的点源和在x_j处的表面元的线画出平面边界312。根据分流边界条件，在平面近似中光子密度是在由所述平面边界定义的半无限平板中在x_i处的点源的解。在这种情况下，把所述样本表面局部地对待为方向与连接在x_i处的源和在x_j处表面元的线垂直的简单的平面边界。特别地，把边界条件简化为平面边界情况，尽管所述边界方向可以改变每个表面元。

使用分流边界条件，此简化格林函数是在半无限平板中点源的解析解：

$G_{\mathrm{ij}}=\frac{1}{2\mathrm{\πD}}\{\frac{\exp (-{\mathrm{\μ}}_{\mathrm{eff}}{r}_{\mathrm{ij}})}{r_{\mathrm{ij}}}-\frac{1}{z_b}\exp (r_{\mathrm{ij}}/z_b)E_1[({\mathrm{\μ}}_{\mathrm{eff}}+\frac{1}{z_b})r_{\mathrm{ij}}\left]\right\}---\left(7\right)$

这里r_ij＝|x_j-x_i|，E₁是一阶指数积分并且

μ_eff＝[3μ_A(μ_A+μ′s)]             (8)

$z_b=\frac{2D}{c}\frac{1+R_{\mathrm{eff}}}{1-R_{\mathrm{eff}}}---\left(9\right)$

在刚刚描述的简化模型中，所述简化的格林函数只依赖于在所述体积元和表面之间的距离。不必要使用诸如简化近似之类的解析形式来定义格林函数。

本发明不依赖于诸如上述近似之类的解析形式。在另一实施例中，查找表可以定义所述格林函数。可以通过先前测量的在样本(或近似为基本上等于当前样本的类似样本)中的光子传输，或通过使用诸如蒙特卡罗或有限元建模之类的技术计算模拟来创建查找表。这种特殊方法用于由诸如动物研究对象之类的不均匀介质组成的样本。在这种情况下，来自方程8的光学属性μ_a和μ_s现在具有空间相关。

上述的平面边界近似对于具有大曲率半径的平滑表面以及对于其中吸收系数不太小(□_a＞0.1cm^-1)的情况用起来最好。上述平面近似技术的优点是计算上方便于解答具有任意复杂边界的扩散方程，例如老鼠。具有更多结构的区域，例如老鼠的头或四肢，可以受益于更精确的边界模型。使用有限元建模代码来计算格林函数是一种获得更精确边界模型的选择。诸如来自PDE Solutions公司的Flex PDE之类的有限元代码例如可以被使用。另一选择将是把平面表面近似依照曲率扩展为一阶，其可以允许继续使用G_ij的解析表达式。

一旦确定了格林函数，通过解出线性方程组来获得所述重构，所述线性方程组使在表面上的光子密度与在所述物体内的源分布相关。然后，处理流程540通过求解所有的内部体积元来进行(546)。更具体地说，在给定了上述建立模型的情况下，所述重构技术求解使在表面上的光子密度与在所述样本内的源分布相关的线性方程组。因此，一旦确定了所述格林函数，可以对每个体积元-表面元对计算以便获得所述线性方程组(公式4)。重构方法的最后的步骤是求解线性方程组(4)，以得到源强度为S_i。返回参考公式(4)，因为□已知，并且如上所述可以确定G_ij，那么重构技术可以求解S_i。典型地没有线性方程组的精确解，因为多个点源只是实际源分布的近似。那么一个适当的重构是所述线性方程组的最优近似解。

在一个实施例中，把线性方程组转换为最优化问题以找到近似解。在这种情况下，使用代价函数来定义所述最优化问题。可以通过源强度之和来表示一个适当的代价函数：

$\mathrm{Cost}=\underset{j}{\mathrm{\Σ}}{S}_j---\left(10\right)$

所述代价函数受到一个或多个线性约束条件的影响。第一个适当的约束条件集是源强度是正定的：

S_i≥0                 (11)

第二个适当的约束条件集是基于所述线性方程组(4)的线性不等式组。它要求遍及所有的所述体积元求和的光子密度的基值小于在每个表面元处的光子密度的测量值：

$\underset{j}{\mathrm{\Σ}}{G}_{\mathrm{ij}}{S}_j\≤{\mathrm{\ρ}}_i---\left(12\right)$

在具体实施例中，通过使得受到约束条件(11)和(12)影响的代价函数(10)最大化来找到源强度S_i的最优解。

然后，处理流程540在解集中求出当前的最优解(550)。可以使用线性分析技术来解决上述最优化问题。一个尤其高效的技术是单纯形(SIMPLEX)法。

可以评价所述解的质量(552)。在一个实施例中，所述估计测量在所观察的光子密度和所计算的光子密度之间的差值。例如，可以使用“χ²”标准：

${\mathrm{\χ}}^2=\underset{j}{\mathrm{\Σ}}{\left[\frac{{\mathrm{\ρ}}_j-\underset{j}{\mathrm{\Σ}}{G}_{\mathrm{ij}}{S}_i}{{\mathrm{\ρ}}_j}\right]}^2---\left(13\right)$

χ²的值测量在所述样本表面上所观察的光子密度ρ_i和所计算的光子密度

之间的差值。

可以把基于表面光发射数据的内部3D光源分布数据的确定嵌入到一个或多个迭代循环中。迭代可以允许体积网格的细分，并且可以允许搜索表面元的最佳采样和最佳加权函数配置。一般说来，通过降低线性方程组的大小来改善重构的收敛性。在一个实施例中，本发明的重构技术反复地改变下列中的一个或多个：体积元配置、表面元配置和/或代价函数配置。

迭代还可以帮助获得解。单纯形法迅速地收敛于解，并且由此具有作为计算上高效的算法的优点。为了帮助单纯形算法找到全局最大值，并因此得到最优解，当系统地改变在所述算法中的关键元时可以重复单纯形计算多次(554)。这可以通过使用改变体积元数目的自适应建网或自适应分格以及通过使用所述表面元的子集来实现。

循环556(图5D)说明了当改变体积元的数目时重复所述求解步骤。典型地，在该问题中减少体积元的数目而维持在所述源附近的体积元的高密度是有益的。这可以通过使用自适应建网来实现。在一个适当的自适应建网应用中，最初把粗略的体积元网格施加到整个样本体积上，并且找到了当前的解(544、546和550)，产生S_i的初始解。接下来精化(即细分)具有源强度大于零(S_i＞0)的体积元，并去除那些源强度等于零(S_i＝0)的体积元。然后，可以重复地迭代解的获得和体积元网格细分，生成局部化于源区域的高密度体积元网格。在每次迭代期间，估计当前解的质量(552)。在具体实施例中，继续该迭代过程直到进一步的细分对χ²的值不会产生显著的下降。

通过改变用于获得该三维表示的表面元的数目N_S可以获得附加的迭代改进(循环555)。使用表面网格的表面元的子集降低了在该问题中的约束条件的数目，这可以简化解的计算。表面元的数目可以用来均匀地样所述表面。在这种情况下，处理流程540迭代对应于在不同的密度处采样所述表面元网格的N_S的不同值，并且使用质量评价(552)来确定在不同N_S值之间的最优解。例如，如果对于小老鼠来说表面元的数目在大约100和大约300个表面元之间，则迭代步长大小在大约10和50之间也可以是适当的。

另一迭代循环改变代价函数来改善解的获取。例如，可以修改所述代价函数以包括加权系数W_i：

$\mathrm{Cost}=\underset{i}{\mathrm{\Σ}}{S}_i/W_i^r,W_i=\underset{j}{\mathrm{\Σ}}{G}_{\mathrm{ij}}---\left(14\right)$

所述加权系数W_i是经由整个表面的光子密度的第i个体积元的基值。指数γ调整内部体积元和那些接近于所述表面的体积元的代价函数的相对基值。当γ＝0时，那么所述内部体积元具有相对较大的权重。当γ＝1时，接近表面的体积元具有较大的权重。当改变γ时可以迭代处理流程540来搜索其中所述源接近和远离所述表面的解。例如，对于γ从0到1的范围可以按照大约0.01到大约0.2来改变步长大小。在具体实施例中，对于γ从0到1的范围可以改变步长大小大约0.05。可以再次使用质量评价(552)来识别所述最优解。

在具体实施例中，把单纯形计算嵌入在由两个最优化循环组成的搜索例行程序中。内循环使用自适应分格来最优化所述体积网格的配置。外循环搜索参数N_S和γ的最优值。使用χ²标准来检验所述搜索算法的每个步骤。把全局解定义为最小化χ²的一个解。

图5E和5F表示用于在仿真老鼠(例如具有嵌入式光学纤维的橡胶老鼠)内的一个点状光源的重构结果。图5E说明了以两个不同的角度在仿真老鼠702a和702b内的重构源分布。示出的表面表示以最大源强度的一半的等值面的轮廓。图5F示出了直接在仿真老鼠702的表面下面的光子密度的图像(例如708a和708b)。第一行706a是所测量的光子密度的对数。第二行706b是根据所述重构源分布而计算的模拟光子密度的对数。第三行706c是第一和第二行的差值。

尽管为了加快处理已经用许多简化来描述处理流程500、510、520和540，但是本发明并不受这些简化计算方法的限制。例如，可以在没有上述许多简化的情况下计算所述格林函数，甚至以增加计算需要为代价。此外，虽然处理流程500、510、520和540基于在图1-2E中说明的系统描述了获得来自所述样本的光发射测量的特定方法，但是本发明的3D重构技术不局限于如何获得2D光发射数据或任何特定设备的使用。实际上，可以使用不同于上述相对于图1-2E的成像设备，例如在共同拥有的待决的专利申请序号为09/905,668中描述的成像设备。此外，可以从独立源获得2D光发射数据并且将其存为在计算机28内的数据，而不必作为经由补充的或局部成像系统成像的结果来生成。

此外，尽管至此已经相对于生物发光源描述了本发明，但是可以使用本发明来获得任何类型光源的3D重构，所述任何类型光源包括一个或多个荧光源。对于荧光源，可以改变在所述表面元和体积元之间的关系(544)。例如，可以把在线性方程组中的格林函数乘以第二格林函数，所述第二格林函数描述了激发光从样本表面到体积元的传输。在特定的近似中，可以把在线性方程组(4)中的格林函数修改为两个格林函数的乘积：

G_ij＝G_i ^EG_ij ^F              (15)

第一格林函数G_i ^E描述了激发光从在样本表面上的激发源到第i个体积元的传输。第二格林函数G_ij ^F描述了荧光从第i个体积元到第j个表面元的传输。可以在均匀介质情况下根据诸如上述简化的近似之类解析表达式或在不均匀介质情况下根据查找表来确定这两个格林函数。所述激发和荧光典型地处于不同的波长，并且从而所述荧光不激发额外的荧光。线性方程组(4)仍然有效，并且可以如上概述使用处理流程540来确定荧光源分布。

本发明还非常适合于处理在输入图像中的噪声。任何实际的成像系统将具有有限的动态范围并且受到噪声源的影响。这可以损害在所述重构中用作输入数据的图像保真度，并且从而可以降低3D表示的质量。为此，本发明还可以包括设计来当所述输入图像是有噪声的和/或具有有限动态范围时改善其性能的技术。

当成像来自嵌入在组织中的源的生物发光时有限的动态范围尤为明显，因为光发射亮度在穿过样本表面时典型地在许多数量级上改变。如果所述成像照相机强加有限的动态范围并把最高亮度区域设置为照相机的上限，那么可能在图像中存在其中发射亮度降到所述动态范围底部之下的区域。将把所述图像的这些区域接收为噪声；并且对应于“固有噪声”。

这里提供的重构技术把表面光发射的图像映射到三维样本上以确定在表面下光子密度的值。由于在所述图像中的噪声，重构技术只映射那些在固有噪声之上的图像的区域。这可以通过为每个图像定义阈值来实现，所述阈值是峰值亮度的一小部分。在一个实施例中，所述阈值与峰值亮度和照相机的动态范围有关。在具体实施例中，所述阈值可以比由峰值亮度除以动态范围所表示的数目大。例如，如果照相机的动态范围是1000并且图像中的峰值亮度是I_P，那么所述阈值的适当值可以是I_P/500。只把在所述阈值之上的图像的区域映射到所述样本上；而不使用那些在所述阈值之下的区域。结果，在完成映射之后会有空的样本表面的区域存在，也就是说，这些区域没有与所述表面网格相关联的光子密度值。典型地，这些将趋向于是最远离生物发光源的表面部分。

处理流程540解决下列最优化问题：根据所计算的光子密度小于在表面处处测量的光子密度的约束条件来最大化总的源强度。然而，处理流程540还可以解决表面空的部分。这是因为当光子密度是未知的时，对于那些接近空的区域的体积元去除在源强度上的约束条件。因此，可能没有充足的信息来确定在那些体积元中的源强度。处理流程540的变型由消除在体积方格中的体积元组成，所述体积元没有足够地被现有的光子密度信息约束。这可以通过把体积元的基值与在空的和非空的表面区域中的光子密度相比较来实现。例如，如果P和Q分别表示由空的和非空的区域组成的表面元集，那么用于去除来自所述体积方格的第i个体积元的标准是：

$\underset{j\&Element;Q}{\mathrm{\&Sigma;}}{G}_{\mathrm{ij}}<\mathrm{\&kappa;}\underset{j\&Element;P}{\mathrm{\&Sigma;}}{G}_{\mathrm{ij}}---\left(16\right)$

常数K可以具有在1-10范围内的值。在形成初始体积方格(542)和每次迭代期间把标准(16)应用于每个体积元，如果使用迭代的话。

除具有有限的动态范围之外，真实图像可以包含噪声。一个最小噪声成分可以通过散粒(泊松)噪声来表示，所述散粒噪声与在每个图像像素中离散的许多光子的采样相关联。当散粒噪声占优时，在所述图像中的信噪比随亮度的平方根而改变。这样我们就可以假定最高亮度的区域，或相当于具有最高光子密度的表面的那些区域应该具有相对最小量的噪声。为了解决此噪声，可以改变解质量的评价(552)。例如，可以修改表达式(13)来解决噪声。更具体地说，可以用系数光子密度ρ_i来加权总和中的每一项以变为：

${\mathrm{\&chi;}}^2=\underset{j}{\mathrm{\&Sigma;}}\frac{1}{{\mathrm{\&rho;}}_j}{[{\mathrm{\&rho;}}_j-\underset{i}{\mathrm{\&Sigma;}}{G}_{\mathrm{ij}}{S}_i]}^2---\left(17\right)$

典型地由适当的处理器或基于计算机的设备来实现本发明的重构技术。参照图6，示例性的计算机系统350包括中央处理器(CPU)352、只读存储器(ROM)354、随机存取存储器(RAM)356、扩展RAM358、输入/输出(I/O)电路360、显示部件362、输入装置364和扩展总线366。计算机系统350还可以选择性地包括诸如磁盘驱动器单元之类的大容量存储单元368或诸如闪速存储器之类的非易失存储器和实时时钟360。

不管计算机系统350的配置是怎样的，其可以使用一个或多个存储器或存储器模块，所述存储器或存储器模块被配置为存储用于获得位于样本内光源的三维表示和这里描述的本发明的其它功能的程序指令。还可以把这种一个或多个存储器配置为存储数据结构、成像数据或这里描述的其它特定的非程序信息。

因为可以使用这种信息和程序指令来实现这里描述的系统/方法，所以本发明涉及包括用于执行这里描述的各种操作的程序指令、状态信息等的机器可读介质。机器可读介质的例子包括但不局限于：诸如硬盘、软盘和磁带之类的磁介质；诸如CD-ROM盘片之类的光学介质；诸如可光读的盘片之类的磁光介质；以及专门配置为存储并执行程序指令的硬件装置，例如只读存储器装置(ROM)和随机存取存储器(RAM)。还可以用通过诸如无线电波、光线路、电线路等之类的适当介质传播的载波来使本发明具体化。程序指令的例子包括例如由编译器生成的机器代码和包含高级编码的文件，所述高级编码可以由计算机使用解释器来执行。

应当记住，尽管为了使论述便于进行这里相当详细地论述了计算机系统350，但是可以使用各种适当的计算机实现技术来实施本发明。一般说来，可以使用任何适当的计算机系统以便获得位于样本内光源的三维表示。此外，可以经由诸如局域网(LAN)、广域网(WAN)或诸如因特网的全球计算机网络之类的计算机网络来实现这里公开的发明的重构技术。在后一种情况中，可以把所述重构技术至少部分地实现为可下载的计算机软件和数据(例如诸如来自Sun Microsystems公司的JAVA^TM小应用程序(applet)的小应用程序)。可下载的计算机软件和数据可以保存在网络上的一个或多个服务器中，其可由有能力并被授权这种访问的任何客户端计算机或末端来访问。在现有技术中网络计算技术和实现是众所周知的，并且为了简便起见这里没有非常详细地进行讨论。

虽然按照几个优选实施例已经描述了本发明，但是存在属于本发明范围的修改、置换和等价物，为了简便起见省去了这些。例如，可以在一个或多个不同的波长获得并分析图像。因此意图是应该参考附加权利要求来确定本发明的范围。

Claims (52)

1.一种用于获得在样本内光源分布的三维表示的方法，所述方法包括：

提供来自发自所述样本表面的光的表面光图像数据，所述光来源于位于在所述样本内光源的分布；以及

根据所述表面光发射数据重构在所述样本内光源分布的三维表示。

2.如权利要求1所述的方法，还包括获得所述样本的表面形貌表示。

3.如权利要求2所述的方法，还包括把所述表面形貌表示划分为一组表面元。

4.如权利要求3所述的方法，其中每个表面元近似为平面。

5.如权利要求3或4所述的方法，还包括在所述样本内创建一组体积元。

6.如权利要求5所述的方法，其中为每个体积元建立模型以便使在它的中心包含一个点光源。

7.如权利要求6所述的方法，其中用一组点光源来近似所述光源分布的三维表示。

8.如权利要求5至7中任何一个所述的方法，还包括把所述表面光图像数据转换为刚好在所述样本表面内的光子密度。

9.如权利要求8所述的方法，其中在一个给定的体积元中的光源发射强度和刚好在表面元内的光子密度之间存在线性关系。

10.如权利要求5至9中任何一个所述的方法，还包括定义一个代价函数和一组约束条件，用于获得光源分布的三维表示的解。

11.如权利要求10所述的方法，其中所述代价函数与在所述样本中每个点源的源强度之和有关，所述约束条件包括下列条件：(i)所述源强度是正定的；并且(ii)由点源分布生成的、在物体表面作为结果的光子密度处处小于所测量的表面光子密度。

12.如权利要求11所述的方法，其中获得三维表示使得受所述约束影响的所述代价函数最大化。

13.如权利要求10至12中任何一个所述的方法，其中在数学上用线性方程系统描述所述代价函数和约束条件，并且使用单纯形法来获得源分布的三维表示的解。

14.如权利要求10至13中任何一个所述的方法，还包括在所述代价函数中包括加权系数，可以改变所述加权系数以生成所述源分布的三维表示的一组解。

15.如权利要求10至14中任何一个所述的方法，还包括改变表面元的数目以生成所述源分布的三维表示的一组解。

16.如权利要求5至15中任何一个所述的方法，还包括改变下列中之一：a)体积元的数目；和b)体积元的配置，以便生成所述源分布的三维表示的一组解。

17.如权利要求5至15中任何一个所述的方法，还包括通过计算每个解的表面光发射并选择使在所计算和所测量的表面发射之间的差值最小化的解来最优化光源分布的三维表示。

18.如权利要求16或17所述的方法，其中改变体积元的数目和改变体积元的配置这二者均包括自适应建网。

19.如权利要求18所述的方法，其中所述自适应建网增加用于描述所述光源的三维表示的体积元的数目。

20.如权利要求19所述的方法，其中所述自适应建网去除具有零光源强度的体积元。

21.如权利要求5至20中任何一个所述的方法，其中用格林函数描述光在所述样本内从给定的体积元到给定的表面元的传输。

22.如权利要求21所述的方法，其中把所述格林函数定义为光在均匀的半空间中扩散的解，所述均匀的半空间具有垂直于把体积元和表面元相连接的线的平面边界。

23.如权利要求21或22所述的方法，其中所述样本内部近似为不均匀的。

24.如权利要求21至23中任何一个所述的方法，其中在查找表中定义所述格林函数。

25.如权利要求21至23中任何一个所述的方法，其中使用蒙特卡罗模拟或有限元建模来计算所述格林函数。

26.如权利要求1至25中任何一个所述的方法，其中所述光源由生物发光或荧光发射组成。

27.如权利要求1至26中任何一个所述的方法，还包括把噪声阈值施加到表面光图像数据上。

28.如权利要求27所述的方法，其中所述噪声阈值与在所述表面光图像数据中的峰值亮度和在所述表面光图像数据中的动态范围之一有关。

29.如权利要求28所述的方法，其中所述噪声阈值与除以在所述表面光图像数据中的动态范围的所述表面光图像数据中的峰值亮度有关。

30.如权利要求27所述的方法，其中把所述表面表示划分为一组表面元，并且去除所有具有在所述噪声阈值之下的表面发射的表面元。

31.如权利要求1至30中任何一个所述的方法，其中所述样本是动物，所述光源发射穿过动物组织的光。

32.如权利要求31所述的方法，其中所述动物组织近似于均匀的。

33.如权利要求1至32中任何一个所述的方法，其中所述样本具有复杂的边界。

34.如权利要求1至33中任何一个所述的方法，还包括生成所述光源的多个可能的三维表示，并且所获得的所述光源的三维表示是最适合于所测量的表面光图像数据的表示。

35.如权利要求1至35中任何一个所述的方法，还包括把所述样本放置在纳入成像室的载物台上，所述成像室与照相机相耦合，所述照相机被配置为捕获所述载物台上样本的图像。

36.如权利要求34所述的方法，还包括：

把所述载物台移动到所述成像室中的第一位置；以及

使用所述照相机从第一位置捕获所述样本的第一图像集。

37.如权利要求35所述的方法，其中第一图像集由荧光图像、结构光图像和照相图像组成。

38.如权利要求34所述的方法，还包括：

把所述载物台移动到在所述成像室中的一个或多个其它位置，其中所述其它位置相对于与所述照相机相关联的固定基准的角度与相对于第一位置的角度不同；以及

使用所述照相机从所述其它位置捕获所述样本的附加的图像集。

39.如权利要求38所述的方法，其中获得所述表面表示包括：根据包括在一个或多个结构光图像内的结构光数据来构建所述样本的形貌表示。

40.如权利要求34所述的方法，其中以一个或多个不同波长来获得所述表面光图像数据。

41.一种用于获得位于样本内光源的三维表示的成像系统，所述成像系统包括：

具有一组壁的成像室，所述壁围住一个内腔，所述成像室包括：

配置为安置照相机的照相机架；

可移动载物台设备，包括一个传送机构和一个载物台，所述载物台被配置成支撑所述内腔之内的样本，所述载物台与所述传送机构耦合以便把所述样本移动到在所述内腔中多个位置之一；

用于把发自所述样本表面的光传送到所述照相机的光传输机构；以及

处理器，被设计或配置成提供来自发自所述样本表面的光的表面光图像数据，以便根据所述表面光图像数据来获得内部光源分布的三维表示，其中所述光来源于在所述样本内光源的分布。

42.如权利要求41所述的系统，其中还把所述处理器配置为根据由所述照相机接收的结构光数据来构建所述样本的表面形貌表示。

43.如权利要求41或42所述的系统，其中所述传送机构能够在二维空间移动所述载物台。

44.如权利要求41至43中任何一个所述的系统，其中所述光传输装置绕固定轴旋转。

45.如权利要求41至44中任何一个所述的系统，还包括配置为把结构光传送到所述样本上的结构光源。

46.如权利要求45所述的系统，其中所述结构光源，在结构光到达所述样本之前，把所述结构光传送到镜子上。

47.如权利要求41至46中任何一个所述的系统，其中所述光源由生物发光或荧光发射组成。

48.一种包括计算机可读介质和经由所述计算机可读介质提供的程序指令的计算机程序产品，所述程序指令包括用于获得在样本内光源分布的三维表示的重构指令，所述重构指令能提供来自发自所述样本表面的光的表面光图像数据，并且根据所述表面光图像数据获得所述光源分布的三维表示，其中所述光来源于在所述样本内光源的分布。

49.如权利要求48所述的计算机程序产品，还包括用于显示所述光源分布的三维表示的指令。

50.如权利要求48或49所述的计算机程序产品，还包括用于把由照相机接收的光图像数据映射到所述样本的表面表示上的指令。

51.如权利要求50所述的计算机程序产品，还包括用于把所述表面表示划分为一组表面元的指令。

52.如权利要求48至51中任何一个所述的计算机程序产品，其中所述光源由生物发光或荧光发射组成。

Images