3-D~2-D registration of CT and MR to X-ray images

Tomazevic D, Likar B, Slivnik T, et al. 3-D/2-D registration of CT and MR to X-ray images[J]. IEEE transactions on medical imaging, 2003, 22(11): 1407-1416.

摘要

图像引导治疗的一个关键部分是术前和术中图像配准，通过配准可以在三维空间中确定患者解剖结构的精确位置和方向。
本文提出了一种新的方法来配准三维图像（CT或MR）到一个或者多个二维图像（X射线）。
配准仅基于二维和三维图像中存在的信息，不需要基准标记、术中X射线图像分割或者数字重建影像的适时构建。
该方法的原创性在于使用术前在三维图像（CT或MR）中定义骨骼的表面法线，使用术中由X射线源与三维表面点的定义的位置上的X射线图像梯度。
配准涉及寻找CT容积或者MR容积的刚性变换，考虑到它们的振幅和方向，提供了表面法线和反投影梯度之间的最佳匹配。
本文通过使用尸体腰椎模型的MR、CT和X射线图像，彻底验证了本文提出的配准方法，通过基准标记建立了尸体腰椎模型的“金标准”配准，并且通过目标配准误差（Target Registration Error，TRE）评估了其准确性。
包含了单个椎骨L1到L5的感兴趣的容积，从不同的起始位置被配准到不同的X射线对，这些位置是围绕“金标准”位置随机地和均匀地选择的。
从平移小于6mm（旋转小于$17^\circ$）或者3mm（$8.6^\circ$）的“金标准”开始，在超过$91%$的试验中快速的CT/X射线或MR/X射线配准都取得了成功（除L1为$82%$）。CT与X射线配准的均方根目标配准误差低于0.5mm，MR与X射线配准的均方根目标配准误差低于1.4mm。

索引

三维/二维图像配准，计算机断层扫描，“金标准”，图像引导治疗，磁共振，脊柱模型，X射线投影

Ch01 介绍

医学成像目前正经历快速发展，重点放在与传统方法的对比，使用成像技术使得外科手术和治疗程序的侵入性越来越小，并且提高特定程序的精确度[^1][^2]。
在影像引导治疗（image-guided therapy, IGT）（手术、放射治疗、放射介入）中，通常三维计算的断层摄影或者磁共振图像作为术前医疗数据用于诊断、计划、模拟、引导或者其他方式辅助外科医生，也可能是辅助机器人执行外科手术或者治疗过程。
该计划是在与术前数据相关的坐标系中构建的，而手术过程是在与病人相关的坐标系中执行的。
术前数据和计划与患者在治疗期间占据的物理空间之间的关系或者空间转换通过配准来建立。
配准是IGT的关键部分，它允许术前图像中定义的任何三维点在患者坐标系中的精确定位，并且因此可以向外科医生提供关于器械与计划轨迹、附近易受损结构或者最终目标的位置信息的相关性。
然而，为了适应临床应用，配准算法必须关注以下几点：配准精度、计算时间、对异常值（例如：由手术器械产生的异常值）的鲁棒性，以及术中数据采集的复杂性和侵入性。

在过去的十年中，已经提出了多种方法用于图像到物理空间（患者）的配准。这些方法的不同之处在于配准所依赖的信息。
最精确的方法是利用术前固定在患者解剖结构上的基准标记点，最常见的是骨骼上[^3][^4][^5]。
使用定位装置在手术中测量标记点的位置，这些标记点在三维CT或者MR图像中很容易检测到。
基于点的配准是通过寻找使得两个空间中的基准点对齐的刚性变换来实现的。
尽管这些方法快速、准确且稳健，但是它们需要将基准标记固定到刚性结构上，这并不总是可能的，或者侵入性太大而不可接受。
基于基准标记的方法的优势在于，配准误差仅取决于基准定位，因此在很大程度上与被配准的特定结构无关[^6]。
因为在这两个空间中基准标记定位可能非常精确，所以通过这种方法进行的配准通常被用作参考或者“金标准”，适用于评估其他配准方法[^7][^8][^9]。

第二类配准方法是基于表面[^9][^10][^11][^12][^13][^14]。这些方法依赖于对解剖结构形状的描述，通常是从手术前的三维图像中提取的皮肤或者骨骼的外表面。
术中模式，用于获取相同结构的形状的相似表示，包括数字化探头、被动式立体视觉和结构光。
然后使用最小化点到表面或者表面到表面距离的配准过程来导出空间变换[^15]。
这些技术的问题是皮肤不是刚性结构，并且包含足够形状信息的骨骼表面必须在手术过程中暴露。

第三类配准方法依赖于校准的X射线设备采集的术中X射线投影。
相对于X射线设备的几何性，在三维CT或MR图像中结构的方向与位置由三维/二维配准来确定。
当前存在两种不同的三维/二维配准方法。
第一种方法是基于特征的方法，需要从术前数据中提取骨骼表面，并且从一个或者多个术中投影图像中提取相同结构的轮廓[^7][^16]。
通过最小化表面模型和连接轮廓点与X射线源的线之间的距离来执行配准。
一旦进行了分割，减少的数据量使得配准更快。不幸的是，术中分割很难自动实现，并且分割中的错误会导致配准中的错误。
第二种方法分别基于三维术前和二维术中图像的体素和像素强度信息。
从CT图像产生模拟的X射线投影图像[^17]，称为数字重建放射图（Digitally Reconstructed Radiograph, DDR），并且通过优化从DDR和X射线图像[^8][^18]计算的相似性度量来估计CT体积相对于X射线图像集的未知姿态。
这些算法需要很少或者不需要分割，但是非常耗时。
通过计算仅包含感兴趣结构的DDR[^19]或者应用更快的DDR生成方法[^20][^21]，可以在一定程度上加快速度。
然而，这些方法不能适用于术前MR图像与X射线图像配准，因为实践中在基于MR的DDR和X射线投影之间没有相关性。

在本文中，我们提出了一种新的方法，用于将三维CT或者MR图像配准到二维X射线图像，最终目标是在手术或者治疗期间估计患者解剖结构的位置和方向。
所提出的方法结合了两种三维/二维方法的优点，既像第一种方法但是所需要的数据量更少，配准速度更快，并且不需要第二种方法的术中分割。
像第一种方法一样，该方法假设在手术前从CT或者MR图像中提取骨结构的表面。该方法使用一种新的准则函数（criterion function, CF）来测量表面的法线和X射线图像的术中分割，从而消除了由术中分割误差引起的配准误差。
我们使用刚性配准，当成像对象本身是刚性的，如椎骨、股骨、骨盆、其他骨骼和头部时，这是合适的，对于许多诊断和治疗目的来说，这些可以被认为是刚性的。
关于应用，我们集中于腰椎的配准，并且提供可以用于椎弓根螺钉放置的结果。

在任何新的配准技术广泛应用于临床之前，一个必要的步骤是对方法进行评估和验证[^22]。
评估配准技术的一个困难是需要一个高度精确的“金标准”，最近，我们设计了一个尸体腰椎体模，获得了体模的三维MR和CT图像以及二维X射线图像，并且通过基准点建立了三维和二维图像之间的精确的“金标准”刚性配准标记[^23]。
我们使用这些数据来研究所提出的方法的捕获范围和配准精度。

Ch02 理论

在本节中，我们首先用公式表示X射线图像生成过程，然后给出X射线图像强度梯度和衰减系数梯度之间的关系，这将在我们的配准方法中使用。

图1：X射线图像生成-几何考虑

X射线图像生成

以下模型描述了X射线图像的生成
$$
I(\bold{p})=I_0\cdot\cos^3\vartheta\cdot e^{-\int_L\mu(\bold{r})dr}
$$
位于探测器平面$U$上的点$\bold{p}$处的图像强度$I(\bold{p})$是通过沿着投影束$L$从X射线源$\bold{r}_s$到点$\bold{p}$对X射线衰减系数$\mu(\bold{r})$进行积分而获得的（参见图1）。该因子$\cos^3\vartheta$描述了X射线光束发散的影响。角度$\vartheta$是指投影光束和线$L_0$之间的角度，该角度源自X射线源，并且垂直于探测器表面。$I_0$是线$L_0$与探测器平面相交处的参考图像强度。如果X射线成像系统校正射线光束发散，则获得更简单的模型：
$$
I(\bold{p})=I_0\cdot e^{-\int_L\mu(\bold{r})dr}
$$
对于具有对数静态响应的X射线传感器，强度和衰减系数$\mu(\bold{r})$之间的关系可以用线性模型表示：
$$
I(\bold{p})=a\int_L\mu(\bold{r})dr+b
$$
或者更加简单：
$$
I(\bold{p})=\int_L\mu(\bold{r})dr
$$

图像强度和衰减系数梯度

我们使用等式（4）中定义的X射线图像生成模型来推导衰减系数梯度和图像强度梯度之间的关系。设$\bold{v_u(p)}=\grad_U I(\bold{p})$是投影平面$U$中的图像强度梯度，设$\bold{v_w(p)}=\grad_W I(\bold{p})$是垂直于投影光束$L$的平面$W$内图像强度梯度（参见图1）。梯度$\bold{v_u(p)}$可以被梯度$\bold{v_w(p)}$来表示，反之亦然（参见附录A）。
$$
\bold{v_u(p)}=(\bold{n}\times\bold{v_w(p)})\times\bold{n}
$$

$$
\bold{v_w(p)}=\frac{(\bold{n}\times\bold{v_n(p)})\times\bold{e}}{\bold{n\cdot e}}
$$

其中，$\bold{n}$是投影平面$U$的单位法向，$\bold{e}$是定义在投影光束$L$上的单位向量。

衰减系数梯度$\bold{v_A(r)}=\grad\mu(\bold{r})$可以被分为两个部分（参见附录B）
$$
\bold{v_A(r)}=\grad\mu(\bold{r})=\bold{v_{Ar}(r)+v_{At}(r)}
$$
第一个分量$\bold{v_{Ar}(r)}$平行于投影光束$L$，第二个分量$\bold{A_{At}(r)}$垂直于投影光束$L$。强度梯度$\bold{v_w(p)}$与梯度$\bold{v_{At}(r)}$之间的关系由下式给出（参见附录C）：
$$
\bold{v_w(p)}=\frac1{|\bold{p-r_s}|}\int_L|\bold{r-r_s}|\bold{v_{At}(r)}dr
$$
其中，$|\bold{p-r_s}|$和$|\bold{r-r_s}|$分别是X射线源$\bold{r_s}$到点$\bold{p}$和$\bold{r}$之间的距离。这个关系表明，强度梯度$\bold{v_w(p)}$仅依赖于垂直于投影光束$L$的衰减系数梯度分量$\bold{v_{At}(r)}$。所获得的关系符合已知的事实，即X射线图像描绘了垂直于投影光束的结构变化，而不是沿着光束的结构变化。从等式（8）中显然可知，如果点$\bold{r}$更靠近X射线源，则分量$\bold{v_{At}(r)}$对$\bold{v_w(p)}$的影响更小，从而获得更小的图像强度梯度$\bold{v_u(p)}$。知道了衰减系数梯度$\bold{v_A(r)}$向点$\bold{p}$投影的方式（参见等式8），我们可以将梯度$\bold{v_w(p)}$向X射线源$\bold{r_s}$反向投影，并且获得沿着光束任意位置$\bold{r}$的反向投影梯度$\bold{v_B(r)}$：
$$
\bold{v_B(r)}=\frac{|\bold{p-r_s}|}{|\bold{r-r_s}|}\bold{v_w(p)}
$$
梯度$\bold{v_B(r)}$是放大的梯度$\bold{v_w(p)}$。点$\bold{r}$越靠近X射线源$\bold{r_s}$，则放大的也越多。通过将等式（6）插入到等式（9）中，我们获得了沿着X射线源反向投影的图像强度梯度$\bold{v_u(p)}$：
$$
\bold{v_B(r)}=\frac{|\bold{p-r_s}|}{|\bold{r-r_s}|}\cdot
\frac{(\bold{n\times v_u(p))\times e}}{\bold{n\cdot e}}
$$

Ch03 方法

刚性3D/2D配准涉及寻找刚性变换，该刚性变换将3D（CT或MR）图像带入到与之对应的2D投影（X射线）图像的最佳可能空间关系。如果满足以下两个假设，则可以实现配准目标。首先，X射线成像设备被校准，即对于每个X射线图像，X射线源的精确位置和探测器平面的位置和方向在参考坐标系中是已知的。其次，CT、MR和X射线图像在几何是相对精确的。在临床环境中，可以按照[^24]中所述进行校准。我们用于刚性3D/2D配准的特征是在术前CT或MR体积和术中X射线图像的强度梯度中发现的骨结构表面的法线。特征选择建立在X射线图像中健壮的强度梯度对应于骨结构的假设边界上。

图2：三维/二维配准中使用的带有曲面法线和强度梯度的几何设置

设$\bold{r_i^{S_v}}$是在三维体积（CT或MR）的坐标系$\bold{S_v}$中定义的点，位于三维结构的表面上，设$\bold{v_{Ai}}$是曲面上点$\bold{r_i^{S_v}}$上的法线向量（参见图2）。相同的点在参考坐标系$\bold{S_{ref}}$中的位置由刚性变换$\bold{T}$给出：
$$
\bold{r_i}=\bold{T(r_i^{S_v})}=\bold{Rr_i^{S_v}+t}
$$
其中，$\bold{R}$和$\bold{t}$分别表示坐标系$\bold{S_v}$相对于$\bold{S_{ref}}$的旋转和平移。矩阵$\bold{R}$是由定义的三个角$\bold{\omega}=(\omega_x,\omega_y,\omega_z)^T$表示绕x轴、y轴、z轴旋转，而向量$\bold{t}=(t_x,t_y,t_z)^T$表示沿着x轴、y轴、z轴平移。因此由六个参数$\bold{q}=(t_x,t_y,t_z,\omega_x,\omega_y,\omega_z)^T$表示变换$\bold{T}$。设$\bold{r_s}$是X射线源在参数坐标系统$\bold{S_{ref}}$中的位置。对于由矢量$\bold{q}$定义的三维图像的给定位置，连接$\bold{r_i}$和$\bold{r_s}$，并且通过单位向量$\bold{e_i}$定义方向的直线$L_i$与探测器平面$U$在点$\bold{p_i=p(r_i)}$相交（参见图2）。在点$\bold{p_i}$处的强度梯度是$\bold{v_{ui}=\grad_U I(p_i)}$。

为了测量曲面法向$\bold{v_{Ai}}$和反向投影梯度$\bold{v_{Bi}}$（等式10）之间的对应关系，我们提出了下面的准则函数：
$$
CF=\frac{\sum_{i=1}^N|\bold{v_{Ai}}|\cdot|\bold{v_{Bi}}|\cdot f(\alpha)}
{\sum_{i=1}^N|\bold{v_{Ai}}|\cdot\sum_{i=1}^N|\bold{v_{Bi}}|}
$$
其中，$f(\alpha)$是依赖于梯度$\bold{v_{Ai}}$和$\bold{v_{Bi}}$之间的角度$\alpha$的权重函数，$N$是曲面上点的数目。

函数$f(\alpha)$的定义如下：
$$
\begin{equation}
f(\alpha)=
\begin{cases}
cos^n\alpha,&0<|\alpha|<90^\circ\
0,&其它
\end{cases}
\end{equation}
$$
倾向于具有相似方向的梯度。参数$n$确定权重对角度$\alpha$的灵敏度（参见图3）。在实践中，$\alpha$小于$90^\circ$时权重函数的计算如下：
$$
f(\alpha)=cos^n\alpha=\frac{|\bold{v_{Ai}\cdot v_{Bi}}|^n}{(|\bold{v_{Ai}|\cdot|v_{Bi}}|)^n}
$$
当多玩一个X射线图像用于配准时，准则函数是各个图像的准则函数的总和。通过找到优化准则函数的一组参数来执行刚性的3D/2D配准。

图3：参数$n$的不同值的加权函数$f(\alpha)$

Ch04 评估

“金标准”

因为很难在临床图像上获得准确的“黄金标准”的配准，所以我们获取了一位80岁女性尸体腰椎切片的“金标准”，包括L1-L5椎骨、椎间盘和几毫米的软组织。脊柱被放入一个塑料管中，并且使用细尼龙绳绑在上面。六个基准标记被牢固地附着在试管的外部。管子里面装满了水来模拟软组织。通过这种方式，获得了逼真的MR、CT和X射线图像。CT图像是使用通用电子调整CT/i扫描仪（General Electric HiSpeed）。轴向切片的片内分辨率为$0.27\times0.27$mm，片间距1mm。对于MR图像（MRI），使用飞利浦Gyroscan NT Intera 1.5-T 扫描仪和T1协议（翻转角$90^\circ$，TR=3220ms，TE=11ms）。轴向切片的片内分辨率为$0.39\times0.39$mm，片间距1.9mm。通过信息最小化方法对MR图像的强度不均匀性进行了回顾性校正[^25]。X射线图像由PIXIUM 400（Trixell）数字X射线探测器获得，该探测器具有$429\times429$mm的大有效表面，$0.143\times0.143$mm的像素大小和14位的动态范围。在图像采集期间，X射线源和探测器平面是固定的，而脊柱模型在转盘上旋转，以模拟具有C形臂的设置。通过将脊柱模型绕其长轴旋转（步长=$20^\circ$），获得了18个X射线图像。图4给出了CT、MR和X射线图像。

图4：CT（左上）、MR（右上）、AP（左下）的轴向切片，和带有基准标记的脊柱模型的侧向（右下）X射线图像

CT或MR到X射线图像的“金标准”配准是通过分别将CT或MR的标记点严格配准到从X射线图像重建的三维标记点而获得的。“金标准”配准的准确性通过估计八个目标的目标配准误差（Target Registration Error，TRE）[^26]获得，每个椎弓根上有四个目标，每个目标是人工定义每隔五个椎骨上。所有TRE的均方根（RMS）对于CT到X射线的配准都低于0.29mm，对于MR到X射线的配准都低于0.42mm，这表明“金标准”是高度精确的。关于脊柱模型、X射线系统校准、“金标准”的计算和“金标准”的验证的详细说明参见[^23]。

实现

包含单个椎骨的立方体子体积在CT和MR图像中被手动定义。使用高斯滤波器（$\sigma=0.5$mm）模糊每个子体积，并且各向同性地重新采样到1mm的分辨率。应用Canny边缘探测器和阈值来自动提取对应于骨结构边界的点的位置，并且估计这些点的表面法线方向。结果是每张CT大约19000个边缘点，每张MR大约32000个边缘点。用高斯滤波器（$\sigma=0.5$mm）模糊X射线图，并且应用Roberts边缘探测器计算强度梯度。然后使用等式（6）得到垂直于投影光束的梯度。等式（14）中用于计算准则函数的参数$n=4$。使用Powel方法完成变换参数$\bold{q}=(t_x,t_y,t_z,\omega_x,\omega_y,\omega_z)^T$的优化[^27]。

实验

为了测度该方法的准确性、速度和捕获范围，从围绕“金标准”配准位置的大范围的起始位置和方向进行配准。如果参数空间被归一化，欧几里德度量可以用于计算起始位置相对于“金标准”位置的位移（在参数空间中）。此外，定义起始点的参数值可以在“金标准”位置附近随机选择，以便在区间$[0,D_e]$内实现位移的均匀分布。对于单个椎骨配准，使用以下假设来归一化参数空间：包含大小为80mm的单个椎骨的体积围绕其0.1rad（弧度）（$5.7^\circ$）的中心旋转导致体积点平均平移2mm。

从18个X射线图像的集合中选择18个X射线图像对，使得形成一对的两个图像之间的旋转角度为$80^\circ$。对于每一对，从妕所述产生的25个起始位置进行配准，将CT和MR配准的$D_e$分别设置为18mm（$51.6^\circ$）和9mm（$25.8^\circ$）。这导致每种模式和脊椎骨有450个配准。通过对MR使用较小的$D_e$值，我们能够更加精确地定义MR/X射线配准的捕获范围，该范围预期小于CT/X射线配准的捕获范围。

为了测量配准前后的配准误差，计算八个目标点（每个椎弓根上有四个）的目标配准误差（TRE），作为配准位置和“金标准”位置的目标点之间的距离：
$$
TRE(r)=|\bold{T_r r-T_g r|}
$$
其中，$\bold{r}$是目标点，$\bold{T_r}$和$\bold{T_g}$分别是通过提出的配准方法和“金标准”配准得到的变换。旋转误差通过将TRE分解为平移和旋转分量来确定[^28]。如果八个目标点的TRE都低于2mm，是认为配准成功。

结果

CT/X射线和MR/X射线配准的结果分别参见表1与表2。

表1：CT/X射线在配准前与配准后及成功配准的比例的RMS、MTR（最大目标配准）、RCE（旋转分量误差）

表2：MR/X射线在配准前与配准后及成功配准的比例的RMS、MTR、RCE

单个椎骨的CT到X射线的成功配准（参见表1），TRE的均方根（RMS）范围从L3的0.2mm到L5的0.5mm，并且最大TRE小于1.2mm。角度误差的RMS范围从L4的$0.35^\circ$到L2的$0.48^\circ$，L2获得的最大值为$3.2^\circ$。如果起始位置距离“金标准”的位置不超过6mm或者$17.2^\circ$，则超过$91%$的配准是成功的。如果起始位置在612mm或者$17.2^\circ\sim34.4^\circ$，则$42%\sim56%$的配准是成功的。如果起始位置在1218mm或者$34.4^\circ\sim51.7^\circ$，则不到$10%$的配准成功。

单个椎骨的MR到X射线的成功配准（参见表2），TRE的均方根（RMS）范围是1.11.4mm，角度误差的RMS范围是$0.95^\circ\sim1.6^\circ$，最大值小于$4.5^\circ$。如果位移$d_e$范围是03mm（即$0^\circ\sim8.6^\circ$），则配准成功的概率超过$82%$，除了L1仅有$53%$。对于初始位移范围是36mm（即$8.6^\circ\sim172.^\circ$），则配准成功的概率为$37%\sim59%$。对于初始位移范围是69mm（即$17.2^\circ\sim25.8^\circ$），则配准成功的概率为$19%\sim36%$。

在Pentium III，733-MHz PC上，X光到CT的平均配准时间为20秒，X光到MR的平均配准时间为32秒。

Ch05 讨论和结论

如果将基于图像的配准用于引导外科手术或者外科手术机器人，精度要求尤其严格。外科医生依据术前构建的计划能够达到的精确度直接依赖于配准的精确度。一般来说，必要的精确度变化很大，这取决于手术的类型和病人的解剖结构。因为我们使用了脊柱模型，并且在两个椎弓根上都选择了目标，所以我们准确性最好在脊柱椎弓根螺钉放置的背景下进行判断。最近，椎弓根螺钉放置的误差范围被定义为螺钉轨迹与沿椎弓根中轴的理想路径之间的允许平移和旋转的偏差[^29]。Rampersaud等人报告了L1~L5最大允许平移（旋转）误差：0.65mm($2.1^\circ$)、0.85mm($2.7^\circ$)、1.5mm($4.9^\circ$)、3mm($9.8^\circ$)、3.8mm($12^\circ$)。所提出的方法因为添加了大量的边缘点到准则函数中，所以对于异常值（即不对应于骨表面的边缘）非常稳健。表1中针对各种X射线视图和起始位置给出的配准结果表明，尽管配准误差采用了与[^29]中备有不同的方式定义，但是仍然在上面给出的误差范围内。如果从起始位置进行配准，则超过$91%$的试验配准成功，起始位置是平移小于6mm或者旋转小于$17^\circ$的“金标准”配准位置。这个捕获范围足够宽，可以通过简单的手动配准技术来达到。在使用本配准方法的临床系统中，简单的可视化工具可以帮助操作者进行初始配准以及估计配准精度，通过这个工具，CT或MR中的骨边缘点被叠加在X射线图像上。每当所提出的3D/2D配准不成功时，配准后获得的姿态与正确的姿态相差足够远，使得所有这些情况都可以使用视觉验证检测到。

重要的是，所提出的配准方法允许MR与X射线图像的配准，但是配准误差大于CT与X射线图像的配准，成功的配准误差分别为1.4mm和$1.6^\circ$，这与椎骨L3L5的配准精度要求相当，但与椎骨L1L2所需要的精度不同[^29]。除此之外，成功配准的更少，在某些情况下即使起始姿势接近“金标准”配准也会失败。这很可能是因为多孔的脊柱浸泡在水中以及MRI序列导致了不能在脊柱的骨骼和周围组织之间产生清晰的组织对比。我们发现从MR图像中提取的大约$40%$的边缘点不对应于骨结构的边界。这个问题存在更多解决方案，一种是使用不同的MRI序列，如：Hoad[^30]等人提出的双回波快速梯度回波序列，以获得骨骼和软组织之间良好的组织对比度；另一种是通过对MR边缘图像的一些后处理来减少不对应于骨骼结构的边缘的数量，与大多数配准方法一样，该方法也容易陷入局部极小值。为了限制这种行为，可以应用全局优化算法而不是局部优化算法[^31]，然而这样需要术中时间，而这个时间是受限的。使用视觉检测，我们观察到质量评估MR到X射线配准比CT到X射线配准更困难。投影到X射线图像上的额外边缘点使得视觉检测更加困难和模糊。

权重函数中的参数值$n$（参见等式14）影响了准则函数的平滑度和最终精度，但是关系更加复杂一些。我们已经观察到，小的$n(n=0,1,2)$导致准则函数靠近正确配准处更加平滑（逼近0.5mm和$0.5^\circ$），同时相比远离正确配准处包含了更多的局部极值。当$n=4$时，CT/X射线配准误差要大于0.5mm。如果$n>5$时，则只有表面法线垂直于光线光线的点会添加到准则函数中，这将导致更加粗糙的准则函数。通过试验不同的$n$值后，我们将其设置为4。这样的选择是平衡了准则函数的平滑性接近或者远离正确的配准。

在IGT系统投入临床使用之前，它必须经过严格的验证，包括单个组件的验证、整个系统的验证以及不确定性如何通过整个IGT过程传播的研究。IGT系统组件验证的先决条件，如三维/二维配准，是验证方法的标准化，它包括验证数据集的设计、相应“基本事实”及其准确性的定义、验证协议和验证指标的设计[^22]。我们使用尸体腰椎体模的MR、CT和X射线图像对我们的配准方法进行了充分地验证，通过基准标记建立了“金标准”配准，并且通过目标配准误差评估了其准确性[^23]。然而，由于其他研究人员使用了不同的模型、验证数据集、验证协议、成功配准的标准和误差度量，因此很难将提出的三维/二维配准方法的结果放在以前发表的工作的上下文中。缺乏标准化的验证方法可能是大多数作者没有对他们的方法进行比较评估的原因[^7][^8][^19][^20][^32]。为了表明进行定量比较确定很困难，我们讨论了这些工作，并尽可能就捕获范围和配准精度进行比较。这些工作中使用的方法基于DRR（与我们的类似），依赖于术中X射线投影，不需要分割，并且已经在腰椎上得到了验证。

在[^32]的实验中使用了一个股骨的CT和X射线图像，该股骨刚性地附着在一个精确加工的校准物体上，该物体包含一个嵌入不锈钢球的网格。这些球用于建立“金标准”三维/二维配准。仅从通过平移或者旋转“金标准”的位置获得起始位置进行了5次记录，随机值在$\pm2$mm或者$\pm2^\circ$的范围内。为了评估配准中的误差，他们计算了100个目标点的三维平均TRE，并且全股骨配置得值2.8mm，另一个股骨配置得值0.8mm。对于初始化估计大于$\pm4$mm或者$\pm4^\circ$，它们通常不会接近“金标准”配准。我们得出以下观察结果：首先，干股骨的配准比在水中的软组织的椎骨的配准更加容易；其次，如果旋转的原点未知，并且如果参数空间没有根据在其上定义目标的物体的尺寸进行归一化，则不可能公平地比较捕获范围；第三，目标点放在哪里很重要。在[^32]中，靶点被随机放置在$1\text{dm}^3$的立方体中，而我们将靶点放置在椎弓根周围，这就是将要执行的结构图像引导手术。

Weese等人[^19]和Penny等人[^8]对一个由聚丙烯包裹的腰椎和骨盆组成的人体模型的CT和前后透视图像进行了实验。精确的“金标准”配准来自于成像前附着在模型上的基准标记。[^19]中的起始估计值与[^8]中的起始估计值略有不同。Weese等人从“金标准”配准开始，并且从$(\omega_{Gi}+\Delta\omega,\omega_{Gi},\omega_{Gi}-\Delta\omega;i=x,y,z)$和$(t_{Gz}+\Delta t,t_{Gz},t_{Gz}-\Delta t)$的所有混合中获得了81个不同的起始估计用于旋转参数和投影平面上的高度，其中$\omega_{G}$是 “金标准”的旋转参数。如果旋转$\Delta\omega$设置为0.1、0.15、0.2弧度（rad），并且变换到我们的归一化参数空间，则它们分别对应于3.4、5.2、6.9mm的位移。[^8]中用于配准的64个起始估计是“金标准”值$\pm\Delta P$，其分别为x、y、z旋转参数值$3.4^\circ$、$7.6^\circ$、$7.8^\circ$，平面内的3.6mm和2.4mm平移参数，平面外的50.8mm平移参数。如果平面外平移被忽略，这些值在我们的标准化参数空间中转化为6mm。因此，在[^8][^19]和我们的实验中的起始估计是可以比较的，但是很难比较我们的方法与[^8][^19]中的方法的成功登记率和准确性。Weese等人认为，如果配准后的相似性度量大于相似性度量中最大值的0.75倍，则配准成功。在[^8]中，故障被定义为在六个自由度中的任何一个自由度上，最终配准的位置距离“金标准”的距离大于$\Delta P$。我们认为这两个标准没有本文中的严格，即每个目标点的TRE都小于2mm。使用[^8]中的标准，我们相当一部分误配将被归类为成功；而[^19]中的相似性度量不是一个好的标准，因为它依赖于度量的行为。小的变化可能对应于从“金标准”位置的大的空间位移。对于基于DRR的方法以及我们的方法来说，重要的是在配准开始时，一些对应的特征至少部分重叠。Weese等人和Penney等人通过从DRR和X射线图像中手动指示的对应点确定平行于投影平面的位移，促进了配准。在[^8]和[^19]中，根据“金标准”与最终旋转和平移参数之间的差异给出了配准精度。参数的准确性或者误差不是一个有用的衡量标准，应该衡量的是TRE[^26]。对于不同的源和物体大小，参数值可能非常相似，但是在TRE中可能有显著的差异。在[^19]中，旋转参数对应于$0.5^\circ$以内，基于投影平面的平行移动对应于0.5mm内，而投影平面上方的高度偏差最大小5mm。如果我们忽略5mm的平面外平移，这些值在我们的标准化参数空间中对应于1mm。Penney等人获得的刚体参数的最小误差值是针对没有添加结构的脊柱模型。如果我们忽略4mm的平面外平移误差，误差值在我们的标准化参数空间中转化化0.7mm。这高于我们获得的CT/X射线配准的RMS（TRE）。

总之，对不同配准方法进行公平地比较只能在有限的程度上进行。进行公正的比较评估需要验证数据集、相应的“基线”定义及其准确性、验证协议和验证指标的设计。我们已经在这个方向上迈出了一步，并且获得了这样的数据，这是公开可用的[^23]。

所提出的三维/二维配准方法既不是母粒的基于强度的方法，例如：为CT/X射线配准开发的基于DRR的方法，也不是基于特征的方法，其需要分割术前和术中图像，这是一个严重的缺点。通过实际结合这两种通用的三维/二维配准方法，所提出的方法继承了它们的优点，即速度、准确性和不需要术中分割，从而获得了通用性和适用性。已经对脊柱模型的图像进行了三维/二维配准，并且给出了不同椎骨的结果以获得配准精度的揭示。进一步的实验应集中在以下问题上：将软组织结构、外科手术器械和植入物引入假体影像将如何影响配准，在特定手术中使用哪些投影和多少X射线图像，以及如何减少MR中不对应于骨结构的边缘。然而，所提出的匹配大量（2万~3万）三维表面法线和反投影二维梯度的方法应该对手术器械和植入物产生的异常值具有鲁棒性。

附录A

假设X射线源位于坐标系的原点，探测器平面$U$垂直于z轴。在球坐标$(r,\Theta,\varphi)$中，点$\bold{p}$处的图像强度$I(\bold{p})$和图像强度梯度$\grad I(\bold{p})$的关系如下：
$$
I(\bold{p})=I(\Theta,\varphi)
$$
并且，
$$
\begin{align}
\grad I(\bold{p})
&=\frac{\partial I}{\partial r}\bold{e_r}+\frac1r\frac{\partial I}{\partial\Theta}\bold{e}_{\Theta}+\frac1{r\sin\Theta}\frac{\partial I}{\partial\varphi}\bold{e}\varphi\notag\
&=\frac1r\frac{\partial I}{\partial\Theta}\bold{e}{\Theta}+\frac1{r\sin\Theta}\frac{\partial I}{\partial\varphi}\bold{e}_\varphi
\end{align}
$$
因此，强度梯度位于平面$W$内，该平面垂直于投影光束$L$：
$$
\grad I(\bold{p})=\grad_W I(\bold{p})
$$
考虑到这些关系：
$$
\begin{align}
x=r\sin\Theta\cos\varphi\notag\
y=r\sin\Theta\sin\varphi
\end{align}
$$
探测器平面$U$上的图像强度梯度：
$$
\grad_U I(\bold{p})=\frac{\partial I}{\partial x}\bold{i}+\frac{\partial I}{\partial y}\bold{j}
$$
在球面坐标系中的定义是：
$$
\grad_U I(\bold{p})=\frac{\partial I}{\partial\Theta}[\frac{\sin\Theta\cos\Theta}{r}\bold{e}_r+\frac{\cos^2\Theta}{r}\bold{e}_\Theta]+\frac1{r\sin\Theta}\frac{\partial I}{\partial\varphi}\bold{e}_\varphi
$$
或者，
$$
\grad_U I(\bold{p})=\grad_W I(\bold{p})-(\grad_W I(\bold{p})\cdot\bold{n})
$$
这种关系可以写成：
$$
\grad_U I(\bold{p})=(\bold{n}\times\grad_W I(\bold{p}))\times\bold{n}
$$
或者，
$$
\grad_W I(\bold{p})=\frac{(\bold{n}\times\grad_U I(\bold{p}))\times\bold{e}_r}{\bold{n}\cdot \bold{e}_r}
$$

附录B

在球坐标系中，衰减系数梯度$\mu(\bold{r})$为：
$$
\begin{align}
\grad_\mu(\tilde{r},\Theta,\varphi)
&=\frac{\partial\mu(\tilde{r},\Theta,\varphi)}{\partial\tilde{r}}+[\frac1{\tilde{r}}\frac{\partial\mu(\tilde{r})}{\partial\Theta}\bold{e}_\Theta+\frac1{\tilde{r}\sin\Theta}\frac{\partial\mu(\tilde{r},\Theta,\varphi)}{\partial\varphi}\bold{e}_\varphi]\
&=\grad_r\mu(\tilde{r},\Theta,\varphi)+\grad_t\mu(\tilde{r},\Theta,\varphi)
\end{align}
$$
梯度由两部分组成：$\grad_r\mu(\bold{r})$在投影光束$L$的方向上，$\grad_t\mu(\bold{r})$垂直于光束。

附录C

设$S$是围绕成像物体的封闭表面，$r_1$和$r_2$分别是从坐标系原点到投影光束进入和离开封闭体积的点的距离。球坐标中的等式（4）给出：
$$
I(\bold{p})=I(\Theta,\varphi)=\int_{r_1(\Theta,\varphi)}^{r_2(\Theta,\varphi)}\mu(\tilde{r},\Theta,\varphi)d\tilde{r}
$$
而$I(\Theta,\varphi)$基于$\Theta$的偏微分为：
$$
\frac{\partial I}{\partial\Theta}
=\int_{r_1(\Theta,\varphi)}^{r_2(\Theta,\varphi)}
\frac{\partial\mu(\tilde{r},\Theta,\varphi)}{\partial\Theta}d\tilde{r}
+\mu(\tilde{r_2}(\Theta,\varphi),\Theta,\varphi)\frac{\partial r_2}{\partial\Theta}-\mu(\tilde{r_1}(\Theta,\varphi),\Theta,\varphi)\frac{\partial r_1}{\partial\Theta}
$$
通过假设$\mu(\bold{r})$是一个光滑函数，并且对于存在于$S$上的点$\mu(\bold{r})=0$，我们得到：
$$
\frac{\partial I}{\partial\Theta}
=\int_{r_1(\Theta,\varphi)}^{r_2(\Theta,\varphi)}
\frac{\partial\mu(\tilde{r},\Theta,\varphi)}{\partial\Theta}d\tilde{r}
$$
同样的，$I(\Theta,\varphi)$基于$\varphi$的偏微分为：
$$
\frac{\partial I}{\partial\varphi}
=\int_{r_1(\Theta,\varphi)}^{r_2(\Theta,\varphi)}
\frac{\partial\mu(\tilde{r},\Theta,\varphi)}{\partial\varphi}d\tilde{r}
$$
通过在等式（17）中插入偏导数$\frac{\partial I}{\partial\Theta}$和$\frac{\partial I}{\partial\varphi}$，得：
$$
\grad_W I(\bold{p})=\frac1r\int_{r_1}^{r_2}\tilde{r}\grad_t\mu(\tilde{r},\Theta,\varphi)d\tilde{r}
$$