张春雷，戴 丽，刘 宇，李 鹤

(东北大学 机械工程与自动化学院，辽宁 沈阳 110819)

随着计算机技术、图像处理技术与临床医学的深度融合，手术导航系统应运而生，其能够实现手术工具对患者靶点的精确定位，并为医生提供稳定、可靠的可视化导航系统，目前已在诸多外科手术中展现出巨大应用价值和市场潜力.患者配准是手术导航系统的关键技术，其目标是获得患者实体空间与其三维虚拟模型所在图像空间的刚性变换关系，以实现术前数据与术中数据的位置统一，配准精度和效率在很大程度上决定了导航系统优劣和手术成败.实现患者配准需借助分别在图像空间和患者空间中描述的至少三组相同的特征基准点[1]，再利用配准算法计算完成.根据基准点的配对情况，患者配准方法可分为成对基准法和非成对基准法.

成对基准法通常需将易在医学影像中辨识和采集的医学标记物植入患者骨骼内或粘于体表，其配准精度通常可达0.5～1.5 mm[2].Chen等[3]在口腔种植系统中将钛钉作为医学标记物，得到的基准配准误差和目标配准误差分为1.12 mm和1.35 mm，但这种侵入式标记物会使患者的痛苦和心理负担倍增.Lin等[4]设计了一种自定义基准标记并将其置于患者头部表面，提出了一种基于光学跟踪系统的实时自动患者配准方法，平均配准误差约为0.7 mm，然而，采用这种方法时，一旦标记发生移位或其重建模型出现缺损，就会造成较大的配准误差.

非成对基准法通常通过点云数据描述患者空间和图像空间中的共同解剖结构特征，两组点云非一一对应关系.Chen等[5]利用颅颌面骨模型表面点云进行配准时的目标误差可达1mm左右.Kong等[6]以人工拾点方式在图像空间中采集骨模型表面点云来大致描述骨模型的局部外廓，但配准精度低、稳定性差，且操作极为耗时.最常用的三维点云配准算法是由Besl等[7]提出的迭代最近点(iterative closet point, ICP)算法，其核心思想是基于最小二乘原理，迭代寻找两组点云的最优刚性变换关系.但ICP算法还存在诸多缺陷，若两组点云初始位姿相差较大，且没有良好的初始变换值，其配准结果将极易陷入局部最优解；最近点搜索过程严重影响计算效率，由此提出的基于多维二叉树的最近点搜索策略[8]和加速ICP算法[9]虽然大幅提升了最近点搜索速度，但这会额外增加创建树结构的时间.此外，聚类ICP算法[10]、鲁棒尺度ICP算法[11]和主成分分析法[12]等基于复杂点云几何特征的配准算法更适于解决规模大、轮廓特征明显的点云匹配问题.因受实际手术条件限制，很难在患者实体上采集和构建可靠且几何特征丰富的大规模点云，因此，传统ICP算法及其变体算法均难以直接应用于手术导航系统的患者配准中，且需要找到一种兼顾精度高、效率高和可操作性好的患者配准方法.

本文结合外科手术导航系统的应用要求和ICP算法的特点，以提升配准精度和可操作性为目标，提出基于三点法初始配准和ICP算法精确配准的患者配准方法.实验结果表明，该方法简便易行、结果稳定，具备良好的配准精度和配准效果.

1 患者配准问题描述及数据获取

手术导航系统主要分为3个子系统，即光学定位子系统、图像子系统和患者子系统.光学定位子系统(O)选用加拿大NDI公司Polaris Vega被动式光学三维定位系统，其工作组件包括位置传感器和被动式定位工具(探针工具和定位刚体).其中，定位工具上固定4个位于同一平面但不共线的反射标记球，位置传感器能追踪并返回每个标记球的形心坐标值或组合工具中心点的位姿信息.患者子系统(P)为治疗对象，将定位刚体与患者保持固定来配合手术执行.图像子系统(I)用于患者的CT/MRI数据和三维模型的处理及可视化，基于3D Slicer医学图像处理平台实现.各子系统的坐标系转换关系如图1所示.

各子系统坐标系转换矩阵的求解过程如下:

(1)

(2)

1) 将定位刚体固定于患者骨骼并一同进行CT扫描，然后在3D Slicer中利用CT数据重建出其三维模型，得到图像空间I，并提取完整的骨模型表面点云MI，样本数通常可达几万至几十万.

2) 在患者实体空间P中，利用探针工具尖端触碰若干次患者骨骼表面，采集得到点云VP，样本数通常仅能达到几十个，这是由于在实际手术中会受到骨骼解剖结构的可暴露范围限制.

由此，即可将患者配准问题转化为密集目标点云MI(图像空间)与稀疏源点云MP(患者空间)之间的配准问题，这两组点云的样本数相差较大，且轮廓和形状等几何特性差异悬殊.

2 配准方法

2.1 基于三点法的初始配准

将定位刚体(固定于患者骨骼)上的反光标记球作为图像空间I和患者空间P中的共同特征基准.

(3)

(4)

在光学定位子系统中，定位刚体上每个标记球的球心坐标均是在位置传感器坐标系O中实时更新和描述的.与上述方法同理，可获得局部坐

(5)

设点F4在图像空间I和患者空间P的齐次坐标分别为fI和fP，则定义初始配准误差eInit为

eInit=|fI-TInit·fP|.

(6)

由于CT图像存在伪影且三维重建时易导致标记球模型发生畸变或缺损，而依赖于成对基准点的配准方法又对基准点拾取精度要求较高，因此仅将上述方法用于初始配准.

2.2 基于ICP算法的精确配准

设MP和MI的齐次坐标矩阵分别为V′和U:

(7)

(8)

利用初始配准矩阵TInit将源点云V′映射至图像空间中，得到新的源点云V，即

V=TInit·V′.

(9)

执行精确配准的目标是获得源点云V相对于目标点云U的转换矩阵TICP，其核心思想是根据设定的目标函数，基于最小二乘原理迭代估计点云之间的最优匹配关系，方法如下:

1) 输入源点云V={vi,i=1,2,…,m}和目标点云U={uj,j=1,2,…,n}，其中m≪n.

2) 目标点云U的抽样:点云U中样本数庞大且分布密集，但其绝大部分与源点云V是不匹配的，从而成为干扰噪声，不仅会降低配准精度，还会严重增加最近点对搜索过程的计算耗时.因此，对目标点云U进行合理抽样简化是十分必要的.

以源点云V中各点vi为中心，建立以r为半径的球区域Si；保留目标点云U中所有位于区域Si中的样本(nr个)，得到更新后的目标点云U(rh)(h为迭代计算次数)，如图3所示.如果r值设置过大，则目标点云的样本仍会较多；如果r值设置过小，则可能会造成目标点云U中没有任何样本位于区域Si内，从而造成关键样本数据丢失，最终影响配准精度.实际上，可将配准误差e作为球形抽样区域半径r的取值依据，即

r=k·e.

(10)

式中，k为可设定的比例系数(k> 0).

3) 源点云V的异常样本剔除:将无目标点云U样本的球形区域所对应的源点云样本视为异常点，并予以剔除，得到更新后的源点云为V(r)，其样本数为mr.若出现该问题时，除r值设置过小的因素外，则很可能是在患者空间中采样的操作失误.

4) 最近点搜索配对与错误剔除:对于源点云V(r)中的每个点，均在目标点云U(r)中遍历搜索与之欧式距离最小的点，以组成最近对应点对.为提升配准精度，需删除距离大于设定阈值d的异常对应点对，并得到更新后的目标点云P={pl}和源点云Q={ql}(l=1,2,…,N)，此时二者样本数一致.

5) 建立目标误差函数E(TICP,h):

(11)

式中:h为迭代计算次数；N为最近对应点对的数量；Rh为旋转变换矩阵；Ph为平移变换向量.

6) 利用奇异值分析法[13]求解Rh和Ph，以使式(11)中E取得最小值.

7) 判断迭代终止条件:如果e小于预设阈值ε或迭代次数h到达预设最大值hmax，则停止迭代，配准完成；否则，重复步骤2)～6)，执行第(h+1)次迭代，直至e<ε或h=hmax.

3 配准实验案例

3.1 实验设置

以猪股骨和猪髂骨为实验对象进行配准实验.本文涉及的配准算法、数据操作及图形可视化均基于Python与3D Slicer交互编程实现，计算机操作环境为Windows7 x64位系统、Intel(R) Core(TM) i3-6100(3.70 GHz)处理器和8 GB内存(RAM).具体方法如下:

首先，将固定有定位刚体的股骨和髂骨经CT扫描，并利用3D Slicer中重建三维模型；分别提取标记球的三维虚拟模型和骨模型表面点云；获得各标记球心在光学位置传感器和图像空间中的位置；执行初始配准；利用探针工具分别在股骨和髂骨表面上的6个特定区域内各采集4个点(如图4所示)，获得24个源点云样本，图像空间中的目标点云(模型点云)样本数分别为46 270和47 712个；分别利用本文提出的方法和传统ICP算法执行精确配准；在患者空间中的每个采样区内，保证采样数相同，仅改变采样位置，共重复执行30组实验.其中，设置关键参数:hmax=150，r=6 mm，d=3 mm.

3.2 点云配准效果验证

以股骨配准为例，通过预设实验数据的方式比较本文方法和传统ICP算法的点云配准效果.如图5所示，预先给定一个确定的刚性变换矩阵Tdef，并将初始目标点云A0变换为点云A1；在点云A1中提取子集源点云B0，以模拟通过探针工具在患者空间采集的点云数据；给定初始配准矩阵T0，分别利用本文方法和传统ICP算法对目标点云A0和源点云B0进行匹配，在本文方法中，经初始配准后(第一次迭代时)的源点云为B1，此时点云A0经抽样后得到点云A0S，其样本数减少至2 473个；利用最终的配准矩阵将点云A1映射至初始目标点云A0，并得到验证点云A10；观察点云A0和点云A10的全局匹配效果，即表示目标点云A0与源点云B0的实际配准效果.

图6为本文方法与传统ICP算法的配准效果.由图6可知，当采用本文方法时，配准后的目标点云与源点云的位置和姿态趋于一致，重合度较高；但利用传统ICP算法时，两组点云的整体匹配效果明显较差，配准结果陷入局部最优.

3.3 结果分析

利用配准后两组点云中各对应点对的距离平均值作为配准误差度量，各组实验结果如图7所示.当采用本文方法对猪股骨和猪髂骨执行配准时，其最大配准误差分别为1.08mm和1.05mm，平均配准误差分别为(0.83±0.10) mm和(0.86±0.09) mm，由此可见，配准误差平均值及标准差均较小，各组结果分布较平稳；此外，配准运算的耗时较短，其平均值分别仅为(0.034±0.003) s和(0.035±0.004) s.然而，当利用传统ICP算法进行配准时，其最大误差分别达4.31 mm和4.18 mm，平均误差分别为(3.35±0.42) mm和(3.28±0.44) mm，且运算耗时均在0.8s左右，显然，传统ICP算法对应的配准误差及各组结果波动范围均较大，且运算耗时明显更长.由此可见，相较于传统ICP算法，本文所提方法在实现准确、稳定和高效的患者配准应用中更具适用性和可行性.

为进一步说明采用本文方法执行患者配准的结果稳定且具有可重复性，分别在股骨表面前侧、旁侧和后侧设置患者空间采样区，保持其他条件设置相同，对每种采样配置各执行10组配准实验，结果如图8所示.股骨前侧、旁侧和后侧采样配置对应的平均配准误差分别为(0.85±0.09) mm,(0.89±0.09) mm和(0.81±0.10) mm，可见，误差平均值均小于0.90 mm，标准差相近均未超过0.10 mm.由实验数据可知，在患者空间中不同的采样配置对使用本文方法所得的配准结果影响较小，从而能在一定程度上表明本文方法具有较好的稳定性和可重复性.

实际上，配准精度会受到多方面因素的影响，例如球形抽样半径r、最近点对的距离阈值d、迭代次数h和目标点云样本密度等参数，以及源点云采样误差、三维模型重建误差和计算误差等因素.

本文通过控制变量试验，仅考察r值对配准精度的影响.保持股骨配准实验的其余设置不变，仅改变r值，得到其与配准误差e的关系，如图9所示.可见，在不同初始配准误差eInit下的e随r的变化趋势基本相同.以eInit=1.5 mm为例进行详细说明，当r<2 mm时，e大致呈递减趋势，此时目标点云中所保留的样本数较少，e更易受到eInit的影响；当2 mmeIint.同时，r值上限还受到配准运算耗时和预期配准精度的约束.若要求缩短配算耗时，则r值应尽量减小；在算法实现中，由于运算耗时很短，则可令r值按照较小增量Δr变化并遍历计算出各自所对应的实际配准误差，最后再从中筛选出满足预设配准精度的r值范围.

4 结 论

1) 采用本文所提方法进行患者配准具有比传统ICP算法更佳的配准精度和稳定性.

2) 合理抽样简化图像空间点云利于降低点云噪声对配准精度产生的不利影响和缩短运算耗时，因此应根据初始配准误差和运算耗时等约束条件合理设定球形抽样区域的半径r.

3) 通过本文方法解决患者配准问题的可操作性强，临床应用潜力大，既无需在图像空间和患者空间执行规模大且繁琐的采样，也利于在不侵入患者体内植入标记物的前提下保证较高的配准精度.