脊柱手术导航中分步式2D_3D图像配准方法

第19卷第9期2007年9月计算机辅助设计与图形学学报JO U RNAL OF COM PU T ER AI DED DESIGN &COM PU T ER GRA PHI CS V ol 19,No 9Sep ,2007收稿日期:2006-12-18;修回日期:2007-04-23 基金项目:上海市科技攻关项目(05DZ19511) 张 翼,男,1976年生,博士研究生,主要研究方向为计算机辅助手术、医学图像处理 王满宁,男,1977年生,博士研究生,讲师,主要研究方向为计算机辅助手术、医学图像处理 宋志坚,男,1960年生,博士,教授,博士生导师,主要研究方向为计算机辅助手术、医学图像处理、模式识别(zjsong@fudan edu cn)脊柱手术导航中分步式2D 3D 图像配准方法张 翼 王满宁 宋志坚(复旦大学数字医学研究中心 上海 200032)(1145@sohu com )摘要 为实现微创脊柱手术导航,提出一种将脊柱术前CT 与其术中X 射线图像配准的方法 首先基于一种最近点迭代法进行快速粗配准;然后提取图像的梯度特征,根据投影变换原理,采用寻找具有CT 最大梯度投影位置的方法进行精配准 在模拟数据及临床标本上进行实验的配准率分别为92%和78% 实验结果表明:该方法鲁棒性强、人工干预少、适合于临床应用关键词 2D 3D 配准;数字X 射线图像;计算机断层扫描;手术导航中图法分类号 T P391 7Multi Step 2D 3D Image Registration in Image Guided Spine SurgeryZhang Yi Wang Manning Song Zhijian(Digital M edical Research Cen ter ,Fudan Univ ersity ,Sha nghai 200032)Abstract A novel registration method of preoperative CT images and intra operative X ray imag es is pro posed for minimally invasive image guided spine surgery Firstly,an Iterative Closest Point algorithm is ap plied for coarse registration Secondly ,gradient feature is ex tracted from imag e Based on the principle of projective transformation,the fine registration is performed by finding the max imum CT gradient projec tions T he correct rate of registration w as of 92%and 78%respectively on simulated data and clinical phan tom research Ex perimental results show that proposed method is suitable for clinical application w ith good robustness and less degree of human interventionKey words 2D 3D registration;dig ital X ray;computed tomography ;image guided surgery 手术导航是指以CT,M RI 等医学影像数据为基础,通过虚拟现实技术,借助光学定位仪跟踪并显示手术器械相对于病变组织的位置关系,从而实现对手术过程的实时引导 该技术应用于脊柱外科,对提高定位精度、减少术中创伤 降低手术失误率有重要意义 图像空间与病人参考空间配准是手术导航的关键步骤 在基于CT 图像的脊柱手术导航中,欲获得上述配准关系,需要在脊柱表面放置标记物,或暴露脊柱选取表面点,这种做法不利于微创治疗的开展[1] C 型臂X 射线机作为脊柱手术中的常用辅助设备,其图像可以利用视觉标定技术与病人进行无创配准,但由于其为2D 图像,提供的导航信息有限 因此,寻找一种CT 图像(3D)与X 射线图像(2D)的配准方法,间接获得CT 与病人的配准关系,是解决该问题的一条有效途径 目前,2D 3D 图像配准方法大致分为3类[2 3]:1)基于几何特征,通过寻找待配准物体上2组特征间的最短空间距离进行配准 按特征来源分为外加特征与内在特征 前者配准精度高,但因需要侵入性植入标记物,临床应用受到限制;后者需要对边界进行精确的分割,在实际应用中误差较大 2)基于图像密度特征,通常是采用DRR(digitally reconstructed radiographs)算法,利用CT数据重建的X射线图像与实际X射线图像密度的相似性进行配准,其精度较高;但由于该算法耗时、捕捉范围有限且需要划分兴趣区域,因此效率不高 3)基于上述2种特征结合的方法,通过寻找2类图像中对应的几何特征,如梯度特征,利用相似性原理进行配准,可同时获得较高的配准精度与速度;该算法的关键在于几何特征及最优化函数的选取,其捕捉范围较小及鲁棒性较弱是目前存在的问题针对上述特点,本文提出一种基于不同特征的分步配准方法 首先根据几何特征,采用最近点迭代法(iterative closest point,ICP)算法进行快速粗配准,将误差限制到一定范围之内,以解决捕捉范围有限的问题 精配准基于一种我们改进后的图像梯度特征的配准算法,根据CT梯度在成像平面的投影强度构造最优化函数,能获得较稳定的局部极值,减小了X射线图像对配准精度的影响,提高了系统的鲁棒性 另外,由于本文算法只对形成图像边缘的光线进行追踪,因此配准速度也得到了大幅提升1 配准方法在C型臂系统中,X射线由焦点发出,穿越病人后在影像增强器上形成X射线图像(如图1所示) CT图像空间到病人空间坐标转换矩阵经由C 型臂系统间接得到T pat i ent CT =T X rayCTT i ntensifi erX rayT se nsorint e nsi fier(T sensorpat i ent)-1 图1 系统组成及局部空间坐标转换关系其中,T int e nsifierX ray为X射线空间到增强器空间的转换矩阵,可利用视觉标定方法获得;T se nso rint e nsi fier 与T sensorpatient分别为增强器及病人到定位仪空间的转换矩阵,由定位仪直接得到;T X rayCT为CT到X射线图像的空间转换矩阵,是本文研究的核心内容 配准前需要获得2~3幅不同角度的X射线图像,根据病人体位调整CT模型的姿态,使两者角度误差小于10,距离误差小于30m m,以满足粗配准的初始条件1 1 基于几何特征的粗配准ICP算法[4]是解决2类点集配准的常用算法,其基本思想是寻找从浮动点集q到参考点集p的一个转换矩阵T,使得2组点集之间对应点的平均距离达到最小,即T^=arg T min!∀p i-T q i∀(1) 该算法从初始的对应关系C0开始迭代 在每次迭代中,根据参考点的位置寻找最近的浮动点,建立起对应关系C i,用点配准法计算出T i,将T i应用于原浮动点集q后,再寻找新的对应关系C i+1 我们对原始ICP算法进行改动,将参考点设置成为活动状态:在X射线图像中的脊柱轮廓上,随机选取200个点,参考点位于焦点到这些点的射线上 浮动点位于CT脊柱表面:采用Marching Cube算法生成脊柱表面三角网格模型,利用三角片折叠的方法进行网格动态化简,将精简后的三角网格的顶点作为浮动点集 配准过程如下:对于每条射线,在浮动点集中查找与之相切的最近点作为q i,根据q i计算出在该射线上与之距离最近的点p i,按照式(1)进行迭代由于CT表面顶点的数目庞大,本文采用构建顶点最小外接球树的方法[5]加快最近浮动点的搜索速度 首先对所有顶点生成一个最小外接球作为根节点,然后将该球的外接盒分成8份,对每份外接盒中包含的顶点再生成一个最小外接球作为子节点 这样逐渐细分,直到外接球的半径小于设定的阈值(如图2所示) 搜索初始时,对每条射线构造一个以外接球为元素的堆,并按照外接球边界到射线的距离进行排序,堆的根一直是具有最短距离的外接球 我们利用顶点所在三角片的法向量进行相切判断,当其中一个与射线方向相同、另一个相反时,则认为该顶点与射线相切 搜索步骤如下:Step1 初始化最近点v^,其到射线的距离设为#;Step2 构造排序堆H,将外接球树的根节点元素插入H;Step3 while H不为空1)从H中取出根节点元素u;2)如果u的距离大于v^到射线的距离则搜索结束,最近点为v^;3)如果u包含子节点则在H中插入u的全部子节点;否则寻找u中与射线相切的最近点v,根据v到射线11559期张 翼等:脊柱手术导航中分步式2D 3D图像配准方法的距离判断是否替换v ^;end while图2 脊柱表面顶点外接球连续细分的结果为准确提取X 射线图像中对应的脊柱轮廓,本文采用Live Wire 交互式分割方法[6](如图3所示) 针对透视图像噪声多、目标轮廓强弱不均的特点,我们使用高斯函数的一种变体构造相邻像素p 和q 间的代价函数C (p ,q )=1-e -(f -m ean)22 2其中,f 为Canny 算子边缘检测后的像素值;mean 和 是通过分割前人工选取图像上局部较弱边缘作为训练区域,并以该训练区域的f 值统计得出图3 L ive Wire 交互式分割方法提取脊柱轮廓1 2 基于密度梯度特征的精配准1 2 1 原理根据X 射线成像原理,对于到达成像平面上任意像素的X 射线强度,有E (u ,v )=E 0e -∃L (p )d r ;其中, 为组织对X 射线的衰减系数,r 为组织的深度,E 0为射线L 发出时的强度,p 为L 上任意一点,p %是L 与成像平面的交点,u ,v 是成像平面坐标(如图4所示)图4 X 射线梯度投影模型由于接收器对X 射线呈对数响应特征,因此其图像灰度I 与衰减系数之间存在线性关系I (p %)=a∃L(p )d r +b(2)根据CT 值V (p )与衰减系数间的线性关系,由式(2)推导出I (p %)&∃LV (p )d r(3)将p 表示成参数形式p ( )=r s + (p %-r s ), ∋[0~1];设r (p %)=∀p %-r s ∀,则式(3)变为I (p %)&r (p %)∃V (p ( ))d(4)在u 方向对式(4)求导得I (p %) u =lim !u (0I (p %+!u)-I (p %)!u & lim !u (0∃{r (p %+!u)V(r s + (p %+!u -r s )- r (p %)V (r s + (p %-r s ))} {!u }d &r (p %)lim !u (0∃V(p ( )+ !u )-V (p ( ))!u d & r (p %)∃ !V (p ( )) u T d &r (p %)∃ !V (p ( ))d uT(5)同理,I (p %) v&r (p %)∃!V(p ( ))d v T(6)将式(5),(6)合并,得到X 射线图像梯度与CT 梯度投影的关系为!I (p %)=I (p %) u I (p %) v&r (p %)∃!V (p ( ))du vT(7)1 2 2 算法实现一般算法是根据CT 梯度投影与X 射线图像梯1156计算机辅助设计与图形学学报2007年度的相似程度进行最优解搜索,如Wein等[7]采用两者互相关性为测度函数;Tomazevic等[8]提出将X 射线图像梯度在射线上的反投影与CT梯度比较的方法 但我们在研究中发现:X射线图像梯度方向会随射线能量的(kV)不同产生变化,因此文献[7 8]方法鲁棒性不高 根据投射原理,X射线图像边缘的产生是由于射线经过了相应投照物体的边缘,因此在配准情况下,这部分形成边缘的射线也必然与CT模型的边缘相切 在2D图像中,位于图像边缘的像素具有最大的梯度值,根据式(7)的梯度对应关系可以得到结论:在以配准点为中心的局部范围内,CT数据场在这部分射线上的梯度投影必然在配准点上取得最大值 我们定义射线的梯度投影函数为GP(T,p%)=T∃ !V(T-1(p()))d u v T;其中T为CT图像到病人空间的变换矩阵,包含3个平移及3个旋转参数(t x,t y,t z,∀x,∀y,∀z);式中省略了r(p%)加快计算速度通过对上述射线GP值的模求和,将配准转化成一个多参数最优化过程,并避免了X射线图像梯度方向变化对配准的影响,即T^=arg T max!p%∋edges∀GP(T,p%)∀本文采用Pow ell多参数优化算法,在每一维内使用Brent搜索算法[9] 初始搜索方向设为共轭的单位向量,Pow ell与Brent算法的收敛值分别设定为10-3与10-2 梯度场采用CT数据场同高斯微分函数卷积( =0 5mm)获得,以减小非骨骼梯度信息对配准的影响 X射线图像边缘由Canny算子在粗配准中得到的脊柱边界内自动提取 提取的边缘可能包含有其他非脊柱的信息,如手术器械等,由于在CT模型中没有与之对应的物体,所以这部分光线的GP值很小,可以设定一个阈值将其过滤 每次迭代时,对经过图像边缘的每条射线以T-1变换到CT图像空间,并按照光线投射法原理对梯度场采样,采样间隔1mm2 实 验C型臂采用GE公司的OEC9800,空间分辨率为1024)768,灰度分辨率12bit CT的扫描参数如下:准直器宽度2mm,螺距1 0,视野20cm)20cm 光学跟踪设备采用NDI公司的Polaris以及配套的跟踪器械,定位误差∗0 35mm 所有程序运行于Dell650工作站(Xeon2 8GHz,1GB RAM) 我们以腰椎为例设计了2步实验进行验证 配准精度应用表面配准误差及T的各分量误差评估,表面配准误差为sT RE=∀T p CT-T g ol d p CT∀其中,T为本文算法得到的坐标转换关系,T g old为标准坐标转换关系,p CT为CT脊柱模型的表面点 实验结果如表1所示,sT RE∗3mm则认为配准成功表1 配准实验结果步骤sT RE mm!x!y!z!#x!#y!#z时间 s成功率 %实验1粗配准5 503 182 663 311 371 690 9421 6精配准1 360 790 650 570 300 240 4742 392实验2粗配准5 272 792 992 912 091 552 5420 7精配准2 281 081 311 061 140 981 2762 178距离误差单位mm,角度误差单位实验1 为排除C型臂系统标定误差的影响,采用临床病人的脊柱CT数据及由它生成的模拟X 射线图像进行配准 首先人工设定T gold及射线投影参数,然后将CT数据以DRR方式生成2~3幅不同角度的模拟X射线图像,最后随机生成CT模型位置进行配准实验 为了进一步检验精配准算法,我们以CT模型经T g ol d的转换后的标准位置为中心,对6个空间变换方向分别进行最优化函数值统计,如图5所示,纵轴为函数值,横轴分别为沿坐标轴轴旋转及平移量 结果表明:在局部范围内该函数具有单一的极值,并出现在标准位置附近,曲线呈现出较好的收敛特性实验2 在C型臂系统上采用人体脊柱标本模拟实际手术 术前将5个人工标记物附着在标本上进行CT扫描;术中将利用标记物点配准方法获得的空间转换关系作为T g ol d,同时X射线图像经不同方位采集后经B样条校正[10];利用直接线性变换法进行系统标定,获得射线投影参数 图6所示为一配准实例,我们将CT投影到X射线图像上以显示两者配准情况11579期张 翼等:脊柱手术导航中分步式2D 3D图像配准方法3 讨 论本文采用基于不同特征的分步2D 3D配准方法,较好地克服了单一配准方法的局限性 整个配准过程采用简单的交互操作,有利于配准效率的提高 由于本文算法只利用了X射线图像边缘的位置信息,因此放宽了对成像条件的限制,方便与不同类型的C形臂系统配套使用 在配准精度方面,参照国外手术导航系统的临床应用情况[11],配准误差一般要求小于5mm,因此本文算法的配准精度处于临床可接受范围内 在配准成功率方面,虽然模拟实际手术的成功率有所下降,但部分原因是由于人工估计初始位置的角度误差较大造成的,而在实际应用中可以采用根据病人体表解剖标志点进行初始位置估计的方法加以避免[5] 同时,由于手术导航系统要求术前采用标记点定位来验证配准精度,因此在配准误差较大时需要重新进行配准,从而避免了错误的术中定位 总之,本文提出的2D 3D图像配准方法在精度与速度上取得了较好的平衡,并且操作简便,能满足脊柱外科手术导航的临床需要参 考 文 献[1]Zhang Ying,Yuan Wen The current status and development ofcomputer assisted spinal surgery[J] Journal of Spine Surgery,2005,3(1):51-53(in Chinese)(张 颖,袁 文 计算机辅助脊柱外科手术的应用现状及进展[J] 脊柱外科杂志,2005,3(1):51-53)[2]Fi tzpatrick J M,West J B,M aurer C R Predicting error inrigi d body point based registration[J] IEEE T ransactions onM edical Imaging,1998,17(5):694-702[3]Tomazevic D,Likar B,Pernus F 3 D 2 D regi stration by integrating2 D informati on in3 D[J] IEEE Transactions on M edical Imaging,2006,25(1):17-27[4]Besl P J,M ckay N D A method for registration of3 D shapes[J] IEEE Trans actions on Pattern Analysis and M achine Intelli gence,1992,14(2):239-256(下转第1165页)1158计算机辅助设计与图形学学报2007年[4]Kaucic Branko,Zalik Borut An overview of vi sibility problemalgorithms in1,5D[C] Proceedings of the10th Internati onal Conference in Central Europe on Computer Graph ics,Visualiza tion and Computer Vision+2002 Plzen,2002:241-274 [5]Kim Young Hoon,Rana Sanjay,Wise Steve Expl oring multiple view shed analysis using terrain features and optimisation techn i ques[J] Computers and Geoscien ces,2004,30(9 10): 1019-1032[6]Lee Jay,Stucky Dan On applying view shed analysis for determining least cost paths on di gital elevation models[J] Interna tional Journal of Geograph i cal Information S ystems,1998,12(8):891-950[7]Zhang Jinfang,Li Lei,Wang Yuxi n T errain visibility analysis[J] Journal of System Simulati on,2005,17(8):1916-1921 (in Chinese)(张金方,李 磊,王宇心 地形可视性分析[J] 系统仿真学报,2005,17(8):1916-1921)[8]Lu Pin,Zhang Jinfang,Lu M in An optimal method for multiple observers sitting on terrai n bas ed on improved simulated an nealing techniques[C] Proceedings of the19th Internati onal Conference on Industrial,Engineering&Oth er Applications of Appli ed Intelligent Systems,Annecy,2006:373-382[9]Li Zhilin,Zhu Qing Di gital elevation model[M] Wuhan:W uhan U niversity Press,2000:5-59(i n Chinese)(李志林,朱 庆 数字高程模型[M] 武汉:武汉大学出版社,2000:5-59)[10]Tomlin C Dana Geographic information systems and cartographic modeling[M] Englew ood Cliffs:Prentice Hall,1990:85-112[11]Chao Qingyu,Lee Jay,M unro Stasiuk M andy J Extensi ons toleast cos t path algorithms for roadway plann i ng[J] International Journal of Geographical Information Scien ce,2003,17(4):361-376[12]Saha A K,Arora M K,Gupta R P,et al GIS based routeplanning in landslide prone areas[J] Internati onal Journal ofGeographical Information S cience,2005,19(10):1149-1175 [13]Wang Ling Intelligent optimization algorithm and i ts application[M] Beijing:T si nghua University Press,2001:17-33(inChinese)(王 凌 智能优化算法[M] 北京:清华大学出版社,2001:17-33)[14]M ineter M J,Dow ers A,Caldw ell D R,et al High throughput computing to enhance invervisibility analysis[OL] (200309 08)[2006 10 25] http: w 2003 [15]Caldw ell D R,M ineter M J,Dow ers S,e t al Analysis and visualizati on of visibi lity surfaces[OL] (2003 09 08)[2006 1025] http: ww w geocomputation org 2003[16]M ineter M J,Dow ers S,Caldw ell D R A multicomputing software environment for ArcInfo intervi sibility analysis[OL](2003 09 08)[2006 10 25] http: w w w geocomputation org2003(上接第1158页)[5]Livyatan H,Yaniv Z,Joskow icz L Gradient based2D 3D rigidregistration of fluoroscopic X ray to CT[J] IEEE T ransactions on M edical Imaging,2003,22(11):1395-1406[6]Chodorow ski A,M attsson U,Langi lle M,et al Colour lesionboundary detection using live w ire[C] Proceedings of S PIE, San Diego,2005,5747:1589-1596[7]W ein W,Roper B,Navab N 2D 3D registrati on based on volume gradients[C] Proceedi ngs of SPIE,S an Diego,2005, 5747:144-150[8]T omazevic D,Likar B,Pernus F Ri gid2D 3D registrati on ofintraoperative digital X ray images and preoperative CT and M R images[C] Proceedings of SPIE,M il an,2002:507-517[9]Press W H,Flannery B P,Teukolosky S A,et al Numericalrecipes in C[M] 2nd ed Cambridge,UK:Cambridge University Press,1992:412-419[10]Livyatan H,Yaniv Z,Joskow icz L Robust automati c C armcalibration for fluoroscopy based navigation:a practical approach[C] Proceedings of International Conference on M edical ImageComputing and Computer As sisted Intervention,Tokyo,2002:60-68[11]Herring J L,Daw ant B M,M aurer C R Surface based registration of CT images to physical space for image guided s urgery ofthe spine:a sens i tivity s tudy[J] IEEE Transactions on M edical Imaging,1998,17(5):743-75211659期鲁 敏等:基于平均视距的最优可视覆盖路径。