颅脑手术可视化中手术场空间与图像空间之间的四点配准算法
作者：李修往 李树祥 周猛 魏民
单位：第一军医大学医学图像全军重点实验室(广州 510515)
关键词：颅脑手术可视化；点配准；最小二乘法；单值分解(SVD)
北京生物医学工程990407 摘 要 手术场空间与图像空间之间的配准是影响颅脑手术可视化系统精度的一个主要因素，给定两个三维空间之间的对应两组点{pi}和{p′i};i=1,2,…，N。两个空间中点之间的关系可由p′i=Rpi+T表示，其中R是旋转矩阵，T是平移矢量。本文提供一种基于矩阵单值分解(SVD)求解R和T的最小二乘解的四点算法，并对算法进行了推导，该算法与迭代法相比，具有速度快、用点少的优点，尤其适用于实际应用的颅脑手术可视化系统。
A Four-Point Registration Algorithm of Image-Space With
, 百拇医药
Surgery-Space in Brain Surgery Visualization System
Li Xiuwang, Li Shuxiang, Zhou Meng, Wei Min
(Key Laboratory of Medical Imaging,the First Military Medical University,Guangzhou 510515)
Abstract
Registration between Surgery-Space and Image-Space is one of the key factors that affect the accuracy of the brain surgery visualization system.Given two sets of 3-D corresponding points{pi} and {p′i},i=1,2,…，N.Relationship between 3-D point sets of two spaces fits p′i=Rpi+T,where R is a rotation matrix,T is a translation vector.A four-point algorithm is presented in this paper for finding the least-squares solution of R and T,which is based on the singular value decomposition (SVD) of matrix,and the derivation of the algorithm is provided too.The algorithm has advantages of less computing time,and requires the least number of registration points,compared with conventional iterative algorithm.Especially,it is applicable in practical brain surgery visualization system.
, 百拇医药
Key words:Brain surgery visulization;Point-based registration;Least-squares;Singular value decompostion(SVD)
0 引 言
颅脑手术可视化技术是90年代国际上兴起的计算机辅助外科的一个重要分支。计算机辅助外科(computer-aided surgery)又称作图像引导外科(image-guideqd surgery)、虚拟外科(virtual surgery)等。它利用病人的CT/MRI等影像数据、重建与显示病灶及其周围的实际三维结构，并将之与手术工具所在的坐标系统一起来，根据实际手术需要对手术工具周围的组织作出相应的显示，对手术方案进行定量化评价，以及对手术过程进行实时导航[1]。它是随着各种影像技术(如CT、MRI、DSA等)的发展，以及计算机技术的进步，在近几年来才发展起来的高新技术。其目的是通过将计算机图像技术引入外科手术，尽量减少对正常组织的损伤，保护关键组织，使手术更快速、更精确、更安全。
, http://www.100md.com
颅脑手术可视化作为计算机辅助外科的一个分支，它强调在不改变医生传统手术习惯的前提下，借助现代影像技术、计算机技术、可视化技术、立体定位技术等，对颅脑外科手术提供帮助。它能够从CT、MRI等设备的二维影像数据精确地重建出颅脑的三维结构，通过交互式的界面，可以从各种角度、用多种方法观察和研究病人大脑的内部结构，找出最佳的手术路径，作出准确的手术方案。在手术过程中，由于手术工具周围的结构可以显示出来，医生可以尽量避开血管和神经，以及重要的功能区。因此可以减小手术创口，扩大手术的适应范围，逐步使过去的一些禁忌区变为适应区[1]。
将实际的病人颅脑数据和在手术过程中病人所在的空间位置联系起来是系统的一个关键环节。通过定位系统可以找出病人的手术场空间和病人图像空间的对应点的空间坐标，作为将两个空间配准的依据。在配准系统中，我们可以利用这两个空间中的已知的对应点的空间坐标求出两个空间之间的坐标变换矩阵，将病人的手术场空间坐标系和病人图像空间坐标系统一起来。这样，如果我们知道了病人的手术场空间中某一点的坐标，利用所求得的坐标变换矩阵，就可以计算出它在病人图像空间中的相应的位置。以CT扫描为例，可以认为病人的手术场空间{pi}和病人图像空间{p′i}之间的坐标变换是刚性变换。两空间存在正交变换矩阵R和平移矢量T，使p′i=Rpi+T,如何已知一定数量的对应点，就可以准确地求出R和T。传统的求R和T的方法是迭代法[2]，时间开销较大，并且需要较多的配准点。为此，本文选用基于矩阵单值分解的最小二乘算法，时间开销是迭代法的三分之一，只需四对对应点的空间坐标，就可以求出它们之间的变换关系。以下就该算法进行详细推导和描述。
, http://www.100md.com
1 病人手术场空间坐标系到病人图像空间坐标系的四点配准算法
1.1 说明和引理
为方便下边算法的推导和描述，此处给出一些说明和一个引理：
说明：{pi}和{p′i}(其中i=1,2,…，N)分别代表病人的手术场空间和病人图像空间的三维坐标点，刚性变换的关系可由p′i=Rpi+T表示，其中R是旋转矩阵，T是平移矢量。Trace( )表示求矩阵的迹(对角线元素之和)，det( )表示求矩阵的行列式值，At中的t表示矩阵的转置[3～4]。
引理：对于任何正定矩阵AAt，和任何标准正交矩阵B，
, 百拇医药
(1)
引理证明：假定ai是A的第i个列向量，则
依许瓦尔兹(Schwarz)不等式[3]
因此，
1.2 算法推导
已知病人手术场空间和病人图像空间的三维坐标点，刚性变换的关系为：
, http://www.100md.com
(2)
其中R是旋转矩阵，T是平移矢量，用最小二乘的方法找到R和T使下式的值最小：
(3)
如果(2)的最小平方解是
和
，那么{p′i}和{p″i=pi+}的中心相同。设{p′i}的中心为p′,{p″i}的中心为p″，有：
, http://www.100md.com
(4)
其中：
(5)
(6)
(7)
令：
(8)
, http://www.100md.com
(9)
有：
(10)
至此，最初的最小平方问题分成了两个部分：
(1)在公式(10)中求使Σ2最小化的
(2)平移向量T依公式(4)和(6)则可以由下式求出：
(11)
现对公式(10)右端进行扩展：
, 百拇医药
这样，使Σ2最小就等价于使下式的值最大：
(12)
其中：
(13)
假定H的SVD结果为：
(14)
其中，U和V是3×3的标准正交矩阵，Λ是3×3的对角矩阵，并且对角线元素非负，设：
(X也是标准正交矩阵)
, http://www.100md.com
(15)
则：
(16)
是对称并且正定的，依照引理，对于任何3×3标准正交矩阵B，均有：
(17)
因此，在所有的3×3的标准正交矩阵中，X使(11)式的F值最大，并且，如果det(X)=+1,X是一个本文要求的旋转矩阵；det(X)=-1,X是一个反射矩阵，这种情况发生在{qi}共面的情况下，在实用的可视化定位系统中，在病人的头皮上用四个不共面的点作为配准点，就可避免这种情况的发生。
, 百拇医药
1.3 算法描述
已知病人的头颅上的四个不共面配准点的手术场空间坐标系的坐标{pi}和CT扫描出的病人图像空间坐标系中对应四点的坐标{p′i}(其中i=1,2,3,4)，下面给出计算和的SVD算法：
(1)由公式(5)和(7)计算出{p′}和{p}，然后再由公式(8)和(9)计算出{qi}和{q′i}；
(2)由公式(13)计算出下面的3×3矩阵：
, 百拇医药
(3)计算H矩阵的SVD： H=UΛVt
(4)计算 X=VUt
(5)计算X的行列式det(X)。如果det(X)=+1,那么=X。
(6)由公式(11)计算平移矢量T。
2 结论与讨论
我们在配准实际病人的手术场空间和病人的图像空间的过程中采用了颅外标记点法，过多的标记点对病人是一种负担，也增加了操作的复杂性。为此，本文提出上述算法，只需四个定标点，就可以求出两个空间的配准矩阵，并且在算法的执行速度上与迭代方法进行了比较。在PentiumⅡ 233MHz微机上，两种方法的时间开销如表1所示，本法的速度是迭代法的近3倍，经过1次后迭代处理，精度也超过纯迭代方法，特别是其用点少的优势，成为实际系统应用的极佳选择。
, 百拇医药
表1 算法运行结果比较 配准点数(个)
所用算法
单值分解(SVD)
迭代法(迭代次数)
4
5.4ms
——
7
4.2ms
11.1ms(12次)
11
3.7ms
, http://www.100md.com
10.5ms(8次)
16
3.9ms
9.4ms(5次)
作者简介：李修往(1972-)，男，第一军医大学生物医学工程系硕士研究生，全军医学图像重点实验室实习研究员，主要研究方向为：计算机图形学和颅脑手术可视化。
参考文献
1 李树祥，等.“九五”国家医学重大项目建议书：颅脑手术可视化技术研究，1995
2 魏 民.颅脑手术可视化的跟踪、立体定位、配准和图像引导技术研究.第一军医大学生物医学工程系硕士论文，1998，1：26～31
3 史荣昌.矩阵分析.第1版，北京：北京理工大学出版社，1996，142～155
4 杨华中，等.数值计算方法与C语言工程函数库.第1版，北京：科学出版社，1996，33～44
(1998-10-05收稿), http://www.100md.com
单位：第一军医大学医学图像全军重点实验室(广州 510515)
关键词：颅脑手术可视化；点配准；最小二乘法；单值分解(SVD)
北京生物医学工程990407 摘 要 手术场空间与图像空间之间的配准是影响颅脑手术可视化系统精度的一个主要因素，给定两个三维空间之间的对应两组点{pi}和{p′i};i=1,2,…，N。两个空间中点之间的关系可由p′i=Rpi+T表示，其中R是旋转矩阵，T是平移矢量。本文提供一种基于矩阵单值分解(SVD)求解R和T的最小二乘解的四点算法，并对算法进行了推导，该算法与迭代法相比，具有速度快、用点少的优点，尤其适用于实际应用的颅脑手术可视化系统。
A Four-Point Registration Algorithm of Image-Space With
, 百拇医药
Surgery-Space in Brain Surgery Visualization System
Li Xiuwang, Li Shuxiang, Zhou Meng, Wei Min
(Key Laboratory of Medical Imaging,the First Military Medical University,Guangzhou 510515)
Abstract
Registration between Surgery-Space and Image-Space is one of the key factors that affect the accuracy of the brain surgery visualization system.Given two sets of 3-D corresponding points{pi} and {p′i},i=1,2,…，N.Relationship between 3-D point sets of two spaces fits p′i=Rpi+T,where R is a rotation matrix,T is a translation vector.A four-point algorithm is presented in this paper for finding the least-squares solution of R and T,which is based on the singular value decomposition (SVD) of matrix,and the derivation of the algorithm is provided too.The algorithm has advantages of less computing time,and requires the least number of registration points,compared with conventional iterative algorithm.Especially,it is applicable in practical brain surgery visualization system.
, 百拇医药
Key words:Brain surgery visulization;Point-based registration;Least-squares;Singular value decompostion(SVD)
0 引 言
颅脑手术可视化技术是90年代国际上兴起的计算机辅助外科的一个重要分支。计算机辅助外科(computer-aided surgery)又称作图像引导外科(image-guideqd surgery)、虚拟外科(virtual surgery)等。它利用病人的CT/MRI等影像数据、重建与显示病灶及其周围的实际三维结构，并将之与手术工具所在的坐标系统一起来，根据实际手术需要对手术工具周围的组织作出相应的显示，对手术方案进行定量化评价，以及对手术过程进行实时导航[1]。它是随着各种影像技术(如CT、MRI、DSA等)的发展，以及计算机技术的进步，在近几年来才发展起来的高新技术。其目的是通过将计算机图像技术引入外科手术，尽量减少对正常组织的损伤，保护关键组织，使手术更快速、更精确、更安全。
, http://www.100md.com
颅脑手术可视化作为计算机辅助外科的一个分支，它强调在不改变医生传统手术习惯的前提下，借助现代影像技术、计算机技术、可视化技术、立体定位技术等，对颅脑外科手术提供帮助。它能够从CT、MRI等设备的二维影像数据精确地重建出颅脑的三维结构，通过交互式的界面，可以从各种角度、用多种方法观察和研究病人大脑的内部结构，找出最佳的手术路径，作出准确的手术方案。在手术过程中，由于手术工具周围的结构可以显示出来，医生可以尽量避开血管和神经，以及重要的功能区。因此可以减小手术创口，扩大手术的适应范围，逐步使过去的一些禁忌区变为适应区[1]。
将实际的病人颅脑数据和在手术过程中病人所在的空间位置联系起来是系统的一个关键环节。通过定位系统可以找出病人的手术场空间和病人图像空间的对应点的空间坐标，作为将两个空间配准的依据。在配准系统中，我们可以利用这两个空间中的已知的对应点的空间坐标求出两个空间之间的坐标变换矩阵，将病人的手术场空间坐标系和病人图像空间坐标系统一起来。这样，如果我们知道了病人的手术场空间中某一点的坐标，利用所求得的坐标变换矩阵，就可以计算出它在病人图像空间中的相应的位置。以CT扫描为例，可以认为病人的手术场空间{pi}和病人图像空间{p′i}之间的坐标变换是刚性变换。两空间存在正交变换矩阵R和平移矢量T，使p′i=Rpi+T,如何已知一定数量的对应点，就可以准确地求出R和T。传统的求R和T的方法是迭代法[2]，时间开销较大，并且需要较多的配准点。为此，本文选用基于矩阵单值分解的最小二乘算法，时间开销是迭代法的三分之一，只需四对对应点的空间坐标，就可以求出它们之间的变换关系。以下就该算法进行详细推导和描述。
, http://www.100md.com
1 病人手术场空间坐标系到病人图像空间坐标系的四点配准算法
1.1 说明和引理
为方便下边算法的推导和描述，此处给出一些说明和一个引理：
说明：{pi}和{p′i}(其中i=1,2,…，N)分别代表病人的手术场空间和病人图像空间的三维坐标点，刚性变换的关系可由p′i=Rpi+T表示，其中R是旋转矩阵，T是平移矢量。Trace( )表示求矩阵的迹(对角线元素之和)，det( )表示求矩阵的行列式值，At中的t表示矩阵的转置[3～4]。
引理：对于任何正定矩阵AAt，和任何标准正交矩阵B，
, 百拇医药
(1)
引理证明：假定ai是A的第i个列向量，则
依许瓦尔兹(Schwarz)不等式[3]
因此，
1.2 算法推导
已知病人手术场空间和病人图像空间的三维坐标点，刚性变换的关系为：
, http://www.100md.com
(2)
其中R是旋转矩阵，T是平移矢量，用最小二乘的方法找到R和T使下式的值最小：
(3)
如果(2)的最小平方解是
和，那么{p′i}和{p″i=pi+}的中心相同。设{p′i}的中心为p′,{p″i}的中心为p″，有：, http://www.100md.com
(4)
其中：
(5)
(6)
(7)
令：
(8)
, http://www.100md.com
(9)
有：
(10)
至此，最初的最小平方问题分成了两个部分：
(1)在公式(10)中求使Σ2最小化的
(2)平移向量T依公式(4)和(6)则可以由下式求出：
(11)
现对公式(10)右端进行扩展：
, 百拇医药
这样，使Σ2最小就等价于使下式的值最大：
(12)
其中：
(13)
假定H的SVD结果为：
(14)
其中，U和V是3×3的标准正交矩阵，Λ是3×3的对角矩阵，并且对角线元素非负，设：
(X也是标准正交矩阵), http://www.100md.com
(15)
则：
(16)
是对称并且正定的，依照引理，对于任何3×3标准正交矩阵B，均有：
(17)
因此，在所有的3×3的标准正交矩阵中，X使(11)式的F值最大，并且，如果det(X)=+1,X是一个本文要求的旋转矩阵；det(X)=-1,X是一个反射矩阵，这种情况发生在{qi}共面的情况下，在实用的可视化定位系统中，在病人的头皮上用四个不共面的点作为配准点，就可避免这种情况的发生。
, 百拇医药
1.3 算法描述
已知病人的头颅上的四个不共面配准点的手术场空间坐标系的坐标{pi}和CT扫描出的病人图像空间坐标系中对应四点的坐标{p′i}(其中i=1,2,3,4)，下面给出计算
和的SVD算法：(1)由公式(5)和(7)计算出{p′}和{p}，然后再由公式(8)和(9)计算出{qi}和{q′i}；
(2)由公式(13)计算出下面的3×3矩阵：
, 百拇医药
(3)计算H矩阵的SVD： H=UΛVt
(4)计算 X=VUt
(5)计算X的行列式det(X)。如果det(X)=+1,那么
=X。(6)由公式(11)计算平移矢量T。
2 结论与讨论
我们在配准实际病人的手术场空间和病人的图像空间的过程中采用了颅外标记点法，过多的标记点对病人是一种负担，也增加了操作的复杂性。为此，本文提出上述算法，只需四个定标点，就可以求出两个空间的配准矩阵，并且在算法的执行速度上与迭代方法进行了比较。在PentiumⅡ 233MHz微机上，两种方法的时间开销如表1所示，本法的速度是迭代法的近3倍，经过1次后迭代处理，精度也超过纯迭代方法，特别是其用点少的优势，成为实际系统应用的极佳选择。
, 百拇医药
表1 算法运行结果比较 配准点数(个)
所用算法
单值分解(SVD)
迭代法(迭代次数)
4
5.4ms
——
7
4.2ms
11.1ms(12次)
11
3.7ms
, http://www.100md.com
10.5ms(8次)
16
3.9ms
9.4ms(5次)
作者简介：李修往(1972-)，男，第一军医大学生物医学工程系硕士研究生，全军医学图像重点实验室实习研究员，主要研究方向为：计算机图形学和颅脑手术可视化。
参考文献
1 李树祥，等.“九五”国家医学重大项目建议书：颅脑手术可视化技术研究，1995
2 魏 民.颅脑手术可视化的跟踪、立体定位、配准和图像引导技术研究.第一军医大学生物医学工程系硕士论文，1998，1：26～31
3 史荣昌.矩阵分析.第1版，北京：北京理工大学出版社，1996，142～155
4 杨华中，等.数值计算方法与C语言工程函数库.第1版，北京：科学出版社，1996，33～44
(1998-10-05收稿), http://www.100md.com