CN102499701B - X射线和荧光双模式活体成像系统的几何校准方法 - Google Patents

Abstract

本发明公开了一种X射线和荧光双模式活体成像系统的几何校准方法，包括：建立双模式活体成像系统的基本坐标系和相关坐标系；用小钢球标记经过视场的多条激发光束，得到激发光束的起点；改变双轴振镜的输入电压，使激发光束依次扫描视场内的多个位置，得到激发光束的方向参数；对视场内的多个磷光小球进行成像，建立空间点与光学投影点的关系；采集对象的X射线投影数据和荧光投影数据，并对X射线投影数据进行重建；生成用于重建荧光团分布的数据；进行荧光团分布重建，通过直接图像叠加完成图像配准。本发明实现更通用的几何校准方法，对双模式活体成像系统的系统结构限制更少，并且对象可以为任意复杂形状，实现来自于不同子系统的图像的配准。

Description

X射线和荧光双模式活体成像系统的几何校准方法

技术领域

本发明属于分子影像技术领域，涉及到成像系统的几何校准和图像处理，特别涉及到适用于复杂形状对象的非接触式X射线和荧光双模式活体成像系统的几何校准方法。

背景技术

传统的医学成像技术如X射线计算机断层成像(CT)、磁共振成像(MRI)和超声成像等主要利用生物体本身的物理特性或生理参数作为成像源。这些物理量或生理量对于与疾病或生理功能相关的细胞或分子没有特异性。而荧光分子断层成像技术能够对活体小动物体内的特异性荧光探针进行整体的三维、定量成像，具有非侵入、无电离辐射，低成本等优点。将X射线计算机断层成像与荧光分子断层成像技术相结合的X射线和荧光双模式活体成像系统，可以在同一实验平台上获得小动物的分子信息和结构信息，可能在疾病早期诊疗、药物研发和基础研究等方面发挥重要作用。近年来发展的非接触式探测技术利用电荷耦合器件(CCD)相机作为探测器，大大提高了光子的空间采样率；而适用于复杂形状对象的非接触式探测技术使成像过程中不需要挤压对象或利用光学参数匹配液简化边界条件，大大简化了实验操作。

发展适用于复杂形状对象的非接触式X射线和荧光双模式活体成像系统，必须获得描述各部件相互空间关系的几何参数集，并对来自于不同子系统的图像进行精确的配准。目前已有多种校准方法可实现对荧光分子断层成像子系统的探测器—CCD相机的校准，如美国专利US7949150Automaticcamera calibration and geo-registration using objects that provide positionalinformation，US6437823Method and system for calibrating digitalcameras，但都需要能够提供多个标记点位置信息的特殊校准模体，并且由于这些方法并非基于线性模型，因此无法采用光线跟踪方法判断探测器的有效性。

近年来，已有多种方法用来解决X射线和荧光双模式活体成像系统的几何校准及图像配准问题，但均对系统结构或对象的形状有较严格的限制。DaSilva等提出了一种几何和光学校准技术，该方法要求样品必须放在圆柱形容器内并浸泡在参数匹配液中，校准及实验过程非常复杂。Schulz等提出了一种荧光分子断层成像子系统的校准方法，但该方法只有在CCD相机的光轴平行于旋转架平面且旋转架平面垂直于旋转轴时才有效。Cao等提出了解析计算和最优化相结合的几何校准方法，并可实现子系统图像的直接融合。但该方法只允许CCD相机在一个方向上存在角度偏移，并且对校准模体的位置有严格的限定。

发明内容

有鉴于此，本发明的目的在于提供一种X射线和荧光双模式活体成像系统的几何校准方法，用于实现更通用的几何校准方法，对双模式活体成像系统的系统结构限制更少，并且对象可以为任意复杂形状，并发展相应的数据处理方法，实现来自于不同子系统的图像的配准。

本发明的实施例提供了一种X射线和荧光双模式活体成像系统的几何校准方法，包括：

建立双模式活体成像系统的基本坐标系和相关坐标系；

在所述基本坐标系与相关坐标系中，用小钢球标记经过视场的多条激发光束，得到激发光束的起点；

改变双轴振镜的输入电压，使激发光束依次扫描视场内的多个位置，得到激发光束的方向参数；

对视场内的多个磷光小球进行成像，建立空间点与光学投影点的关系；

采集对象的X射线投影数据和荧光投影数据，并对X射线投影数据进行重建；

利用得到的所述激发光束的起点、方向参数及空间点与光学投影点的关系，直接根据CT重建结果，生成用于重建荧光团分布的数据；

进行荧光团分布重建，通过直接图像叠加完成图像配准。

本发明建立了一种适用于复杂形状对象的非接触式X射线和荧光双模式活体成像系统的几何校准方法，得到荧光分子断层成像子系统的各个部件在双模式活体成像系统的基本坐标系内的位置和方向等几何参数，可以精确地描述双模式活体成像系统的结构；以CT重建结果为基础构建和提取重建荧光团分布的算法所需要的原始数据，因此子系统的重建结果在空间上是自然对应的，通过在三维空间的直接叠加即可实现图像配准，配准过程并无任何误差引入。本发明提供的方法适用于更广泛的系统结构和任意形状的对象，不仅可以获得精确的双模式活体成像系统结构的几何参数集，并且提供了一种包括几何校准，数据处理和图像配准的完整方法。

附图说明

图1为本发明实施例提供的适用于复杂形状对象的非接触式X射线和荧光双模式活体成像系统的几何校准方法的流程图；

图2为本发明实施例提供的双模式活体成像系统的基本坐标系和相关坐标系；

图3为本发明实施例提供的根据CT重建结果生成用于荧光团分布重建数据的方法流程图；

图4为基于本发明实施例提供的方法得到的组织模型切片上的探测器分布图和所有光源及探测器的分布图；

图5为本发明实施例中激发光图像和荧光图像中探测器的光学投影点的位置分布图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面结合附图对本发明作进一步的详细描述。

本发明实施例提供的X射线和荧光双模式活体成像系统的几何校准方法，首先建立了双模式活体成像系统的基本坐标系，并实现了对双模式活体成像系统的几何校准，将得到的几何参数集用于荧光分子断层成像原始数据的处理过程，最终得到的荧光团分布重建结果可通过直接图像叠加的方式实现与CT重建结果的配准。本发明实施例提供的方法适用于更广泛的系统结构和任意形状的对象，不仅可以获得精确的双模式活体成像系统结构的几何参数集，并且提供了一种包括几何校准、数据处理和图像配准的完整方法。

本发明实施例所述的适用于复杂形状对象的非接触式X射线和荧光双模式活体成像系统的几何校准方法，基于如下结构的双模式活体成像系统：CT子系统与荧光分子断层成像子系统具有共同的视场，载物台可带动对象进行360°旋转；荧光分子断层成像子系统利用双轴振镜实现激发光在样品表面的扫描，通过CCD相机采集从样品表面溢出的激发光信号和荧光信号。对激发光和探测器的校准是利用不同的模体独立进行的，相互之间并无任何关系。模体是对所用的小钢球和磷光小球的统称。

图1是本发明实施例提供的X射线和荧光双模式活体成像系统的几何校准方法的流程图，其实施步骤如下：

步骤101、建立X射线和荧光双模式活体成像系统的基本坐标系和相关坐标系。这里的相关坐标系是指以下三个坐标系：荧光分子断层成像子系统的坐标系O_f-X_fY_fZ_f、物理图像坐标系O₁-XY和像素坐标系O₂-UV。

采用针孔模型来描述双模式活体成像系统中的CCD相机，即CCD相机可以被模拟为一个成像面和投影中心的组合。为了描述对象的光学投影关系，建立如图2所示的四个相关坐标系：

在图2(a)中，O_c-X_cY_cZ_c为CT子系统的坐标系，将该坐标系视为双模式活体成像系统的基本坐标系，后面将要介绍的几何校准过程主要在该坐标系内实现。其中载物台的旋转轴被定义为Z_c轴，经过X射线管焦点且垂直于Z_c轴的轴被定义为X_c轴，垂直于X_cZ_c平面且经过X_c轴和Z_c轴交点的轴被定义为Y_c轴。O_f-X_fY_fZ_f为荧光分子断层成像子系统的坐标系，其原点O_f为投影中心，Z_f轴定义为沿着CCD相机光轴的方向，X_f轴和Y_f轴分别与成像面上像素行和列的方向平行。

在图2(b)中，图像坐标系O₁-XY和O₂-UV均定义在成像面上。其中O₁-XY为物理图像坐标系，其原点O₁定义在CCD相机的光轴与成像面的交点，该点通常位于图像中心处；X轴和Y轴分别定义为沿着成像面的水平和垂直方向，即像素的行和列方向。O₂-UV为像素坐标系，其原点O₂定义为图像最左下角的像素，U轴表示像素行方向，V轴表示像素列方向。

步骤102、用小钢球标记经过视场的多条激发光束，得到激发光束的起点。小钢球对X射线的吸收很强，在X射线投影图像中可以产生很好的对比度，在很多CT校准文章中都用小钢球作为样品。

首先为双轴振镜提供一对输入电压值，使激发光束经过视场，并在该光束的传播路径上固定两个半径为0.4mm的小钢球来标记该路径；然后改变输入电压使光束偏转，再次固定两个小钢球来标记偏转后的光束路径；如此循环操作，一共标记三条以上的激发光束。

然后用CT子系统对所有的钢球同时进行扫描及重建，并计算小钢球的坐标。只要确定了小钢球的坐标，即确定了所标记的光束路径。所有的激发光束均可视为由空间内的某一固定点发出，光束的方向随着双轴振镜的输入电压变化而改变，因此被标记的光束路径的交点即为所有激发光束的共同起点O_ex，利用最优化方法计算得到该起点的坐标。本实施例采用的最优化方法是Nelder-Mead单纯形法。

步骤103、改变双轴振镜的输入电压，使激发光束依次扫描视场内的多个位置，得到激发光束的方向参数。

激发光束方向D_ex的确定原理如下：如果用和θ分别表示激发光束方向与基本坐标系的X_c轴和Z_c轴正方向的夹角，则光束的单位方向向量可表示为：

其中，和θ₀分别是输入电压为(0,0)时光束方向与X_c轴和Z_c轴正方向的夹角；和Δθ分别为单位输入电压引起的X轴和Y轴振镜的偏转角；V_x和V_y分别为X轴振镜和Y轴振镜的输入电压。假设激发光束最终照射在视场内某一点P_ex，则光束的单位方向向量还可表示为：

$D_{\mathrm{ex}}=\frac{\stackrel{\→}{O_{\mathrm{ex}}{P}_{\mathrm{ex}}}}{|\stackrel{\→}{O_{\mathrm{ex}}{P}_{\mathrm{ex}}}|}---\left(2\right)$

在视场内任意多个位置放置小钢球，调整双轴振镜的输入电压使激发光束正好照射在小钢球上，记录下此时的电压值；然后利用CT子系统对所有的小钢球进行扫描和重建，得到小钢球的坐标；最后将记录的输入电压及小钢球的坐标作为已知量，代入公式(1)和(2)，利用最优化方法求得θ₀、和Δθ的值。

步骤104、对视场内的多个磷光小球进行成像，建立空间点与光学投影点的关系。

描述CCD相机的针孔模型需要以下几个关键参数：焦距f，像素的物理尺寸dx和dy，像主点的像素坐标(u₀,v₀)。但在大多数情况下，只关心投影点在投影图像上的像素坐标(u,v)，而不考虑投影点在空间中的位置。在这种情况下可以设焦距f＝1，将dx'＝dx/f，dy'＝dy/f称为归一化像素尺寸，该假设并不影响空间点和光学投影点的对应关系。

CT子系统的坐标系与荧光分子断层成像子系统的坐标系均为三维笛卡尔坐标系，因此它们之间的关系可以用旋转矩阵R与平移向量t来描述。旋转矩阵R取决于荧光分子断层成像子系统的坐标系相对于CT子系统的坐标系绕X_c、Y_c和Z_c轴的旋转角度α，β和γ，关系式如下：

$R=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\α}& -\sin \mathrm{\α}\\ 0& \sin \mathrm{\α}& \cos \mathrm{\α}\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\β}& 0& \sin \mathrm{\β}\\ 0& 1& 0\\ -\sin \mathrm{\β}& 0& \cos \mathrm{\β}\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\γ}& -\sin \mathrm{\γ}& 0\\ \sin \mathrm{\γ}& \cos \mathrm{\γ}& 0\\ 0& 0& 1\end{array}\right]---\left(3\right)$

平移向量t＝(t_x,t_y,t_z)^T，表示从CT子系统的坐标原点O_c到荧光分子断层成像子系统的坐标原点O_f的向量。因此空间点P在CT子系统的坐标系与荧光分子断层成像子系统的坐标系下的齐次坐标如果分别是(x_c,y_c,z_c,1)^T与(x_f,y_f,z_f,1)^T，则存在如下关系：

$\left[\begin{array}{c}{x}_f\\ y_f\\ z_f\\ 1\end{array}\right]=\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]=M_2\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]---\left(4\right)$

公式4中，0＝(0,0,0)^T，M₂为4×4矩阵。

根据图2中的几何关系，P点在荧光分子断层成像子系统的坐标系下的齐次坐标(x_f,y_f,z_f,1)^T与它的投影点P′的齐次物理图像坐标(x,y,1)^T具有如下关系：

$z_f\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{cccc}-1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{x}_f\\ y_f\\ z_f\\ 1\end{array}\right]---\left(5\right)$

而投影点P′的齐次物理图像坐标(x,y,1)^T与齐次像素坐标(u,v,1)^T有如下关系式：

$\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}-\frac{1}{{\mathrm{dx}}^{\text{'}}}& 0& u_0\\ 0& -\frac{1}{{\mathrm{dy}}^{\text{'}}}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]---\left(6\right)$

结合公式(4)，(5)和(6)，得到以下关系式：

$\begin{array}{c}{z}_f\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}-\frac{1}{{\mathrm{dx}}^{\text{'}}}& 0& u_0\\ 0& -\frac{1}{{\mathrm{dy}}^{\text{'}}}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}-1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]\\ =\left[\begin{array}{cccc}\frac{1}{{\mathrm{dx}}^{\text{'}}}& 0& u_0& 0\\ 0& \frac{1}{{\mathrm{dy}}^{\text{'}}}& v_0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]=M_1{M}_2\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]=M\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]\end{array}---\left(7\right)$

其中M为3×4矩阵，称为投影矩阵，它将空间点的坐标转化为投影点在图像上的像素坐标。M₁称为CCD相机的内部参数矩阵，M₂称为CCD相机的外部几何参数矩阵。

如果已知多个空间点的坐标及其光学投影点的像素坐标，即可利用最小二乘法计算投影矩阵M。CCD相机的几何参数根据以下公式求解：

r₃＝m₃₄m₃

$\begin{array}{c}{u}_0=(a_x{r}_1^T+u_0{r}_3^T)r_3=m_{34}^2{m}_1^T{m}_3,v_0=(a_y{r}_2^T+v_0{r}_3^T)r_3=m_{34}^2{m}_2^T{m}_3\\ {\mathrm{dx}}^{\text{'}}=\frac{1}{m_{34}^2|m_1\×m_3|},{\mathrm{dy}}^{\text{'}}=\frac{1}{m_{34}^2|m_2\×m_3|}\\ r_1=\frac{m_{34}}{a_x}(m_1-u_0{m}_3),r_2=\frac{m_{34}}{a_y}(m_2-v_0{m}_3)\\ t_z=m_{34},t_x=\frac{m_{34}}{a_x}(m_{14}-u_0),t_y=\frac{m_{34}}{a_y}(m_{24}-v_0)\end{array}---\left(8\right)$

其中，m_i ^T(i＝1～3)为投影矩阵M中第i行的前三个元素组成的行向量；m₃₄为投影矩阵M中第3行第4列的元素；r_i ^T(i＝1～3)为旋转矩阵R的第i行；t_x,t_y,t_z分别为平移向量t的三个分量；求得r₁、r₂、r₃后，即得到R矩阵。采用最优化方法，即可求得α，β和γ的值，至此就得到了CCD相机的全部几何参数值。

由以上叙述可知，CCD相机校准的关键之一是获得多个空间点的坐标及其光学投影点的像素坐标，这可以通过构建特殊的双模式校准模体——磷光小球来实现。该校准模体采用三个半径为0.4mm的小钢球，表面均匀镀一层磷光粉，制作成半径约0.5mm的磷光小球来模拟多个点光源，固定在视场内不同的位置。这些磷光小球在X射线投影图像和光学投影图像上都能产生很好的对比度。首先利用CT子系统进行扫描，旋转360°过程中采集400幅X射线投影图，然后利用荧光分子断层成像子系统从同样的初始角度开始进行光学扫描，旋转360°过程中采集200幅光学投影图。X射线投影图经CT算法重建后可计算磷光小球在初始角度的坐标，其他角度下的坐标可根据旋转角度推出。由光学投影图可得到磷光小球的投影点的像素坐标。由于磷光小球并不是理想的点光源，它在投影图像上的投影为一个覆盖超过一个像素的光斑，因此将投影图像经全局阈值处理并二值化后，再计算投影光斑质心的像素坐标。

步骤105、采集对象的X射线投影数据和荧光投影数据，并利用CT算法对X射线投影数据进行重建。对象是指要成像的物体，任何能发射荧光信号的小动物都可以作为对象，如体内有荧光团的小鼠等

X射线投影数据通常为400幅X射线投影图，载物台每旋转0.9°采集一幅；采集荧光投影数据的角度间隔通常为18～36°，每个角度下对数十至上百个激光光源位置采集光学投影图像，包括激发光图像，荧光图像，以及在相同实验条件但光源关闭情况下采集的背景噪声图像。

步骤106、利用得到的几何参数集，直接根据CT重建结果，生成用于重建荧光团分布的数据。几何参数集是指前面得到的所有几何参数，在步骤102-104中求出，包括：O_ex的三维坐标，θ₀、和Δθ的值，α，β和γ的值，(dx',dy')，(u₀,v₀)，平移向量t。

生成用于重建荧光团分布的数据的过程按图3所示的流程图进行：

步骤1061、计算在初始旋转角度下样品的边界体素坐标。

对每一幅CT切片进行二维中值滤波，滤除图像中的椒盐噪声，并保持图像边缘轮廓的细节。确定一个阈值T，用来区分CT切片中的空气和生物组织，利用该阈值将所有CT切片转化为二值化切片。为了减少后面介绍的光线求交测试的计算时间，可先进行像素合并后再进行二值化操作，或先进行二值化操作后再进行等间隔取样，形成一幅新的二值化切片，进行B×B的像素合并或沿行和列方向进行间隔为B个像素的等间隔取样，效果是类似的。假设原CT切片的大小为M×N，则进行本步骤的操作后，二值化切片的大小为U×V，其中U＝M/B,V＝N/B。对得到的每一幅二值化切片，利用Canny边缘检测器进行边缘检测，并对检测得到的边缘进行8连通测试，其中面积最大的边缘即为样品边界上的体素。8连通测试是图像处理中一种通用的测试像素连通性的方法，任何一个中心像素周围的8个相邻像素中，如果某一个像素与该中心像素具有相同的值，则认为该像素与中心像素连通。求出样品边界体素的坐标，其计算公式如下：

x_c＝(u-U/2-0.5)BS

y_c＝(v-V/2-0.5)BS (9)

z_c＝(-n_sli+N_sli/2+0.5)S

其中(u,v)为边界体素在切片上的行号和列号，n_sli为切片编号，S为像素尺寸。对每一幅CT切片进行相同的操作，求出所有在样品边界上的体素的坐标。该坐标值对应于初始旋转角度，即采集第一幅X射线投影图像时的角度，在数据采集过程中，随着载物台的旋转，边界体素的坐标随之变化。

步骤1062、计算光源的位置和方向。

采用光线跟踪法计算光源的位置和方向，即计算激发光束与样品的边界体素的交点坐标。根据对激发光的几何校准结果，结合实验过程中记录的双轴振镜的输入电压值，即可由公式(1)确定激发光束的起点和方向。在扫描过程中，样品绕Z_c轴旋转360°，因此对于某一边界体素，其z_c坐标值始终不变，改变的是x_c和y_c分量。在步骤1061中已经计算出在初始旋转角度下样品的边界体素坐标，可求得其对应的矢量的长度和方向：

$\begin{array}{c}\mathrm{\ρ}=\sqrt{{x_c}^2+{y_c}^2}\\ \mathrm{\θ}=\arctan \left(\frac{y_c}{x_c}\right)\end{array}---\left(10\right)$

在扫描过程中，该边界体素形成的轨迹即为该矢量在Z_c＝z_c平面顺时针旋转形成的圆，在转过Δθ角度后，该边界体素的坐标为：

x_c'＝ρcos(θ+Δθ)

y_c'＝ρsin(θ+Δθ)

z_c'＝z_c (11)

利用每一个边界体素的坐标构建一个小的立方体，其边长为1.5BS，即二值化切片中像素尺寸的1.5倍。计算得到激发光束与所有的边界体素形成的立方体的交点。如果交点不存在，则说明激发光没有照射在样品上，为无效光源；如果交点存在，取最近的交点作为光源位置，并将激发光束的方向作为光源的方向。对扫描过程中的所有角度进行同样的计算过程，得到所有光源的位置和方向。

步骤1063、建立组织模型。

首先确定组织模型的体素尺寸S_v。为了方便随后对CT切片的操作，令S_v为CT切片中像素尺寸S的整数倍，即S_v＝B₂S,B₂∈Z。然后确定重建区域，首先确定重建区域在Z_c方向对应的第一幅和最后一幅CT切片n_min和n_max，令n_max-n_min＝(N₂-1)B₂，这里N₂为组织模型沿Z_c方向的切片层数。从而，选定层号为{n_min,n_min+B₂,n_min+2B₂,...,n_max}的CT切片作为建立组织模型的依据。对于选定的每一层CT切片，先进行中值滤波，去除图像中的噪声；然后利用图像分割算法将切片分割为空气，软组织，骨骼等不同类型，用不同整数值表示不同组织类型。将图像分割后的切片进行等间隔采样，采样间隔为B₂，形成组织模型中的切片。

步骤1064、计算探测器的坐标和方向。

在非接触式光学成像中，CCD相机记录的是从样品边界出射的光强，因此可直接将探测器视为位于样品边界上，必须建立一种合理和简便地选择探测器的方法。

对组织模型中的一张切片，首先要确定切片上可作为探测器的所有体素。由于样品边界上所有体素都可以作为探测器，因此可利用前面所述的canny边缘检测器确定样品的边缘。但实际实验中发现该边缘检测方法并不太适用，这首先是因为在样品被区分为多种组织(如软组织和骨骼等)的情况下，样品内部不同组织类型的交界面同样会被认为是有效的边缘；其次，可作为探测器的边界体素可能为软组织，也可能为骨骼等其他组织类型，因此判断边缘的最佳阈值会随着组织类型而变化；再次，由于样品内部也可能存在空气，因此样品内部也存在空气与组织的交界面；最后，某些骨骼在切片中会从样品边界深入到内部，因此通过连通性测试排除面积较小的边缘的方法也不能完全排除所有内部的体素。因此，本实施例采用逐行或逐列判断的方法来确定可选探测器。

对切片中的每一行像素，根据每个像素的值判断其组织类型，如果全是空气，则该行像素均不能作为探测器，如果存在非空气像素，则取具有最小列号和最大列号的非空气像素作为可选探测器。对切片中的每一列像素，根据每个像素的值判断其组织类型，如果全是空气，则该列像素均不能作为探测器，如果存在非空气像素，则取具有最小行号和最大行号的非空气像素作为可选探测器。将沿行方向和沿列方向的可选探测器的并集作为该层切片上的可选探测器。可选探测器的数目通常大于实际需要的探测器数目，为了减少探测器的数目，需要对探测器进行等间隔的筛选，筛选原则如下：首先计算本层切片上所有可选探测器的坐标，计算公式为：

x_d＝(u_d-U/2-0.5)S_v

y_d＝(v_d-V/2-0.5)S_v

z_d＝(-n+N_sli+0.5)S (11)

其中u_d和v_d分别为可选探测器的行号和列号，并得到所有可选探测器坐标的平均值(重心)O_d(x_Od,y_Od)，作为该层上可选探测器的中心。然后计算从O_d到每个可选探测器的向量与X_c轴正方向的夹角，并按该夹角的大小对可选探测器进行排序，然后按设定的间隔B₃选择探测器。对n_min＜n＜n_max且间隔为B₃的所有切片进行相同操作，得到所有探测器的坐标和方向。

对大鼠头部成像的实例中，图4(a)为一张典型的组织模型切片，其中黑色表示空气，灰色表示软组织，白色表示骨骼等其他组织，用圆圈标示出的为选定的探测器；图4(b)为对感兴趣区域的光源和探测器分布，十字表示光源，圆圈表示探测器。

步骤1065、生成用于重建荧光团分布的原始数据文件。

根据采集投影图像时的旋转角度调整所有探测器的坐标。假设在采集图像时载物台的位置相对于初始位置旋转了Δθ，则其坐标的调整方法与边界体素坐标的调整方法相同，可利用公式(11)计算。根据公式(7)可由所有调整后的探测器坐标计算出探测器在投影图像上的像素坐标(u,v)。然后，利用光线跟踪方法，以像素(u,v)为起点，经过投影中心(其位置已经在CCD相机校准过程中求出)生成一条光线。该光线与样品边界必然有交点，最近的交点即为该像素对应的探测器。将所有探测器对应的像素的灰度值作为该探测器记录的光强值，并记录对应的光源编号，探测器编号。图5(a)和图5(b)为典型的大鼠头部激发光图像和荧光图像，圆点表示探测器在图像上的投影点。

步骤1066、进行荧光团分布重建，通过直接图像叠加完成图像配准。

荧光团分布重建得到荧光团在小动物体内的三维分布，是低分辨率的图像，通过插值，可以得到与重建区域内包含的CT重建结果相同大小的图像，可直接将相同位置的CT切片和荧光切片叠加。由于在本实施例提供的方法中，荧光分子断层成像子系统和CT子系统具有共同的基本坐标系，并且用于荧光团分布重建算法的体素模型是以计算机断层成像子系统的CT重建结果为依据建立的，因此本方法产生的的原始数据在空间上是配准的，通过直接叠加图像即可得到配准后的双模式图像，避免了传统的重建后配准引入的额外工作量和计算时间。

本实施例提出的双模式活体成像系统的几何校准方法中，对激发光和探测器的校准是利用不同的模体独立进行的，相互之间并无任何关系。因此本方法适用于任何以双轴振镜为扫描器件或以CCD相机为光学探测器的X射线和光学双模式成像系统，例如将CT与扩散光断层成像结合的双模式成像系统，将CT与生物发光断层成像结合的双模式成像系统等。

总之，以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。

Claims (12)

1.一种X射线和荧光双模式活体成像系统的几何校准方法，其特征在于，包括：

建立双模式活体成像系统的基本坐标系和相关坐标系；

在所述基本坐标系与相关坐标系中，用小钢球标记经过视场的多条激发光束，得到激发光束的起点；

改变双轴振镜的输入电压，使激发光束依次扫描视场内的多个位置，得到激发光束的方向参数；

对视场内的多个磷光小球进行成像，建立空间点与光学投影点的关系；

采集对象的X射线投影数据和荧光投影数据，并对X射线投影数据进行重建；

利用得到的所述激发光束的起点、方向参数及空间点与光学投影点的关系，直接根据CT重建结果，生成用于重建荧光团分布的数据；

进行荧光团分布重建，通过直接图像叠加完成图像配准。

2.根据权利要求1所述的几何校准方法，其特征在于，所述建立双模式活体成像系统的基本坐标系和相关坐标系包括：

建立CT子系统的坐标系O_c-X_cY_cZ_c，作为所述双模式活体成像系统的基本坐标系：原点为O_c，载物台的旋转轴被定义为Z_c轴，经过X射线管焦点且垂直于Z_c轴的轴被定义为X_c轴，垂直于X_cZ_c平面且经过X_c轴和Z_c轴交点的轴被定义为Y_c轴；

建立荧光分子断层成像子系统的坐标系O_f-X_fY_fZ_f，原点O_f为投影中心，Z_f轴定义为沿着CCD相机光轴的方向，X_f轴和Y_f轴分别与成像面上像素行和列的方向平行；

在成像面上建立图像坐标系O₁-XY和O₂-UV：其中O₁-XY为物理图像坐标系，原点O₁在CCD相机的光轴与成像面的交点，X轴和Y轴分别为沿着成像面的水平和垂直方向；O₂-UV为像素坐标系，原点O₂为图像最左下角的像素，U轴表示像素行方向，V轴表示像素列方向。

3.根据权利要求2所述的几何校准方法，其特征在于，所述得到激发光束的起点具体包括：

为双轴振镜提供一对输入电压值，使激发光束经过视场，并在该光束的传播路径上固定两个所述小钢球来标记该路径；

改变输入电压使光束偏转，再次固定两个所述小钢球来标记偏转后的光束路径；

循环执行以上步骤，共标记三条以上的激发光束；

用CT子系统对所有的钢球同时进行扫描及重建，并计算小钢球的坐标；

被标记的光束路径的交点为所有激发光束的共同起点O_ex，利用最优化方法计算得到该起点的坐标。

4.根据权利要求3所述的几何校准方法，其特征在于，所述得到激发光束的方向参数具体包括：

激发光束的单位方向向量为：

其中，和θ分别表示激发光束方向与基本坐标系的X_c轴和Z_c轴正方向的夹角，和θ₀分别是输入电压为(0,0)时光束方向与X_c轴和Z_c轴正方向的夹角，和Δθ分别为单位输入电压引起的X轴和Y轴振镜的偏转角，V_x和V_y分别为X轴振镜和Y轴振镜的输入电压；

假设激发光束最终照射在视场内一点P_ex，则光束的单位方向向量为：

$D_{\mathrm{ex}}=\frac{\stackrel{\→}{O_{\mathrm{ex}}{P}_{\mathrm{ex}}}}{|\stackrel{\→}{O_{\mathrm{ex}}{P}_{\mathrm{ex}}}|}---\left(2\right)$

在视场内任意多个位置放置小钢球，调整双轴振镜的输入电压使激发光束正好照射在小钢球上，记录下此时的电压值；然后利用CT子系统对所有的小钢球进行扫描和重建，得到小钢球的坐标；最后将记录的输入电压及小钢球的坐标作为已知量，代入公式(1)和(2)，利用最优化方法求得θ₀、和Δθ的值。

5.根据权利要求4所述的几何校准方法，其特征在于，所述建立空间点与光学投影点的关系具体包括：

描述CCD相机针孔模型的参数为：归一化像素尺寸(dx',dy')和像主点的像素坐标(u₀,v₀)，其中dx'＝dx/f，dy'＝dy/f，(dx,dy)为实际像素尺寸，f为焦距；

CT子系统的坐标系与荧光分子断层成像子系统的坐标系之间的关系用旋转矩阵R与平移向量t来描述；旋转矩阵R取决于荧光分子断层成像子系统的坐标系相对于CT子系统的坐标系绕X_c、Y_c和Z_c轴的旋转角度α，β和γ，关系式如下：

$R=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\α}& -\sin \mathrm{\α}\\ 0& \sin \mathrm{\α}& \cos \mathrm{\α}\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\β}& 0& \sin \mathrm{\β}\\ 0& 1& 0\\ -\sin \mathrm{\β}& 0& \cos \mathrm{\β}\end{array}\right]\left[\begin{array}{ccc}\cos \mathrm{\γ}& -\sin \mathrm{\γ}& 0\\ \sin \mathrm{\γ}& \cos \mathrm{\γ}& 0\\ 0& 0& 1\end{array}\right]---\left(3\right)$

平移向量t＝(t_x,t_y,t_z)^T，表示从CT子系统的坐标原点O_c到荧光分子断层成像子系统的坐标原点O_f的向量；空间点P在CT子系统的坐标系与荧光分子断层成像子系统的坐标系下的齐次坐标如果分别是(x_c,y_c,z_c,1)^T与(x_f,y_f,z_f,1)^T，则存在如下关系：

$\left[\begin{array}{c}{x}_f\\ y_f\\ z_f\\ 1\end{array}\right]=\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]=M_2\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]---\left(4\right)$

其中，0＝(0,0,0)^T，M₂为4×4矩阵；

P点在荧光分子断层成像子系统的坐标系下的齐次坐标(x_f,y_f,z_f,1)^T与它的投影点P′的齐次物理图像坐标(x,y,1)^T具有如下关系：

$z_f\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]=\left[\begin{array}{cccc}-1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{x}_f\\ y_f\\ z_f\\ 1\end{array}\right]---\left(5\right)$

而投影点P′的齐次物理图像坐标(x,y,1)^T与齐次像素坐标(u,v,1)^T有如下关系：

$\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}-\frac{1}{{\mathrm{dx}}^{\text{'}}}& 0& u_0\\ 0& -\frac{1}{{\mathrm{dy}}^{\text{'}}}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}x\\ y\\ 1\end{array}\right]---\left(6\right)$

结合公式(4)，(5)和(6)，得到以下关系：

$\begin{array}{c}{z}_f\left[\begin{array}{c}u\\ v\\ 1\end{array}\right]=\left[\begin{array}{ccc}-\frac{1}{{\mathrm{dx}}^{\text{'}}}& 0& u_0\\ 0& -\frac{1}{{\mathrm{dy}}^{\text{'}}}& v_0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}-1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]\\ =\left[\begin{array}{cccc}\frac{1}{{\mathrm{dx}}^{\text{'}}}& 0& u_0& 0\\ 0& \frac{1}{{\mathrm{dy}}^{\text{'}}}& v_0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]=M_1{M}_2\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]=M\left[\begin{array}{c}{x}_c\\ y_c\\ z_c\\ 1\end{array}\right]\end{array}---\left(7\right)$

其中M为3×4投影矩阵，M₁为CCD相机的内部参数矩阵，M₂为CCD相机的外部几何参数矩阵；

如果已知多个空间点的坐标及其光学投影点的像素坐标，利用最小二乘法计算投影矩阵M；CCD相机的几何参数根据以下公式求解：

r₃＝m₃₄m₃

$\begin{array}{c}{u}_0=(a_x{r}_1^T+u_0{r}_3^T)r_3=m_{34}^2{m}_1^T{m}_3,v_0=(a_y{r}_2^T+v_0{r}_3^T)r_3=m_{34}^2{m}_2^T{m}_3\\ {\mathrm{dx}}^{\text{'}}=\frac{1}{m_{34}^2|m_1\×m_3|},{\mathrm{dy}}^{\text{'}}=\frac{1}{m_{34}^2|m_2\×m_3|}\\ r_1=\frac{m_{34}}{a_x}(m_1-u_0{m}_3),r_2=\frac{m_{34}}{a_y}(m_2-v_0{m}_3)\\ t_z=m_{34},t_x=\frac{m_{34}}{a_x}(m_{14}-u_0),t_y=\frac{m_{34}}{a_y}(m_{24}-v_0)\end{array}---\left(8\right)$

其中，m_i ^T，i＝1～3为投影矩阵M中第i行的前三个元素组成的行向量；m₃₄为投影矩阵M中第3行第4列的元素；r_i ^T，i＝1～3为旋转矩阵R的第i行；t_x,t_y,t_z分别为平移向量t的三个分量；求得r₁,r₂,r₃后，得到R矩阵；采用最优化方法，能够求得α，β和γ的值。

6.根据权利要求5所述的几何校准方法，其特征在于，所述生成用于重建荧光团分布的数据具体包括：

计算在初始旋转角度下样品的边界体素坐标；

计算光源的位置和方向；

建立组织模型；

计算探测器的坐标和方向；

生成用于重建荧光团分布的原始数据文件。

7.根据权利要求6所述的几何校准方法，其特征在于，所述计算在初始旋转角度下样品的边界体素坐标具体包括：

对每一幅CT切片进行二维中值滤波，滤除图像中的椒盐噪声，并保持图像边缘轮廓的细节；

确定一个阈值T，用来区分CT切片中的空气和生物组织，利用该阈值将所有CT切片转化为二值化切片；

进行像素合并后再进行二值化操作，或先进行二值化操作后再进行等间隔取样，形成一幅新的二值化切片，进行B×B的像素合并或沿行和列方向进行间隔为B个像素的等间隔取样；

假设原CT切片的大小为M×N，则进行本步骤的操作后，二值化切片的大小为U×V，其中U＝M/B,V＝N/B；

对得到的每一幅二值化切片，利用Canny边缘检测器进行边缘检测，并对检测得到的边缘进行8连通测试，其中面积最大的边缘为样品边界上的体素；求出样品边界体素的坐标，计算公式如下：

x_c＝(u-U/2-0.5)BS

y_c＝(v-V/2-0.5)BS (9)

z_c＝(-n_sli+N_sli/2+0.5)S

其中(u,v)为边界体素在切片上的行号和列号，n_sli为切片编号，S为像素尺寸；对每一幅CT切片进行相同的操作，求出所有在样品边界上的体素的坐标；该坐标值对应于初始旋转角度，在数据采集过程中，随着载物台的旋转，边界体素的坐标随之变化。

8.根据权利要求7所述的几何校准方法，其特征在于，所述计算光源的位置和方向具体包括：

采用光线跟踪法计算光源的位置和方向，计算激发光束与样品的边界体素的交点坐标；根据对激发光的几何校准结果，结合双轴振镜的输入电压值，由公式(1)确定激发光束的起点和方向；根据初始旋转角度下样品的边界体素坐标，求得其对应的矢量的长度和方向：

$\begin{array}{c}\mathrm{\ρ}=\sqrt{{x_c}^2+{y_c}^2}\\ \mathrm{\θ}=\arctan \left(\frac{y_c}{x_c}\right)\end{array}---\left(10\right)$

在扫描过程中，该边界体素形成的轨迹为该矢量在Z_c＝z_c平面顺时针旋转形成的圆，在转过Δθ角度后，该边界体素的坐标为：

x_c'＝ρcos(θ+Δθ)

y_c'＝ρsin(θ+Δθ)

z_c'＝z_c (11)

利用每一个边界体素的坐标构建一个小的立方体，边长为1.5BS；计算得到激发光束与所有的边界体素形成的立方体的交点；如果交点不存在，则说明激发光没有照射在样品上，为无效光源；如果交点存在，取最近的交点作为光源位置，并将激发光束的方向作为光源的方向；对扫描过程中的所有角度进行同样的计算过程，得到所有光源的位置和方向。

9.根据权利要求8所述的几何校准方法，其特征在于，所述建立组织模型具体包括：

首先确定组织模型的体素尺寸S_v；令S_v为CT切片中像素尺寸S的整数倍，S_v＝B₂S,B₂∈Z；

然后确定重建区域：首先确定重建区域在Z_c方向对应的第一幅和最后一幅CT切片n_min和n_max，令n_max-n_min＝(N₂-1)B₂，N₂为组织模型沿Z_c方向的切片层数；选定层号为{n_min,n_min+B₂,n_min+2B₂,...,n_max}的CT切片作为建立组织模型的依据；

对于选定的每一层CT切片，先进行中值滤波，去除图像中的噪声；然后利用图像分割算法将切片分割为不同类型，用不同整数值表示不同组织类型；所述不同类型包括空气，软组织，骨骼；将图像分割后的切片进行等间隔采样，采样间隔为B₂，形成组织模型中的切片。

10.根据权利要求9所述的几何校准方法，其特征在于，所述计算探测器的坐标和方向具体包括：

采用逐行或逐列判断的方法来确定可选探测器；

对切片中的每一行像素，根据每个像素的值判断其组织类型，如果全是空气，则该行像素均不能作为探测器，如果存在非空气像素，则取具有最小列号和最大列号的非空气像素作为可选探测器；

对切片中的每一列像素，根据每个像素的值判断其组织类型，如果全是空气，则该列像素均不能作为探测器，如果存在非空气像素，则取具有最小行号和最大行号的非空气像素作为可选探测器；

将沿行方向和沿列方向的可选探测器的并集作为该层切片上的可选探测器；

对探测器进行等间隔的筛选，筛选方法为：首先计算本层切片上所有可选探测器的坐标，计算公式为：

x_d＝(u_d-U/2-0.5)S_v

y_d＝(v_d-V/2-0.5)S_v

z_d＝(-n+N_sli+0.5)S (12)

其中u_d和v_d分别为可选探测器的行号和列号，并得到所有可选探测器坐标的平均值O_d(x_Od,y_Od)，作为该层上可选探测器的中心；然后计算从O_d到每个可选探测器的向量与Xc轴正方向的夹角，并按该夹角的大小对可选探测器进行排序，然后按设定的间隔B₃选择探测器；对n_min<n<n_max且间隔为B₃的所有切片进行相同操作，得到所有探测器的坐标和方向。

11.根据权利要求9所述的几何校准方法，其特征在于，所述生成用于重建荧光团分布的原始数据文件具体包括：

根据采集投影图像时的旋转角度调整所有探测器的坐标；假设在采集图像时载物台的位置相对于初始位置旋转了Δθ，则其坐标的调整方法与边界体素坐标的调整方法相同，利用公式(11)计算；根据公式(7)由所有调整后的探测器坐标计算出探测器在投影图像上的像素坐标(u,v)；然后，利用光线跟踪方法，以像素(u,v)为起点，经过投影中心生成一条光线；该光线与样品边界最近的交点为该像素对应的探测器；将所有探测器对应的像素的灰度值作为该探测器记录的光强值，并记录对应的光源编号，探测器编号。

12.根据权利要求1至11中任意一项所述的几何校准方法，其特征在于，在所述几何校准方法中，对激发光和探测器的校准是利用不同的模体独立进行的，相互之间并无任何关系。