CN102488528A - 一种层析成像几何参数的校准方法 - Google Patents

Abstract

本发明提供了一种层析成像几何参数的校准方法，涉及层析成像技术领域。该方法包括以下步骤：采集层析成像中的全角度投影图像，对所述采集到的全角度投影图像分别进行负对数运算，将每个角度的负对数运算后的图像进行叠加，得到叠加图像，利用所述叠加图像建立与几何参数相关联的目标函数，利用单纯形-模拟退火算法对所述目标函数进行全局最小值优化。本发明不需要任何精密的模体就可以实现多个几何参数的标定；由于利用了叠加图像，因此对噪声不敏感，抗噪声能力强；此外，本发明采用的单纯形-模拟退火全局最优化算法，可以避免陷入局部极小值，保证了多个几何参数同时求解的精度与准确度。

Description

一种层析成像几何参数的校准方法

技术领域

本发明涉及层析成像技术领域，具体涉及一种层析成像几何参数的校准方法，尤其适合于锥形束X射线层析成像***的几何参数校准。

背景技术

锥形束X射线层析成像(cone beam computed tomography，CBCT)是临床医学、生物医学等领域非常重要的诊断与研究工具，它包括乳腺CT、牙科CT、四肢CT、小动物CT等。CBCT***的几何参数对获取高分辨率、低伪影的重建图像起着决定性作用，几何参数的些许偏差将会带来成像质量的严重下降。因此，需要对CBCT***的几何参数进行校准。

现有的CBCT***几何参数的校准方法可以分为两类：

一类是离线校准，即在对被测样品扫描之前，先对专用模体成像得到***的几何参数，然后保证***固定不动，再进行样品试验。这类方法中有些需要由几十个钢球构成的已知相对几何结构的精密模体(美国专利：US_2005_0117708_A1；美国专利：US_2006_0245628_A1；中国发明专利：ZL 200510045796.3；Physics In Medicine and Biology 54(2009)1633-1660；Physics In Medicine and Biology 54(2009)7239-7261；Medical Physics 37(2010)3844-3854；Medical Physics 38(2011)2829-2840；Medical Physics 31(2004)3242-3266；Medical Physics 32(2005)968-983)，有些需要精密加工的图形模具(美国专利：US_2007_0041508_A1；欧洲专利：WO_2010_050970_A1)、圆杆(中国发明专利：ZL 200910023137.8)、线模体(中国发明专利：ZL 200910087131.7)、已知相对垂直距离的两个滚珠模型(中国发明专利：ZL 200910188615.0；中国发明专利：ZL200910079277.7；中国发明专利：ZL 200610066252.X；Physics In Medicine andBiology 45(2000)3489-3508；Medical Physics 33(2006)1695-1706；IEEETransactions on Information Technology in Biomedicine 15(2011)655-660，)。该类方法对模体的要求较高，这些模体的加工与制作精度关乎着几何校准的准确度，实际操作中很容易引入各种机械或者测量误差。

另一类是在线校准，即利用被测样品自身的投影图计算出CBCT***的几何参数，不需要制作任何模体。这是因为在实际***中，由于机械结构受环境的影响会逐渐出现偏差，致使***的几何参数在动态的变化，特别在高分辨***中，微小的***几何参数的变化都会带来性能的急剧下降；此外，还有部分***需要频繁调整视场，放大倍率等，这使得离线校准的方法繁琐并且不可靠。中国发明专利ZL 200810136662.6通过对投影图进行二值化处理后求图像重心的方法得到了旋转中心位置；中国发明专利ZL200810112193.4利用正弦图找出物体边缘的投影点在探测器阵列中的位置，从而求出了旋转中心位置；文章Medical Physics 36(2009)48-58利用正弦图求解出了旋转轴位置和平面内的倾角；文章Physics In Medicine and Biology53(2008)3841-3861利用正弦图求得与几何参数相关联的代价函数，然后采用单纯形算法求得旋转轴位置及探测器的两个倾角，但该方法不适合大锥角的***，计算结果对感兴趣区域的选择较敏感，操作复杂，采用的单纯形算法易陷入局部极小从而得不到正确的结果；文章Medical Physics 38(2011)4934-4945是把几何参数代入到CT重建算法中，通过比较重建出的图像边缘的锐化程度得到了旋转轴位置、探测器的两个倾角及射线源到探测器的距离，该方法在求解过程中，是采用了穷举的方法，没有做优化，而且每一次调整几何参数都需要重建一次图像，过程复杂耗时。现有的在线校准方法的整体问题就是目标函数的建立非常困难、求解复杂，且精度较离线校准稍低。

发明内容

有鉴于此，本发明的目的在于提供一种层析成像几何参数的校准方法，用于在不需要精密模体的情况下，实现高精度校准。

本发明提供了一种层析成像几何参数的校准方法，包括以下步骤：

采集层析成像中的全角度投影图像，对所述采集到的全角度投影图像分别进行负对数运算，将每个角度的负对数运算后的图像进行叠加，得到叠加图像，利用所述叠加图像建立与几何参数相关联的目标函数，利用单纯形-模拟退火算法对所述目标函数进行全局最小值优化。

本发明利用了层析成像全角度投影图的叠加图像的对称性与几何参数的数学关系，采用了单纯形-模拟退火全局最优化算法精确求解出了图像重建所需要的关键几何参数，并利用求出的几何参数重建出高分辨率、低伪影的正确结果。本发明不需要任何精密的模体就可以实现多个几何参数的标定；由于利用了叠加图像，因此对噪声不敏感，抗噪声能力强；此外，本发明采用的单纯形-模拟退火全局最优化算法，可以避免陷入局部极小值，保证了多个几何参数同时求解的精度与准确度。

附图说明

图1为本发明实施例提供的层析成像几何参数校准***的结构图；

图2为本发明实施例提供的层析成像几何参数校准方法的流程图；

图3为本发明实施例提供的建立目标函数的流程图；

图4为本发明实施例中偏移探测器上像素点坐标与理想探测器平面上像素点坐标之间的对应关系示意图。

具体实施方式

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

本发明实施例提供的层析成像几何参数校准***的结构及待求几何参数如图1所示，射线源焦点1到探测器2之间的距离为d，射线源焦点1在探测器2上的投影为(u₀，v₀)，其中u₀为旋转轴位置，v₀为中平面位置，探测器2绕y轴旋转的角度记为η(即探测器的倾斜角skew)，探测器2绕z(v)轴旋转的角度记为θ(即探测器的滚转角roll)。射线源包括X射线，γ射线，可见光与近红外光。

本发明实施例提供的层析成像几何参数校准方法的流程如图2所示，该方法包括以下步骤：

步骤201、定义右手笛卡尔坐标系x-y-z。其中z轴是旋转台的旋转轴，y轴垂直于旋转轴且穿过射线源焦点，x轴垂直于y-z平面且穿过y轴与z轴的交点，探测器平面上横向定义为u轴，纵向定义为v轴。

步骤202、采集层析成像中的投影图像，即物体或者机架每旋转一个角度，探测器就采集一幅投影图像，直至等角度旋转完360度。

步骤203、对层析成像中采集到的全角度投影图像分别进行负对数运算。全角度是指投影图像覆盖了0～360度的旋转角度。运算方法由公式1给出：

$p\left(\mathrm{\φ}\right)=-\ln \left(\frac{I}{I_0\left(\mathrm{\φ}\right)}\right)---\left(1\right)$

其中，p(φ)是在投影角度为φ时的负对数运算结果，I是入射X射线强度，I₀(φ)是投影角度为φ的投影图上各点测量的出射X射线强度。

步骤204、将每个角度的负对数运算后的图像进行叠加，得到叠加图像。

$\mathrm{SOP}=\underset{\mathrm{\φ}=0}{\overset{2\mathrm{\π}}{\mathrm{\Σ}}}p\left(\mathrm{\φ}\right)---\left(2\right)$

其中，SOP为叠加图像。

步骤205、得到叠加图像后就开始建立与几何参数相关联的目标函数。目标函数的建立过程实际上就是对偏移的叠加图像绕不同的轴旋转并重新投影得到理想叠加图的过程。建立目标函数的流程如图3所示：

步骤2051、建立空间中点的坐标与探测器平面上像素点之间的对应关系。

如图4所示，P₀是没有任何几何偏移的理想的探测器采集到的叠加图像，P₁和P₂是探测器偏移后采集到的叠加图像。P₀，P₁和P₂三个叠加图像都有相同的矩阵大小M×N。S_0i，S_1i和S_2i分别表示叠加图像上的任意一个像素，它们在图像矩阵中的位置分别是(u_0m，v_0n)，(u_1m，v_1n)和(u_2m，v_2n)，它们在x-y-z坐标系中的坐标分别是(x_0i，y_0i，z_0i)，(x_1i，y_1i，z_1i)和(x_2i，y_2i，z_2i)，其中1≤m≤M，1≤n≤N，1≤i≤M×N。射线源焦点坐标为S(0，d，0)。

因为P₀和P₁都位于x-z平面内，所以P₀和P₁的平面方程都为y＝0，可知y_0i＝y_1i＝0。所以，探测器P₀上的每个像素(u_0m，v_0n)在x-y-z坐标系中的坐标位置的计算公式如下式所示：

$\left[\begin{array}{c}{x}_{0i}\\ z_{0i}\end{array}\right]=(\left[\begin{array}{c}{u}_{0m}\\ v_{0n}\end{array}\right]-\left[\begin{array}{c}{u}_0\\ v_0\end{array}\right])\×s_p---\left(3\right)$

其中，s_p是像素尺寸。

步骤2052、将叠加图像绕y轴旋转角度η，也就是把叠加图像的每一个像素对应的三维坐标与平面内的旋转矩阵相乘。

绕y轴旋转叠加图像P₀，旋转角度为η，得到叠加图像P₁，如图3(a)所示。P₁中的任意一个像素的坐标(x_1i，y_1i，z_1i)可以通过P₀中对应的像素坐标(x_0i，y_0i，z_0i)与绕y轴的旋转矩阵相乘得到，如下面的公式所示：

$\left[\begin{array}{c}{x}_{1i}\\ y_{1i}\\ z_{1i}\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}{c}{x}_{0i}\\ y_{0i}\\ z_{0i}\end{array}\right]=R_y\×\left[\begin{array}{c}{x}_{0i}\\ y_{0i}\\ z_{0i}\end{array}\right]---\left(4\right)$

其中R_y就是从P₀转换到P₁的旋转矩阵。

步骤2053、将叠加图像绕z轴旋转角度θ，即根据射线源到探测器的距离以及射线源焦点在探测器上的投影坐标对叠加图像进行CBCT的重新投影。

绕z(v)轴旋转叠加图像P₁，旋转角度为θ，得到叠加图像P₂，如图3(b)所示。可知P₂所在平面的平面方程为：

y_2i＝kx_2i，k＝tan(θ)(5)

如图3(b)所示，P₂上的任意一点的坐标S_2i就是所在的平面方程与通过点S和点S_1i的直线的交点，可以通过求解如下所示的方程得到点S_2i的坐标：

$\left\{\begin{array}{c}{y}_{2i}={\mathrm{kx}}_{2i}\\ \frac{x_{2i}-x_{1i}}{x_{2i}}=\frac{y_{2i}-y_{1i}}{y_{2i}-d}=\frac{z_{2i}-z_{1i}}{z_{2i}}\\ y_{1i}=0\end{array}\right.---\left(6\right)$

解方程求解出的点S_2i的坐标(x_2i，y_2i，z_2i)为：

$\left\{\begin{array}{c}{x}_{2i}=\frac{{\mathrm{dx}}_{1i}}{d+{\mathrm{kx}}_{1i}}\\ y_{2i}=\frac{{\mathrm{kdx}}_{1i}}{d+{\mathrm{kx}}_{1i}}\\ z_{2i}=\frac{{\mathrm{dz}}_{1i}}{d+{\mathrm{kx}}_{1i}}\end{array}\right.---\left(7\right)$

其中，点S_2i对应的P₂中的像素位置(u_2m，v_2n)可以通过下式得到：

$\left[\begin{array}{c}{u}_{2m}\\ v_{2n}\end{array}\right]=\left[\begin{array}{c}{x}_{2i}\\ z_{2i}\end{array}\right]\÷s_p+\left[\begin{array}{c}{u}_0\\ v_0\end{array}\right]---\left(8\right)$

至此，理想叠加图像P₀上的任意一点(u_0m，v_0n)的像素值均可以由偏移叠加图像上对应的点(u_2m，v_2n)及其周围的点通过双线性插值得到。

步骤2054、根据推导出的新的叠加图像，即绕Y轴和Z轴旋转后得到新的叠加图像，再计算新的叠加图像关于旋转轴对称的最小二乘误差，得到最终的以几何参数为变量的多参数目标函数。

P₀关于u₀对称的误差的计算公式如下：

$\begin{array}{cc}\mathrm{err}=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{u_0}{\mathrm{\Σ}}}{[h(m,n)-h(m,2u_0-n)]}^2& u_0\≤N/2\end{array}---\left(9\right)$

其中，err为误差值。

在理想条件下(即没有任何的几何参数误差)，叠加图像P₁和P₂与P₀是重合的，叠加图像P₀关于u₀完全对称，误差err等于0。但由于几何参数偏移的存在，误差err是大于0的。由上面的数学推导可知该误差是与***的几何参数紧密相关。计算误差err所代入的几何参数越接近真实的几何参数，那么对应的误差err就越小，因此求解真实几何参数的过程实际上就是求解误差err在全局最小值时对应的变量，它们的关系可由下面的数学式子表述，也即目标函数：

f(η，θ，u₀，v₀，d)＝minimize(err)    (10)

步骤206、利用单纯形-模拟退火算法对多参数的目标函数进行全局最小值优化，求解出探测器的倾斜角(skew)、滚转角(roll)、旋转轴位置、中平面位置以及射线源到探测器的距离这几个关键的几何参数，从而完成了层析成像***的几何校准。初始条件设置如下：

|u₀-u_in|≤40pixels

|v₀-v_in|≤40pixels

|η-η_in|≤4degrees    (11)

|d-d_in|≤100pixels

|θ-θ_in|≤4degrees

其中(u_in，v_in，η_in，d_in，θ_in)表示输入的理想的初始几何参数。

该算法结合了单纯形快速收敛与模拟退火可以跳出局部极小的优点，可以快速高精度的寻找到全局最优解。

利用本实施例，本申请人基于实验的CBCT做了大量的计算机仿真实验、水模实验、离体骨头实验和小动物活体实验，实验结果均验证了本实施例的有效性。本方法也可同样适用于工业CT中，包括千伏特(kV)和兆伏特(MV)的工业CT。此外，本方法还可以扩展到单光子发射层析成像(SPECT)和正电子发射层析成像(PET)中。

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

Claims (8)

1.一种层析成像几何参数的校准方法，其特征在于，包括以下步骤：

采集层析成像中的全角度投影图像，对所述采集到的全角度投影图像分别进行负对数运算，将每个角度的负对数运算后的图像进行叠加，得到叠加图像，利用所述叠加图像建立与几何参数相关联的目标函数，利用单纯形-模拟退火算法对所述目标函数进行全局最小值优化。

2.根据权利要求1所述的校准方法，其特征在于，该方法进一步包括以下步骤：

设定右手笛卡尔坐标系x-y-z，z轴为旋转台的旋转轴，y轴垂直于旋转轴且穿过射线源焦点，x轴垂直于y-z平面且穿过y轴与z轴的交点，探测器平面上横向设定为u轴，纵向设定为v轴。

3.根据权利要求1所述的校准方法，其特征在于，所述层析成像所使用的射线源具体包括：

X射线、γ射线、可见光或近红外光。

4.根据权利要求1所述的校准方法，其特征在于，所述全角度投影图像为该投影图像覆盖了0～360度的旋转角度。

5.根据权利要求1至4中任意一项所述的校准方法，其特征在于，所述建立目标函数的步骤具体包括：

建立空间中点的坐标与探测器平面上像素点之间的对应关系；

将所述叠加图像绕y轴旋转角度η；

再将所述叠加图像绕z轴旋转角度θ；

根据获得的绕y轴和z轴旋后新的叠加图像，计算所述新的叠加图像关于旋转轴对称的最小二乘误差，得到最终的以几何参数为变量的多参数目标函数。

6.根据权利要求5所述的校准方法，其特征在于，所述目标函数为：

f(η，θ，u₀，v₀，d)＝minimize(err)

其中，η，θ，u₀，v₀，d分别为探测器的倾斜角skew、探测器的滚转角roll、旋转轴位置、中平面位置以及射线源焦点到探测器的距离，err为误差值。

7.根据权利要求6所述的校准方法，其特征在于，所述对目标函数进行全局最小值优化的方法具体包括：

设置初始条件为：

|u₀-u_in|≤40pixels

|v₀-v_in|≤40pixels

|η-η_in|≤4degrees

|d-d_in|≤100pixels

|θ-θ_in|≤4degrees

其中，(u_in，v_in，η_in，d_in，θ_in)表示输入的理想的初始几何参数。

8.根据权利要求7所述的校准方法，其特征在于，该方法适用于千伏特kV和兆伏特MV的工业CT，单光子发射层析成像SPECT以及正电子发射层析成像PET。

Images