CN1582863A - 一种神经外科手术导航系统中脑组织变形校正的方法 - Google Patents

Abstract

一种神经外科手术导航系统中脑组织变形校正的方法，属于医学图像处理及应用领域。本发明首先采用基于MRI的三维自动分割算法，获得目标组织(脑组织)，随后将分割出来的脑组织网格化，在线弹性理论的基础上通过对每一网格单元赋予相应的生物力学属性，建立脑组织的物理模型。借助三维激光扫描设备，通过跟踪算法来跟踪裸露脑皮层的运动，将其作为边界条件并结合物理模型进行有限元计算，获得整个脑组织任意位置的变形，最后采用一种插回算法更新术前三维数据场用于指导手术。该方法实施简单，精度可靠，便于临床应用，可集成在现有导航系统中，从而大幅度提高导航系统精度。

Description

一种神经外科手术导航系统中脑组织变形校正的方法

技术领域

本发明属医学图像处理及应用领域，涉及一种外科手术导航系统精度校正方法，具体涉及一种神经外科手术导航系统中脑组织变形校正的方法。

背景技术

临床外科手术中脑组织变形是影响神经外科手术导航系统精度的重要因素。目前的解决方法包括术中影像校正如术中MRI、术中US等和物理模型校正两种方法。其中术中影像校正是精度最高的方法，但缺点是费用昂贵，且容易造成术中感染，故当前研究主要集中在基于物理模型的校正方法上。基于物理模型的方法能够通过脑组织的生物力学属性约束脑组织的运动，计算量小、精度可靠、实施简单，便于临床应用，因此，它是当前国际上的研究热点。Yale大学的ImageProcessing and Analysis课题组在2001年11月的IEEE Workshop onMathematical Methods in Biomedical Image Analysis上提出一种线弹性物理模型，他们将表面位移作为边界条件，并采用双立体相机获取边界条件，这种获取边界条件的方法需要改变现有的导航设备来固定两个相机，同时两个相机的校正、同步使得术中操作更为复杂，因此不能方便地应用于临床。Vanderbilt大学的Biomedical Modeling实验室在《IEEE Transactions on medical imaging》(Vol.18，No.10：866-874，1999.)发表的“Model-Updated ImageGuidance：Initiai Clinical Experience with Gravity-Induced BrainDeformation”一文中提出一种固结理论模型，但由于该模型需要将体力作为边界条件，即需要在术中获得脑脊液的位置，然而脑脊液在临床上是难以测量的，使得该模型目前也无法实现临床应用。

发明内容

本发明的目的是为临床应用提供一种实施简单，操作灵活，无需更改现有导航设备的外科手术导航系统精度校正方法，具体涉及一种神经外科手术导航系统中脑组织变形校正的方法。

本发明采用基于线弹性理论的物理模型，同时为了方便有效的获取边界条件，采用三维激光成像设备LRS(Laser Range Scanner)实时获取术中裸露的大脑皮层，通过跟踪算法跟踪脑皮层运动，从而获得用于驱动模型的边界条件。将该边界条件与线弹性物理模型相结合，可以有效地对术中的脑组织变形进行校正。本发明方法实施简单，操作灵活，无需更改现有的导航设备，便于临床应用。

本发明的技术方案是：首先采用基于MRI的三维自动分割算法，获得目标区域如脑组织，随后将分割出的脑组织网格化，在线弹性理论的基础上通过对每一网格单元赋予相应的生物力学属性，建立脑组织的物理模型。借助三维激光扫描设备，通过跟踪算法跟踪裸露脑皮层的运动，将其作为边界条件并结合物理模型进行有限元计算，获得整个脑组织任意位置的变形，最后采用一种插回算法更新术前三维数据场用于指导手术。

本发明方法通过下述步骤实现，

1.三维自动分割算法

解决脑组织变形的第一步是分割出目标区域如脑组织，本发明采用一种基于灰度直方图原理结合形态学特征的自动分割算法。所述算法首先采用高斯曲线拟合直方图曲线，根据拟合高斯曲线自动判定门限值，门限化后的结果可将脑组织与大部分的颅骨相分离。所采用的门限化判定公式为：

                    t_s＝t_BF+4/5(u_GM-t_BF)

其中t_BF是背景与前景之间的门限，u_GM为脑白质灰度均值，分别对应高斯曲线的第一和第二个峰值。在门限化后，脑组织和颅骨之间还会存在细微的连接，此时可采用形态学中的腐蚀算法，去除所存在的细微连接；通过设定种子点寻求最大连通域可将脑组织从皮肤、骨骼中分离出来；最后针对分割出的脑组织采用形态学中的膨胀算法，恢复腐蚀操作中所损失的脑组织。

2.网格化

针对分割出的脑组织，采用四面体单元进行离散化。线弹性模型的特点在于每一单元都是常应变单元，本发明方法针对在脑组织边界、脑室边界处相对内部组织变化较大的现象，采用一种具有多分变率的网格划分算法，边界处单元密度大，而内部单元密度小，能精确表达表面，较好的模拟出边界处的变形情况，而对于内部变化相对较小的区域，单元则划分的较大些，减少了计算量。

所述网格划分算法结合类八叉树和MT(Marching Tetrahedron)算法，采用case表提高算法效率。算法组织如下述，

1)采用类八叉树算法将三维图像空间划分为由均匀的六面体单元所组成的网格；

2)针对每一六面体单元将其划分为五个四面体单元；

3)对2)所获得均匀四面体网格，在边界处加以细化，获得多分变率网格；

4)采用类MT算法对网格进行切割，去除背景，获得最终的仅包含脑组织的多分变

率网格。

3.物理建模

针对网格化后的脑组织，通过对每一单元赋予相应的生物力学属性建立脑组织的物理模型。由于脑组织变形过程缓慢而且形变小，应变与应力成线性关系，本发明将其模拟为基于线弹性理论的弹性体。模型的数学表达如下：

${\▿}^2{u}_x+\frac{1}{1-2v}\frac{\∂}{\∂x}(\frac{{\∂u}_x}{\∂x}+\frac{{\∂u}_y}{\∂y}+\frac{{\∂u}_z}{\∂z})+\frac{F_x}{\mathrm{\μ}}=0$

${\▿}^2{u}_y+\frac{1}{1-2v}\frac{\∂}{\∂y}(\frac{{\∂u}_x}{\∂x}+\frac{{\∂u}_y}{\∂y}+\frac{{\∂u}_z}{\∂z})+\frac{F_y}{\mathrm{\μ}}=0$

${\▿}^2{u}_z+\frac{1}{1-2v}\frac{\∂}{\∂z}(\frac{{\∂u}_x}{\∂x}+\frac{{\∂u}_y}{\∂y}+\frac{{\∂u}_z}{\∂z})+\frac{F_z}{\mathrm{\μ}}=0$

其中u为位移矢量，F为力矢量。v为泊松比，μ＝E/2(1+v)，E为弹性模量。

4.边界条件

物理模型完成后，采用合适的边界条件驱动模型。术中最易观察到的是裸露的大脑皮层，所以可通过获得脑皮层变形前和变形后的位置计算出表面位移。将该表面位移作为驱动模型的边界条件。本方法只需引入三维成像设备LRS，在开颅后，手术前进行扫描后即可移开，操作简单，无需更改现有导航设备且不占用术中空间，所述设备能够提供丰富的几何和纹理信息用于跟踪皮层的运动。脑皮层变形前的初始位置可通过刚体配准算法获得，

本发明采用一种简单的坐标转换来实现上述刚体配准算法。首先借助于术前的点配准(PBR)算法实现图像空间到参考架空间的变换，然后借助于Polaris实现参考架空间到Polaris空间以及Polaris空间到被跟踪器空间的变换，最后借助于校准器，实现被跟踪器空间到LRS空间的变换。经过上述一系列变换，最终可以实现图像空间到LRS空间的变换，从而可以获得脑皮层在LRS空间的初始位置。

变形后脑皮层位置可通过基于活动轮廓模型的非刚体面配准技术实现。所述模型所对应的能量函数定义如下述，

$E\left(v\right)={\∫\∫w}_{10}(r,s){\left|\frac{\∂v(r,s)}{\∂r}\right|}^2+w_{01}(r,s){\left|\frac{\∂v(r,s)}{\∂s}\right|}^2+P\left(v(r,s)\right)\mathrm{drds}$

其中v(r，s)＝(v₁(r，s)，v₂(r，s)，v₃(r，s))为双参数曲面方程，w₁₀(r，s)为弹性系数函数，w₀₁(r，s)为刚性系数函数，P(v(r，s))为势能函数。

引入时间参数t，该能量函数的极值满足欧拉Euler-Lagrange方程：

其中v₀(s，r)为脑皮层初始位置，可由刚体配准获得。F(v)为图像力，定义如下：

$F\left(v\right)=\▿\left(P\left(v\right)\right)=\▿\left(\frac{1}{1+{(\▿I)}^2}\right)$

其中I为三维激光扫描设备对裸露脑皮层扫描所成的三维表面。

Euler-Lagrange方程可采用FEM(Finite Element Method)求解，从而获得各个节点的变形后位置，结合初始位置，可以获得脑皮层表面节点位移，将该表面位移作为边界条件驱动物理模型。

5.有限元计算

根据获得的脑皮层位移，结合物理模型，采用有限元计算出任意节点处的位移，再结合形函数就可获得脑组织任意位置的变形。

6.更新原始三维数据场

采用插回算法更新术前三维数据场。从变形后的网格单元出发，寻找出单元内的显示坐标点，利用形函数获得该点在未变形前的位置，再利用三线性插值获得该点的灰度值。针对所有变形后的单元进行上述处理，就可以用变形后的三维数据场更新术前三维数据场。

本发明引入三维激光扫描设备获取边界条件，并结合线弹性物理模型预测整个脑组织变形情况，既保证了模型预测精度，又解决了边界条件难以测量这一难题，从而可以在临床上实施，大幅度提高手术导航系统的精度。

附图说明

图1为解决脑组织变形的流程图。

图2为基于256×256×48MRI数据场的三维自动分割算法结果。

图3为针对256×256×48MRI数据场多分辨率网格化的结果，

其中，共划分18485个四面体单元，5410个节点；左图为两维断层图像，中间为该断层的网格化结果，右图为该三维数据场的网格化结果。

图4为刚体配准算法，采用坐标变换实现。

图5为边界条件，其中，左图为裸露脑皮层，用于跟踪术中脑皮层运动；中间为不加约束的脑皮层；右图为脑皮层底部，根据临床经验可认为固定不动。

图6为三维变形结果，其中，左图是未变形前的脑组织，右图是重力所引起的变形后的脑组织。

图7为两维变形结果，其中，左图是未变形前的脑组织断层，右图是变形后的相应断层。

图8为位移场的三维可视化结果，其中，箭头表示位移方向，颜色表示位移大小，颜色由红到篮表示位移由小到大。

图9为3D Digital公司的三维激光扫描设备，型号为Model 200，其中200型有一个彩色摄像机，可以将彩色照片和数据点云叠加，有效测量范围是200mm至750mm。仪器在距离被测物体450mm处的测量精度是±125um。

图10为脑皮层的扫描结果

具体实施方式

实施例1

1.针对256×256×48的三维MRI数据场采用三维自动分割算法，获得脑组织。门限值分别对应拟合高斯曲线的第一和第二个峰值，腐蚀元素采用半径为5个象素的球状元素，膨胀元素采用半径为6个象素的球状元素。

2.采用多分辨率网格化算法，将所分割出的脑组织离散为18485个四面体，节点个数为5410。边界处最大四面体为7.5×7.5×7.5mm³(用四面体外接六面体大小衡量)，内部最大四面体为15×15×15mm³；

3.对每一单元设置脑组织生物力学属性参数。杨氏模量＝3Kpa，泊松比＝0.45；

4.采用坐标变换实现刚体配准，获得脑皮层在LRS空间的初始位置。跟踪器固定在LRS上，Polaris通过跟踪器跟踪LRS，实现LRS在Polaris监控范围内任意移动。

5.在开颅后，手术前利用LRS进行一次扫描，LRS距离病人开颅位置的距离控制在450mm左右，扫描时间5～7秒，分辨率512(水平)×500(垂直)。根据设备提供的变形后的三维表面信息以及由步骤4获得的变形前三维表面信息，结合非刚体面配准技术，跟踪脑皮层运动，获得表面节点位移。

6.采用有限元法将物理模型方程化为矩阵形式：Au＝b(A为总刚矩阵，u为待求位移矢量，b为节点力矢量)，采用共扼梯度法求解含有5410(节点数)×3(维数)＝16230个方程的线性方程组。引入由步骤5所获得的边界条件，消除矩阵A的奇异性，求出位移矢量u，获得脑组织内部任意节点处的位移，再结合有限元中的形函数就可以插值出脑组织内部任意位置的变形。

7.由变形后的数据场，根据形函数计算出所有对显示有贡献的坐标点(整数坐标点)变形前的位置，利用三线性插值计算出该点的灰度值。最后对所有这些整数坐标点组成的三维数据场采用经典的Raycasting算法加以可视化，用于指导手术。

在奔IV2.6G、1G内存、Win2000操作系统下运行，需要2分钟左右。

Claims (8)

1.一种神经外科手术导航系统中脑组织变形校正的方法，其特征是包括下

述步骤：(1)采用基于MRI的三维自动分割算法，获得目标组织；

(2)采用多分辨率网格算法对上述目标组织进行网格化；

(3)对每一网格单元赋予相应的生物力学属性，建立目标组织的物理模型；

(4)采用三维激光扫描设备，通过脑皮层跟踪算法，获得边界条件；

(5)结合边界条件与物理模型进行有限元计算，获目标组织任意位置的变形；

(6)采用插回算法用计算出的变形更新术前三维数据场。

2.根据权利要求1所述的神经外科手术导航系统中脑组织变形校正的方法，其特征是，所述步骤中目标组织是脑组织。

3.根据权利要求1所述的神经外科手术导航系统中脑组织变形校正的方法，其特征是，其中步骤(1)的三维自动分割算法为，采用高斯曲线拟合直方图曲线，根据拟合高斯曲线自动判定门限值，所述的门限化判定公式为，

                   t_s＝t_BF+4/5(u_GM-t_BF)

其中t_BF是背景与前景之间的门限，u_GM为脑白质灰度均值，分别对应高斯曲线的第一和第二个峰值。

4.根据权利要求1所述的神经外科手术导航系统中脑组织变形校正的方法，其特征是，其中步骤(2)的多分辨率网格算法为，用类八叉树算法将三维图像空间划分为六面体单元的网格后，将每一六面体单元划分为五个四面体单元，对所获得的四面体网格边界处细化，得多分变率网格，采用类MT算法切割网格，去除背景，获得最终多分变率网格。

5.根据权利要求1所述的神经外科手术导航系统中脑组织变形校正的方法，其特征是其中步骤(3)的物理模型的数学表达为，

${\▿}^2{u}_x+\frac{1}{1-2v}\frac{\∂}{\∂x}(\frac{{\∂u}_x}{\∂x}+\frac{{\∂u}_y}{\∂y}+\frac{{\∂u}_z}{\∂z})+\frac{F_x}{\mathrm{\μ}}=0$

${\▿}^2{u}_y+\frac{1}{1-2v}\frac{\∂}{\∂y}(\frac{{\∂u}_x}{\∂x}+\frac{{\∂u}_y}{\∂y}+\frac{{\∂u}_z}{\∂z})+\frac{F_y}{\mathrm{\μ}}=0$

${\▿}^2{u}_z+\frac{1}{1-2v}\frac{\∂}{\∂z}(\frac{{\∂u}_x}{\∂x}+\frac{{\∂u}_y}{\∂y}+\frac{{\∂u}_z}{\∂z})+\frac{F_z}{\mathrm{\μ}}=0$

其中u为位移矢量，F为力矢量。μ＝E/2(1+v)，v为泊松比，E为弹性模量。

6.根据权利要求1所述的神经外科手术导航系统中脑组织变形校正的方法，其特征是，其中步骤(4)所述的脑皮层跟踪算法为，

首先采用刚体配准算法获取脑皮层初始位置，算法借助于坐标变换实现，然后采用基于活动轮廓模型的非刚体面配准技术跟踪脑皮层运动，所述活动轮廓模型对应的能量函数定义是，

$E\left(v\right)={\∫\∫w}_{10}(r,s){\left|\frac{\∂v(r,s)}{\∂r}\right|}^2+w_{01}(r,s){\left|\frac{\∂v(r,s)}{\∂s}\right|}^2+P\left(v(r,s)\right)\mathrm{drds}$

其中v(r，s)＝(v₁(r，s)，v₂(r，s)，v₃(r，s))为双参数曲面方程，w₁₀(r，s)为弹性系数函数，w₀₁(r，s)为刚性系数函数，P(v(r，s))为势能函数；

引入时间参数t，所述能量函数的极值满足欧拉Euler-Lagrange方程，

其中v₀(s，r)为脑皮层初始位置，由刚体配准算法获得；

其中F(v)为图像力，定义是，

$F\left(v\right)=\▿\left(P\left(v\right)\right)=\▿\left(\frac{1}{1+{(\▿I)}^2}\right)$

其中I为三维激光设备扫描成像。

7.根据权利要求1所述的神经外科手术导航系统中脑组织变形校正的方法，其特征是，其中步骤(5)以脑皮层运动为边界条件，采用有限元方法求解偏微分方程形式的模型方程，获得所有单元节点位置的变形，再结合形函数计算出脑组织任意位置的变形。

8.根据权利要求1所述的神经外科手术导航系统中脑组织变形校正的方法，其特征是，其中步骤(6)所述的插回算法为，从变形后的网格单元出发，寻找出单元内的显示坐标点(整数坐标点)，利用形函数获得该点在未变形前的位置，再利用三线性插值获得该点的灰度值，对变形后的单元进行上述处理后，可用变形后的三维数据场更新原始三维数据场。

Images