CN101692971B - 非接触式光学断层成像方法 - Google Patents

Abstract

本发明公开了一种非接触式光学断层成像方法，解决了现有技术中无法实现的大孔径镜头下非接触式光学断层成像问题。其实现步骤为：先用多模态光学分子影像系统采集多角度光学图像和生物体表面三维形状及解剖结构信息；然后重建生物体表面三维能量，包括：(a)光学图像的预处理和灰度值到能量值的转换；(b)光学图像与生物体表面三维形状信息的配准；(c)虚拟探测面位置和像素大小的确定以及CCD相机到虚拟探测面的能量映射；(d)计算可见度因子；(e)重建表面三维能量分布；最后反演体内光源并三维立体化显示。本发明具有能实现大孔径镜头下表面三维能量重建和大孔径镜头下非接触式断层成像的优点，可用于非接触式光学分子成像领域。

Description

非接触式光学断层成像方法

技术领域

本发明属于光学分子成像技术领域，涉及非接触式光学断层成像技术，特别是一种由非接触式CCD相机测量的多角度二维光学图像到生物体内部光源间接信息的成像方法。

背景技术

近年来光学分子成像技术凭借其高灵敏特性和低廉的价格以及易操作性在小动物研究中得到了广泛的应用。为得到直观的三维可视化成像效果，光学分子成像技术新发展起来的光学断层成像技术成为新的应用研究热点，其主要目标是根据生物体表面分布的光学信号反演生物体内的光源间接信息从而获得生物体内组织病变信息。在光学断层成像技术中，激发荧光断层成像和生物发光断层成像是两种主要的成像模态。

光学断层成像技术重建生物体内的光源依据的关键数据是生物体表面的光学分布信息，而要得到生物体表面的光学分布信息根据探测器的测量方法不同，主要有接触式和非接触式两种。接触式测量方式采用光纤探头接触被测物体表面获得透射扩散光通量，其优势在于能够直接采集到生物体表面的光学信号信息；但是由于光纤数量有限，采样只是得到在被测物体表面的有限点导致空间分辨率很低。基于高灵敏CCD的接触式测量方式是采用固定几何形状的成像腔同时采用光纤传输，这种方法可以得到边界信息，但是同样存在空间分辨率低的问题。非接触式测量方式采用高灵敏度的CCD和镜头采集生物体表面的光学信号，由于CCD的高分辨率可以捕获生物体表面大规模的光学数据，达到亚毫米的空间分辨率，同时大大降低了接触式采集装置的复杂性。

在非接触式光学断层成像中，光的传输包括两个部分：首先光从被测生物体内部光源处出发，在生物体内传输到生物体与空气的分界面；此后从生物体表面逃逸，在自由空间中传输，最后被探测器接收。因此，非接触式光学断层成像重建生物体内部光源时，首先要获得复杂生物体表面的光学信号分布信息，然后再根据该信息采用合适的重建算法进行体内反演，获取生物体内光源。然而，由于光逃逸出生物体表面在自由空间中传输的朗伯源特性以及自由空间中复杂光学系统的介入，使得2D探测器数据到3D复杂生物体表面的能量重建是一个复杂而又棘手的问题。2003年，Ripoll首次提出了自由空间中的近红外光传输理论，讨论了近红外光从生物体表面透射光通量到非接触式探测器面上光强之间的映射关系，参见：Jorge Ripoll.“Free-space propagation of diffuse light：Theory and Experiments”[J].Physical Review Letters，91(10)，2003。随后，Ripoll等把该理论应用于FMT的研究，他们综合了光在生物体内部和自由空间的传输过程建立了统一的模型，并实现了荧光光源的重建，参见：Nikolaos Deliolanis.“Free-space fluorescence molecular tomographyutilizing 3600geometry projections”[J].Optics Letters，32(4)，2007。但是，在这些研究中，自由空间中光的传输理论仍旧局限在正向问题，即从生物体表面到探测器面的映射关系。白净等人在Ripoll工作的基础上，考虑了光学系统中镜头的空间透视效应，提出了一种简化的自由空间光传输模型，但是同样也只能用来描述光从生物体表面到探测器面的映射问题，参见：姚俊杰.“非接触式近红外荧光断层成像中光在自由空间中的传播模型及验证”[J].红外与毫米波学报，27(5)，2008。由于生物体表面某一点的出射光对探测面上的很多点上的光强都有贡献，使得生物体表面上的光源点与探测面上的接收单元不仅存在一对一的关系，也可能存在一对多的关系，也就是说探测面上每个像素的光强是物体表面所有源点的加权求和。因此，并不能根据已有的、从生物体表面到探测面的能量映射关系的简单求逆，来确定从探测器面到生物体表的空间和能量映射关系。2007年，Ripoll等人结合小孔径镜头成像情况，简化2003年提出的自由空间光传输模型为针孔光传输模型，建立起表面点与探测器平面点之间的一一对应关系，从而可以直接逆向求解表面能量分布，参见：Heiko Meyer.“Noncontact optical imaging in mice with full angularcoverage and automatic surface extraction”[J].Applied Optics，46(17)，2007。然而，这种方法只适用于小孔径镜头的情况，对于稍大孔径镜头或者大孔径镜头都是不适用的，因而无法实现大孔径镜头情况下的非接触式光学断层成像问题。在光学分子成像过程中，由于生物体表面逸出的光信号极其微弱，需要采用大孔径镜头增加采集的光信号的强度。因此，实现一种大孔径镜头下由非接触式CCD相机测量的多角度二维光学图像到生物体内部光源信息的成像方法就成为非接触式光学断层成像技术中亟待解决的问题。

发明内容

本发明的目的在于克服上述已有非接触式光学断层成像技术中存在的不足，提出一种新的非接触式光学断层成像方法，以实现大孔径镜头下非接触式CCD相机测量的多角度二维光学图像到生物体内部光源信息的断层成像过程。

实现本发明目的技术关键是结合光学系统的简化理论，利用朗伯余弦定律和基本辐射亮度定理，进行任意形状生物体表面的三维能量分布重建，同时结合光源重建算法反演内部光源信息，实现大孔径镜头下的非接触式光学断层成像过程。整个成像步骤包括如下：

A.数据采集步骤：利用非接触式光学断层成像系统采集多角度的二维光学图像信息，同时利用微计算机断层成像或微核磁共振成像获得生物体的三维表面形状和内部解剖结构信息；

B.表面三维能量重建步骤：

(B1)对采集的多角度光学图像数据依次进行背景噪声去除、感兴趣区域提取和坏点补偿，以及将光学图像数据灰度值转换到对应的能量值；

(B2)将步骤A中使用的非接触式光学断层成像系统的成像镜头简化为薄透镜模型，并根据薄透镜成像规律计算CCD相机对应的虚拟探测面的位置和每个像素的实际尺寸；

(B3)将采集的二维光学图像信息与生物体表面三维形状信息进行配准，以获得各个角度二维光学图像数据中心点相对于生物体表面三维形状信息坐标系的实际物理坐标；

(B4)利用光学成像系统中像的基本辐射亮度定理，同时结合镜头的透过率，将步骤B1中获得的CCD相机上的能量分布映射到对应的虚拟探测面上；

(B5)离散步骤A中得到的生物体三维表面，记录其表面面片总数为N，N≥500，并记录各个表面面片的表面法向信息和三个顶点坐标信息；

(B6)取第i个表面面片的中心点，1≤i≤N，判断其与各个角度虚拟探测面上的所有像素点之间的空间位置对应关系，同时采用表格搜索策略计算其与CCD相机上任意一个像素点之间的可见度因子；

(B7)根据可见度因子的具体取值，应用下述公式计算第i个表面面片上的能量重建值：

$P\left(r\right)=\frac{1}{\mathrm{\π}}{\∫}_{\mathrm{\Ω}}E\left(r_d\right)T(r_d,t)\mathrm{d\Ω}$

式中，r是第i个表面面片的中心点；

P(r)是表面点r处重建的功率；

E(r_d)是CCD相机点r_d处的功率密度；

dΩ是CCD相机空间；

T(r_d，r)定义为光在自由空间的传输函数：

$T(r_d,r)=\frac{1}{\mathrm{\τ}\left(f\right)}\mathrm{\ξ}(r_d,r)\frac{1}{{|r-r_d-\frac{l}{\cos \mathrm{\θ}}{s}_{r_d\→r_{\mathrm{vd}}}|}^2}\cos {\mathrm{\θ}}_s\cos {\mathrm{\θ}}_d\mathrm{dS}\left(r\right)$

其中，τ(f)是非接触式光学断层成像系统中镜头的透过率，ξ(r_d，r)是表面点与CCD相机像素点之间的可见度因子，l是真实实验中的物像距离，cosθ是CCD相机上像素点r_d到虚拟探测面上像素点r_vd之间的连线与主光轴的夹角余弦，

是CCD相机上像素点r_d到虚拟探测面的像素点r_vd的单位方向向量，cosθ_s是表面点的法向向量与表面点到虚拟探测面点之间连线的夹角余弦，cosθ_d是虚拟探测面点的法向向量与虚拟探测面点到表面点之间连线的夹角余弦，dS(r)是表面点r处的微元面积；

(B8)判断i＝N是否成立，如果成立则完成整个表面能量分布的重建，转向步骤C；否则i＝i+1并转到步骤B6，继续重建；

C.体内光源反演以及三维显示步骤：利用步骤B重建的生物体表面三维能量分布和步骤A获得的生物体解剖结构信息，反演生物体内的光源信息，并将反演的光源信息在生物体解剖结构中进行三维显示。

本发明与现有技术相比具有如下优点：

本发明由于将非接触式光学断层成像系统中的镜头简化为薄透镜模型，能够在大孔径镜头下进行多角度二维光学图像到任意形状生物体表面的三维能量分布重建；

本发明由于将大孔径镜头下生物体表面三维能量分布重建方法与体内光源反演算法相结合，可以实现大孔径镜头下的非接触式光学断层成像过程；

本发明由于采用薄透镜空间的大小限定透镜可见度因子的取值，通过在表面三维能量分布重建过程中调整薄透镜空间的大小，可以实现不同孔径镜头下的非接触式光学断层成像。

附图说明

图1是本发明的非接触式光学断层成像过程图；

图2是本发明非接触式光学断层成像过程中的表面三维能量分布重建过程图；

图3是采用非接触式光学断层成像系统的CCD相机采集的小鼠多角度光学图像；

图4是经过本发明对图3处理后的小鼠多角度光学图像；

图5是采用微计算机断层成像扫描获得的小鼠表面形状和内部解剖结构模型；

图6是本发明方法利用图4的小鼠多角度光学图像和图5的小鼠表面形状重建的小鼠表面三维能量分布图；

图7是本发明利用图6获得的小鼠表面三维能量分布进行小鼠体内光源重建的结果图。

具体实施方式

下面将结合附图详细描述本发明的重建方法，应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

参照图1，本发明的非接触式光学断层成像包括以下步骤：

步骤1，采集多角度光学图像和生物体表面形状及解剖结构数据。

(1.1)搭建多模态光学分子影像系统，该系统包括非接触式光学断层成像系统和微计算机断层成像系统两个子系统。其中非接触式光学断层成像系统由高性能的CCD相机和成像镜头构成，用于采集多角度的二维光学图像；微计算机断层成像系统由X-射线发射管和X-射线探测器构成，用于采集生物体的表面形状和解剖结构信息。这两个子系统在同一个光学平台上垂直正交放置，并共用一套控制设备，包括旋转台和平移台。在旋转台和平移台上放置一个生物体支架，用于固定成像小动物。

(1.2)用搭建的系统采集多角度光学图像和生物体表面形状及解剖结构数据：

首先，通过软件控制平移台调整生物体固定支架，使两个子系统采集的数据初始位置相同；

然后，将软件控制旋转台旋转3个90°，使非接触式光学断层成像系统的CCD相机能够采集到四个角度的二维光学图像；采集完四个角度的二维光学图像后，控制旋转台继续旋转90°，使生物体固定支架回到光学图像初始采集的位置；

最后，通过软件控制旋转台旋转小角度，一般为1°，微计算机断层成像系统采集360°的生物体投影数据，采集的投影数据经过计算机断层成像三维重建就可以得到生物体的三维表面形状信息和内部解剖结构信息。

步骤2，重建生物体表面三维能量分布。

参照图2，该重建步骤的具体实现如下：

(2.1)对四个角度的光学图像数据进行预处理。

在非接触式光学成像技术中，采用高灵敏度的光学检测仪器，如高性能CCD相机，采集并定量检测生物体内发射的光学信号，再将其转换成图像。进行光学数据采集时，一般采取制冷措施使CCD相机工作在黑暗的低温环境，以尽可能地减小背景噪声。但是当采集的光学信号较弱时，CCD相机的暗电流噪声和ADC读出噪声仍是不可忽视的背景噪声，如果任其存在将会对采集的光学图像的使用产生重大影响。由于现有制造工艺的缺陷，CCD相机的探测器上不可避免地存在一些坏点；这些坏点表现为与周围像素点的灰度值明显不符的亮点或暗点，如果不进行处理也将会影响光学图像的质量。此外，用CCD相机进行光学数据采集时，采集到的光学信号并不是分布在相机的整个探测器平面；由于CCD相机耦合成像镜头一起使用，因此在CCD相机的探测器平面上会存在一个采集物体的成像区域，这个成像区域称之为感兴趣区域。利用非接触式CCD相机采集的多角度光学图像进行生物体表面的三维能量重建，首先需要对光学图像中存在的背景噪声和坏点进行处理，同时还要对感兴趣区域进行标记，从而提高重建的效率。所以，多角度光学图像的预处理包括以下三个步骤：

①去除背景噪声。从信号处理的角度来看，背景噪声是直接叠加到测量的光学信号上的，所以可以通过背景相减去除，即：

I_P＝I_O-I_B

其中I_P是去除背景噪声之后的光学图像；I_O是CCD相机采集的原始光学图像；I_B是CCD相机采集光学图像I_O时的背景噪声，实验中取10次多次背景噪声的平均值；

②感兴趣区域提取。由于在非接触式光学断层成像系统中，CCD相机都是与镜头一起使用，因此感兴趣区域一般只占CCD相机探测器的一小部分。为了提高数据的利用效率，CCD相机采集的光学图像一般都是提取出感兴趣区域后再做下一步处理使用。经过第一步去除背景噪声处理后的光学图像可以通过下式进行感兴趣区域的提取：

I_ROI＝I_PΘW

其中I_ROI是进行感兴趣区域提取后的光学图像；Θ定义一种矩阵运算，运算结果为矩阵对应元素的乘积；W是用于提取感兴趣区域的窗函数矩阵，其任一元素w_st通过下式来确定：

$w_{\mathrm{st}}=\left\{\begin{array}{cc}1& (s,t)\∈{\mathrm{\Ω}}_{\mathrm{ROI}}\\ 0& (s,t)\∉{\mathrm{\Ω}}_{\mathrm{ROI}}\end{array}\right.$

其中Ω_ROI是成像体的成像区域；

③坏点补偿。通过邻域内的相邻点做插值运算来补偿光学图像中的坏点，本应用实例中采用四邻域相邻点的均值补偿坏点，但不局限于四邻域补偿。设D为光学图像数据I_ROI中的坏点集合，则坏点补偿公式为：

$I(s,t)=\left\{\begin{array}{cc}{I}_{\mathrm{ROI}}(s,t)& (s,t)\∉D\\ [I_{\mathrm{ROI}}(s-1,t)+I_{\mathrm{ROI}}(s,t-1)+I_{\mathrm{ROI}}(s,t+1)+I_{\mathrm{ROI}}(s+1,t)]/4& (s,t)\∈D\end{array}\right.$

其中I是经过背景噪声去除、感兴趣区域提取和坏点补偿预处理后的光学图像。

(2.2)预处理后的光学图像灰度值到能量值的转换。

步骤(2.1)虽然得到了经过预处理后的光学图像，但是此时的光学图像I仍然是图像的灰度表示，不能直接用于生物体表面三维能量的重建。考虑CCD相机的一些固有性质，通过下式实现未过曝光情况下光学图像灰度值到能量值的转换：

其中，I是光学图像的灰度值；P是经灰度到能量转换后的功率值；

是CCD相机的特性参数，定义为一定拍摄条件下，CCD相机的满井电子数与图像可能达到的最大灰度值的比值；Q是CCD相机的特性参数，定义为CCD相机的量子效率；t_p是光学图像数据拍摄时的曝光时间；E₀是单个光子携带的能量，通过下式计算：

$E_0=\frac{\mathrm{hc}}{\mathrm{\λ}}$

其中h是普朗克常量，其值为h＝6.626196×10^-34J·s；c是光在空气中的传播速度，其值为c＝3×10⁸m/s；λ是CCD相机采集的光学信号的中心波长。

(2.3)计算虚拟探测面的位置和像素尺寸。

在非接触式光学成像系统中，CCD相机耦合的镜头并不是由简单的一个透镜组成的，而是由一组透镜组合而成的。德国科学家Gauss已经证实了对于这样的透镜组合可以用由两个反射面组成的高斯厚透镜等效，这两个反射面称为高斯厚透镜模型的主平面。根据几何和光度学映射关系，高斯厚透镜可以进一步地简化为定位在第一个主平面的薄透镜模型。高斯厚透镜模型和简化后的薄透镜模型都遵循透镜的成像规律，从而满足如下的透镜成像公式：

$\frac{1}{u}+\frac{1}{v}=\frac{1}{f}$

其中，u是光学成像实验中的物距，v是光学成像实验中的像距，f是成像系统的焦距。但是，在实际的光学成像实验中，并不能直接地测量到物距和像距的数值，而是测量的物像距离l，其与物距和像距的关系表述为：

l＝u+v

因此，在已知成像系统焦距和物像距离的情况下，根据上面两式就可以计算成像实验中的物距和像距，同时也可以得到成像实验中的成像比t：

$t=\frac{v}{u}$

虚拟探测面定义为CCD相机相对于薄透镜模型的焦平面，则虚拟探测面上的像素点与CCD相机上的像素点之间满足如下对应关系：

$r_{\mathrm{vd}}=r_d+\frac{l}{\cos \mathrm{\θ}}{s}_{r_d\→r_{\mathrm{vd}}}$

其中，r_vd是虚拟探测面上的像素点坐标，r_d是CCD相机上的像素点坐标；cosθ是两个像素点之间的连线与主光轴的夹角余弦；

是CCD相机上像素点r_d到虚拟探测面上像素点r_vd的单位方向向量。根据上式就可以得到虚拟探测面和CCD相机的轴上中心点之间的关系：

r_vd(k)＝r_d(k)+l(n_d(k))

其中，k＝x，y，z代表坐标系的x，y，z三个轴向；n_d是CCD相机的单位法向，以指向成像体方向为正。例如，k＝x表示虚拟探测面和CCD相机位于坐标系的x轴上，并且|n_d(k)|＝1(若位于x正半轴，则n_d(k)＝-1；若位于x负半轴，则n_d(k)＝1)，上式描述的就是虚拟探测面和相机的中心x坐标之间的关系。在光学成像实验中，CCD相机的轴上中心点坐标可以通过物像距离和生物体三维表面轮廓得到。

根据真实实验中的成像比和CCD相机上每个像素的大小，计算虚拟探测面的像素大小：

${\mathrm{ps}}_{\mathrm{vd}}=\frac{{\mathrm{ps}}_d}{t}$

其中，ps_vd是虚拟探测面的像素大小，ps_d是CCD相机的像素大小。

(2.4)配准基于标记点的多角度光学图像到计算机断层成像数据重建的生物体结构信息坐标系

由于光学图像自身不存在绝对坐标系，并且光学断层成像重建时采用的生物体结构信息来源于计算机断层成像数据的重建，其采用的坐标系是计算机断层成像数据重建后的坐标系。因此，四个角度光学图像的中心点相对于计算机断层成像重建数据坐标系的实际物理坐标必须先确定，才能进行下一步的三维表面能量重建。

以x正半轴记为0°，旋转台逆时针旋转3个90°采集四个角度光学图像为例，基于标记点的多角度光学图像到计算机断层成像数据重建的生物体结构信息坐标系的配准可以通过如下步骤实现：

首先，读取光学图像和计算机断层成像重建后数据中的标记点，分别记录为：

0°： ${\mathrm{mark}}_0^{\mathrm{OT}}=[{\mathrm{mark}}_{0,y}^{\mathrm{OT}},{\mathrm{mark}}_{0,z}^{\mathrm{OT}}],$ ${\mathrm{mark}}_0^{\mathrm{CT}}=[{\mathrm{mark}}_{0,x}^{\mathrm{CT}},{\mathrm{mark}}_{0,y}^{\mathrm{CT}},{\mathrm{mark}}_{0,z}^{\mathrm{CT}}];$

90°： ${\mathrm{mark}}_{90}^{\mathrm{OT}}=[{\mathrm{mark}}_{90,x}^{\mathrm{OT}},{\mathrm{mark}}_{90,z}^{\mathrm{OT}}],$ ${\mathrm{mark}}_{90}^{\mathrm{CT}}=[{\mathrm{mark}}_{90,x}^{\mathrm{CT}},{\mathrm{mark}}_{90,y}^{\mathrm{CT}},{\mathrm{mark}}_{90,z}^{\mathrm{CT}}];$

180°： ${\mathrm{mark}}_{180}^{\mathrm{OT}}=[{\mathrm{mark}}_{180,y}^{\mathrm{OT}},{\mathrm{mark}}_{180,z}^{\mathrm{OT}}],$ ${\mathrm{mark}}_{180}^{\mathrm{CT}}=[{\mathrm{mark}}_{180,x}^{\mathrm{CT}},{\mathrm{mark}}_{180,y}^{\mathrm{CT}},{\mathrm{mark}}_{180,z}^{\mathrm{CT}}];$

270°： ${\mathrm{mark}}_{270}^{\mathrm{OT}}=[{\mathrm{mark}}_{270,x}^{\mathrm{OT}},{\mathrm{mark}}_{270,z}^{\mathrm{OT}}],$ ${\mathrm{mark}}_{270}^{\mathrm{CT}}=[{\mathrm{mark}}_{270,x}^{\mathrm{CT}},{\mathrm{mark}}_{270,y}^{\mathrm{CT}},{\mathrm{mark}}_{270,z}^{\mathrm{CT}}];$

其中，mark_k ^OT是光学图像的标记点，mark_k ^CT是计算机断层成像重建后数据中的标记点，k＝0，90，180，270表示四个不同的角度。

其次，根据计算机断层成像重建数据的分辨率，计算标记点在计算机断层成像重建数据坐标系中的实际物理坐标。计算机断层成像重建数据的分辨率记为 $R^{\mathrm{CT}}=[R_x^{\mathrm{CT}},R_y^{\mathrm{CT}},R_z^{\mathrm{CT}}],$ 则四个角度的标记点的实际物理坐标可由如下方式确定：

${\mathrm{mark}}_k={\mathrm{mark}}_k^{\mathrm{CT}}\mathrm{\Θ}{R}^{\mathrm{CT}}=[{\mathrm{mark}}_{k,x}^{\mathrm{CT}}\×R_x^{\mathrm{CT}},{\mathrm{mark}}_{k,y}^{\mathrm{CT}}\×R_y^{\mathrm{CT}},{\mathrm{mark}}_{k,z}^{\mathrm{CT}}\×R_z^{\mathrm{CT}}]$

其中，mark_k是k角度上标记点的实际物理坐标；Θ定义一种向量运算，运算结果为向量对应元素的乘积。

接着，根据标记点在光学图像数据中的位置以及各个角度光学图像数据中心点的位置，计算两者之间的像素距离dist_k，计算公式如下：

${\mathrm{dist}}_k={\mathrm{centr}}_k^{\mathrm{OT}}-{\mathrm{mark}}_k^{\mathrm{OT}}$

其中，centr_k ^OT是k角度上光学图像数据中心点的位置。

最后，依据光学图像数据相对于仿体的分辨率，计算四个角度光学图像数据中心点相对于计算机断层成像重建数据坐标系的实际物理坐标；四个角度光学图像数据中心点的实际物理坐标记为centr_k＝[centr_k，x，centr_k，y，centr_k，z]，其中k的意义同上，并通过如下过程确定中心点的实际物理坐标：

0°：x坐标centr_0，x通过步骤(2.3)得到，centr_0，y和centr_0，z通过下式计算

${\mathrm{centr}}_{0,y}={\mathrm{mark}}_{0,y}-{\mathrm{dist}}_{0,y}\frac{R_y^{\mathrm{OT}}}{t},$ ${\mathrm{centr}}_{0,z}={\mathrm{mark}}_{0,z}-{\mathrm{dist}}_{0,z}\frac{R_z^{\mathrm{OT}}}{t}$

90°：y坐标centr_90，y通过步骤(2.3)得到，centr_90，x和centr_90，z通过下式计算

${\mathrm{centr}}_{90,x}={\mathrm{mark}}_{90,x}-{\mathrm{dist}}_{90,x}\frac{R_x^{\mathrm{OT}}}{t},$ ${\mathrm{centr}}_{90,z}={\mathrm{mark}}_{90,z}-{\mathrm{dist}}_{90,z}\frac{R_z^{\mathrm{OT}}}{t}$

180°：x坐标centr_180，x通过步骤(2.3)得到，centr_180，y和centr_180，z通过下式计算

${\mathrm{centr}}_{180,y}={\mathrm{mark}}_{180,y}+{\mathrm{dist}}_{180,y}\frac{R_y^{\mathrm{OT}}}{t},$ ${\mathrm{centr}}_{180,z}={\mathrm{mark}}_{180,z}-{\mathrm{dist}}_{180,z}\frac{R_z^{\mathrm{OT}}}{t}$

270°：y坐标centr_270，y通过步骤(2.3)得到，centr_270，x和centr_270，z通过下式计算

${\mathrm{centr}}_{270,x}={\mathrm{mark}}_{270,x}+{\mathrm{dist}}_{270,x}\frac{R_x^{\mathrm{OT}}}{t},$ ${\mathrm{centr}}_{270,z}={\mathrm{mark}}_{270,z}-{\mathrm{dist}}_{270,z}\frac{R_z^{\mathrm{OT}}}{t}$

其中， $R^{\mathrm{OT}}=[R_x^{\mathrm{OT}},R_y^{\mathrm{OT}},R_z^{\mathrm{OT}}]$ 是光学图像的分辨率，t是光学成像实验中的成像比。这样，根据上述过程和步骤(2.3)能够得到最终的四个角度光学图像相对于计算机断层成像数据重建的生物体结构信息坐标系的中心物理坐标。

(2.5)建立CCD相机和虚拟探测面之间的能量映射关系。

对于理想的光学成像系统，像点和物点之间满足像的基本辐射亮度，即像点的辐射亮度等于物点的辐射亮度；所以，对于理想的光学断层成像系统，有如下表述：

L′(r_vd)＝L(r_d)

其中，L′(r_vd)是虚拟探测面上像素点的辐射亮度，L(r_d)是CCD相机上像素点的辐射亮度；r_vd是虚拟探测面上的点的坐标，r_d是CCD相机上的点的坐标，两者满足步骤203中描述的关系，即：

$r_{\mathrm{vd}}=r_d+\frac{l}{\cos \mathrm{\θ}}{s}_{r_d\→r_{\mathrm{vd}}}$

进一步地，可以将像点与物点辐射亮度之间的关系转换为辐射功率之间的关系，即：

P′(r_vd)＝P(r_d)

上式为不考虑能量损失的情况下，光学断层成像系统中虚拟探测面上点与CCD相机上点之间的功率对应关系。但是，实际的光学成像系统，能量损失总是存在的；如果不考虑光的传输损失，而只考虑成像系统带来的损失，则可以把损失归结为光学成像系统中薄透镜模型的透过率；所以，对于实际的光学断层成像系统，虚拟探测面上点与CCD相机上点之间的功率对应关系可以表述为：

$P^{\′}\left(r_{\mathrm{vd}}\right)=\frac{P\left(r_d\right)}{\mathrm{\τ}\left(f\right)}$

其中，τ(f)是光学断层成像系统中薄透镜模型的透过率。

(2.6)离散计算机断层成像重建的生物体表面，得到表面面片总数信息N、每个面片的表面法向信息以及三个顶点坐标信息；初始化当前取出面片编号i＝1。

(2.7)计算第i个面片的中心点r。

根据步骤(2.6)得到的每个面片的表面法向信息以及三个顶点坐标信息，结合几何学中三角形中心计算方法，计算当前取出面片的中心点r和表面法向n_r。

(2.8)计算当前取出面片与CCD相机上任意一个像素点之间的可见度因子。

根据步骤(2.7)已经得到当前取出面片的中心点r，通过下式计算其与CCD相机上任意一个像素点r_d之间的可见度因子：

ξ(r_d，r)＝α(r_vd，r)β(f，r)

其中，ξ(r_d，r)是从CCD相机看向生物体三维表面轮廓的可见度因子，其由虚拟可见度因子和透镜可见度因子两部分组成；虚拟可见度因子α(r_vd，r)是指从虚拟探测面看向生物体三维表面轮廓的可见度因子，其值与生物体三维表面轮廓、虚拟探测面的位置以及r_vd和r两个点的位置有关，可通过如下公式确定：

$\mathrm{\α}(r_{\mathrm{vd}},r)=\left\{\begin{array}{cc}1& (r_{\mathrm{vd}}\∈{\mathrm{\Ω}}^{\′})\∩(s_{r\→r_{\mathrm{vd}}}\∩S=\{r\left\}\right)\\ 0& \mathrm{else}\end{array}\right.$

其中，Ω′是虚拟探测面空间；S是生物体表面轮廓，r是生物体表面上的点并满足r∈S；

是从表面点r指向虚拟探测面上点r_vd的方向曲线；符号∩表示两者的交集。

透镜可见度因子β(f，r)是指从薄透镜模型看向生物体表面轮廓的可见度因子，其值除了与生物体表面轮廓、虚拟探测面的位置以及r_vd和r两个点的位置有关外，还与薄透镜模型的位置和大小有关，可通过如下公式确定：

其中，

与

意义相同，Ω_f是薄透镜模型所在空间。

r_d遍历四个角度CCD相机的所有像素点，就得到当前取出表面面片中心点r与四个角度CCD相机上所有像素点之间的可见度因子向量，从而建立了当前取出表面面片中心点与各个角度虚拟探测面上的所有像素点之间的空间位置对应关系。

(2.9)利用朗伯余弦定律和基本辐射亮度定理重建当前取出表面面片上的能量。

利用步骤(2.8)已经得到的当前取出表面面片中心点与四个角度CCD相机上所有像素点之间的可见度因子向量，通过下式计算当前取出表面面片上的能量重建值：

$P\left(r\right)=\frac{1}{\mathrm{\π}}{\∫}_{{\mathrm{\Ω}}^{\′}}{E}^{\′}\left(r_{\mathrm{vd}}\right)\mathrm{\ξ}(r_d,r)\frac{\cos {\mathrm{\θ}}_s\cos {\mathrm{\θ}}_{\mathrm{vd}}\mathrm{dS}\left(r\right)}{{|r-r_{\mathrm{vd}}|}^2}d{\mathrm{\Ω}}^{\′}$

其中，P(r)是当前取出表面面片中心点r处重建的功率；E′(r_vd)是虚拟探测面上点r_vd处的功率密度，其满足E′(r_vd)dΩ′＝P′(r_vd)；ξ(r_d，r)是表面点与CCD相机上像素点之间的可见度因子；cosθ_s是表面点法向与表面点和虚拟探测面上点之间连线的夹角余弦；cosθ_vd是虚拟探测面法向与表面点和虚拟探测面上点之间连线的夹角余弦；dS(r)是表面点r处的微元面积；dΩ′是虚拟探测面空间。结合虚拟探测面与CCD相机上各像素点之间的位置和能量对应关系，上述积分式也可以整理成如下的形式：

$P\left(r\right)=\frac{1}{\mathrm{\π}}{\∫}_{\mathrm{\Ω}}E\left(r_d\right)T(r_d,r)\mathrm{d\Ω}$

其中，E(r_d)是CCD相机上点r_d处的功率密度，其满足E(r_d)dΩ＝P(r_d)；dΩ是探测平面空间；T(r_d，r)定义为光在自由空间中的传输函数：

$T(r_d,r)=\frac{1}{\mathrm{\τ}\left(f\right)}\mathrm{\ξ}(r_d,r)\frac{1}{{|r-r_d-\frac{l}{\cos \mathrm{\θ}}{s}_{r_d\→r_{\mathrm{vd}}}|}^2}\cos {\mathrm{\θ}}_s\cos {\mathrm{\θ}}_d\mathrm{dS}\left(r\right)$

其中cosθ_d与cosθ_vd的意义和取值相同。

在四个角度的虚拟探测面空间或探测平面空间上求解上述的积分式，就可以得到当前取出表面面片上重建的能量值。

(2.10)判断表面三维能量重建是否完成。

判断当前取出的第i个表面面片是否为最后一个表面面片，即i＝N是否成立；如果当前取出表面面片是最后一个表面面片，则转到步骤3；否则设置i＝i+1，转向步骤(2.7)，重复上述的重建过程。

步骤3，体内光源反演以及三维显示步骤。

利用步骤2重建的生物体表面三维能量分布和步骤1获得的生物体解剖结构信息，反演生物体内的光源信息，并将反演的光源信息在生物体解剖结构中进行三维显示。

本发明的效果可以通过以下实验进一步说明：

在实验中，用一只29g的BALB/C小白鼠作为研究对象，分别进行非接触式光学断层成像和微计算机断层成像。

首先，把取自红色荧光棒(Glowproducts，Canada)的混合液体注入到小型胶囊作为小鼠体内的荧光光源。该混合溶液在充分混合后，会发出波长在650nm左右的荧光。

其次，在荧光光源植入到小鼠体内后，将小鼠用支架固定在旋转台上；通过旋转3个90°，使用科学级冷、背照式薄化CCD相机(PrincetonInstruments/Acton 2048B，Roper Scient-ific，Trenton，NJ)和光圈值f＝2.8的大孔径成像镜头采集小鼠前、右、后和左四个角度的光学信号，如图3所示，其中图3(a)是小鼠前视图，图3(b)是小鼠右视图，图3(c)是小鼠后视图，图3(d)是小鼠左视图。由于实验时荧光光源植入到靠近小鼠腹部的位置，因此穿透小鼠组织达到小鼠右侧和背部的荧光信号极其微弱，与小鼠腹部透射的信号相差很大，所以对小鼠四个角度的光学图像进行最大值归一化时，小鼠的右视图和后视图几乎看不到光信号，如图3(b)和(c)所示。同时，用橡皮泥球标记小鼠支架的四个角度平面，采集标记后小鼠的白光图像，用于与微计算机断层成像数据的配准。

最后，调整好初始角度，用与非接触式光学断层成像系统在同一水平台呈垂直角度的微计算机断层成像系统对小鼠进行扫描，得到全角度的计算机断层成像投影数据，用于小鼠体表和内部结构信息的重建。由于微计算机断层成像系统探测器大小的限制，实验中我们只采集了小鼠躯干部分的计算机断层成像数据，如图4所示，。

完成多角度光学图像和微计算机断层成像数据的采集后，根据本发明所描述的方法进行小鼠表面的三维能量分布重建和小鼠体内植入的荧光光源的重建。

首先，对CCD相机采集的四个角度的光学图像依次进行背景噪声去除、感兴趣区域提取和坏点补偿预处理操作以及灰度到能量的转换，其中灰度到能量的转换系数为9.037×10^-22，实验条件下，CCD相机的满井电子数与图像可能达到的最大灰度值的比值CCD相机在实验光波长下的量子效率约为0.94，采集荧光图像时的曝光时间为6min。经预处理和灰度到能量转换后的光学图像如图5所示，其中图5(a)是小鼠前视图，图5(b)是小鼠右视图，图5(c)是小鼠后视图，图5(d)是小鼠左视图。与图3采集的光学图像类似，由于实验时荧光光源植入到靠近小鼠腹部的位置，因此穿透小鼠组织达到小鼠右侧和背部的荧光信号极其微弱，与小鼠腹部透射的信号相差很大，所以对小鼠四个角度的光学图像进行最大值归一化时，小鼠的右视图和后视图几乎看不到光信号，如图5(b)和(c)所示。

其次，根据成像参数和标记点信息进行虚平面位置的确定以及光学图像到计算机断层成像重建数据的配准，配准后四个角度光学图像的中心点坐标如表1所示。

表1多角度光学图像到计算机断层成像重建数据的配准结果

视图	光学图像中心点坐标(mm)
		前视图	(37.0，24.5085，28.3083)
右视图	(24.3087，11.0，27.5908)
		后视图	(11.0，24.5086，28.6902)
左视图	(24.3086，38.0，28.8683)

接着，进行小鼠躯干的表面能量重建。重建时，首先利用Amira将小鼠体表离散成4000个面片，然后利用本发明中的方法重建这些面片上的能量分布，重建的结果如图6所示。

最后，利用图6所示重建的小鼠表面能量进行了小鼠体内光源的重建。光源重建时采用的重建算法是研究小组开发的自适应hp有限元光源重建算法，重建中采用的小鼠各个组织的光学参数如表2所示，重建光源在小鼠躯干部分的整体显示效果如图7所示。重建的光源与真实光源的中心距离误差为1.1976mm，其中真实光源的中心坐标由CT重建数据获得。此外，由于实验中的荧光光源没有进行标定，所以这里给出的表面重建结果和内部光源重建结果都只是定性的描述。

表2小鼠各个组织的光学参数，单位：mm^-1

	脂肪	心脏	肺	肝脏	肾
						μa	0.0086	0.1382	0.4596	0.8291	0.15498
μ′s	1.2584	1.0769	2.2651	0.7356	2.5329

从实验结果看，本发明的非接触式光学断层成像方法能够很好地解决大孔径镜头下任意形状生物体表面的三维能量分布重建，与生物体内部光源重建算法相结合，能够很好地实现大孔径镜头下的非接触式光学断层成像，重建光源与真实植入光源的距离误差小。如实验中采用的镜头光圈值为f＝2.8，采用自适应hp有限元光源重建算法重建的光源与真实植入光源的距离误差仅为1.1976mm。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (6)

1.一种非接触式光学断层成像方法，包括：

A.数据采集步骤：利用非接触式光学断层成像系统采集多角度的二维光学图像信息，同时利用微计算机断层成像或微核磁共振成像获得生物体的三维表面形状和内部解剖结构信息；

B.表面三维能量重建步骤：

(B1)对采集的多角度光学图像数据依次进行背景噪声去除、感兴趣区域提取和坏点补偿，以及将光学图像数据灰度值转换到对应的能量值；

(B2)将步骤A中使用的非接触式光学断层成像系统的成像镜头简化为薄透镜模型，并根据薄透镜成像规律计算CCD相机对应的虚拟探测面的位置和每个像素的实际尺寸；

(B3)将采集的二维光学图像信息与生物体表面三维形状信息进行配准，以获得各个角度二维光学图像数据中心点相对于生物体表面三维形状信息坐标系的实际物理坐标；

(B4)利用光学成像系统中像的基本辐射亮度定理，同时结合镜头的透过率，将步骤B1中获得的CCD相机上的能量分布映射到对应的虚拟探测面上；

(B5)离散步骤A中得到的生物体三维表面，记录其表面面片总数为N，N≥500，并记录各个表面面片的表面法向信息和三个顶点坐标信息；

(B6)取第i个表面面片的中心点，1≤i≤N，判断其与各个角度虚拟探测面上的所有像素点之间的空间位置对应关系，同时采用表格搜索策略计算其与CCD相机上任意一个像素点之间的可见度因子；

(B7)根据可见度因子的具体取值，应用下述公式计算第i个表面面片上的能量重建值：

$P\left(r\right)=\frac{1}{\mathrm{\π}}{\∫}_{\mathrm{\Ω}}E\left(r_d\right)T(r_d,r)\mathrm{d\Ω}$

式中，r是第i个表面面片的中心点；

P(r)是表面点r处重建的功率；

E(r_d)是CCD相机点r_d处的功率密度；

dΩ是CCD相机空间；

T(r_d，r)定义为光在自由空间的传输函数：

$T(r_d,r)=\frac{1}{\mathrm{\τ}\left(f\right)}\mathrm{\ξ}(r_d,r)\frac{1}{{|r-r_d-\frac{l}{\cos \mathrm{\θ}}{s}_{r_d\→r_{\mathrm{vd}}}|}^2}\cos {\mathrm{\θ}}_s\cos {\mathrm{\θ}}_d\mathrm{dS}\left(r\right)$

其中，τ(f)是非接触式光学断层成像系统中镜头的透过率，ξ(r_d，r)是表面点与CCD相机像素点之间的可见度因子，l是真实实验中的物像距离，cosθ是CCD相机上像素点r_d到虚拟探测面上像素点r_vd之间的连线与主光轴的夹角余弦，

是CCD相机上像素点r_d到虚拟探测面的像素

点r_vd的单位方向向量，cosθ_s是表面点的法向向量与表面点到虚拟探测面点之间连线的夹角余弦，cosθ_d是虚拟探测面点的法向向量与虚拟探测面点到表面点之间连线的夹角余弦，dS(r)是表面点r处的微元面积；

(B8)判断i＝N是否成立，如果成立则完成整个表面能量分布的重建，转向步骤C；否则i＝i+1并转到步骤B6，继续重建；

C.体内光源反演以及三维显示步骤：利用步骤B重建的生物体表面三维能量分布和步骤A获得的生物体解剖结构信息，反演生物体内的光源信息，并将反演的光源信息在生物体解剖结构中进行三维显示。

2.根据权利要求1所述的方法，其中步骤(B1)所述的将光学图像数据灰度值转换到对应的能量值，是在下述未过曝光的条件下进行：

式中，I是光学图像数据的灰度值；

P是光学图像数据经灰度到能量转换后的功率值；

是CCD相机的特性参数，定义为一定拍摄条件下，CCD相机的满井电子数与图像可能达到的最大灰度值的比值；

Q是CCD相机的特性参数，定义为CCD相机的量子效率；

t_p是光学图像数据采集时的曝光时间；

E₀是单个光子携带的能量。

3.根据权利要求1所述的方法，其中步骤(B2)所述的计算CCD相机对应的虚拟探测面的位置，是通过下述的CCD相机像素点与虚拟探测面像素点之间的位置对应关系计算：

$r_{\mathrm{vd}}=r_d+\frac{l}{\cos \mathrm{\θ}}{s}_{r_d\→r_{\mathrm{vd}}}$

式中，r_vd是虚拟探测面上的像素点坐标；

r_d是CCD相机上的像素点坐标。

4.根据权利要求1所述的方法，其中步骤(B2)所述的计算虚拟探测面的像素大小，是通过如下公式计算：

$p{s}_{\mathrm{vd}}=\frac{p{s}_d}{t}$

式中，ps_vd是虚拟探测面的像素大小；

ps_d是CCD相机的像素大小；

t是光学系统的成像比，通过真实实验中测得的像距与物距的比值计算获得。

5.根据权利要求1所述的方法，其中步骤(B3)所述的将采集的二维光学图像信息与生物体表面三维形状信息进行配准，按照如下步骤进行：

(5a)读取光学图像数据和三维生物体表面数据中的标记点，并记录标记点在各自坐标系中的位置；

(5b)根据三维生物体表面数据的分辨率，计算标记点在三维生物体表面数据坐标系中的实际物理坐标；

(5c)计算光学图像数据中标记点到中心点的距离；

(5d)根据光学图像数据的分辨率，计算四个角度光学图像数据中心点相对于三维生物体表面数据坐标系的实际物理坐标。

6.根据权利要求1所述的方法，其中步骤(B6)所述的采用表格搜索策略计算表面点与CCD相机上任意一个像素点之间的可见度因子，是通过如下公式计算：

ξ(r_d，r)＝α(r_vd，r)β(f，r)

式中，ξ(r_d，r)是从CCD相机看向生物体表面的可见度因子；

α(r_vd，r)是从虚拟探测面看向生物体表面的虚拟可见度因子，通过下述公式计算：

$\mathrm{\α}(r_{\mathrm{vd}},r)=\left\{\begin{array}{cc}1& (r_{\mathrm{vd}}\∈{\mathrm{\Ω}}^{\′})\∩(s_{r\→r_{\mathrm{vd}}}\∩S=\{r\left\}\right)\\ 0& \mathrm{else}\end{array}\right.$

其中，Ω′是虚拟探测面空间，S是生物体表面轮廓，r是生物体表

面上的点并满足r∈S，

是从表面点r指向虚拟探测面点r_vd的方向曲线，符号∩表示两者的交集；

β(f，r)是从薄透镜看向生物体表面的透镜可见度因子，通过下述公式计算：

其中，

与

意义相同，Ω_f是薄透镜空间。

Images