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

一种神经外科手术导航系统中脑组织变形校正的方法专利名称：一种神经外科手术导航系统中脑组织变形校正的方法技术领域：本发明属医学图像处理及应用领域，涉及一种外科手术导航系统精度校正方法，具体涉及一种神经外科手术导航系统中脑组织变形校正的方法背景技术： 临床外科手术中脑组织变形是影响神经外科手术导航系统精度的重要因素目前的解决方法包括术中影像校正如术中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.10866-874，1999.)发表的“Model-Updated ImageGuidanceInitiai Clinical Experience with Gravity-Induced BrainDeformation”一文中提出一种固结理论模型，但由于该模型需要将体力作为边界条件，即需要在术中获得脑脊液的位置，然而脑脊液在临床上是难以测量的，使得该模型目前也无法实现临床应用发明内容本发明的目的是为临床应用提供一种实施简单，操作灵活，无需更改现有导航设备的外科手术导航系统精度校正方法，具体涉及一种神经外科手术导航系统中脑组织变形校正的方法本发明采用基于线弹性理论的物理模型，同时为了方便有效的获取边界条件，采用三维激光成像设备LRS(Laser Range Scanner)实时获取术中裸露的大脑皮层，通过跟踪算法跟踪脑皮层运动，从而获得用于驱动模型的边界条件将该边界条件与线弹性物理模型相结合，可以有效地对术中的脑组织变形进行校正。

本发明方法实施简单，操作灵活，无需更改现有的导航设备，便于临床应用本发明的技术方案是首先采用基于MRI的三维自动分割算法，获得目标区域如脑组织，随后将分割出的脑组织网格化，弹性理论的基础上通过对每一网格单元赋予相应的生物力学属性，建立脑组织的物理模型借助三维激光扫描设备，通过跟踪算法跟踪裸露脑皮层的运动，将其作为边界条件并结合物理模型进行有限元计算，获得整个脑组织任意位置的变形，最后采用一种插回算法更新术前三维数据场用于指导手术本发明方法通过下述步骤实现，1.三维自动分割算法解决脑组织变形的第一步是分割出目标区域如脑组织，本发明采用一种基于灰度直方图原理结合形态学特征的自动分割算法所述算法首先采用高斯曲线拟合直方图曲线，根据拟合高斯曲线自动判定门限值，门限化后的结果可将脑组织与大部分的颅骨相分离所采用的门限化判定公式为ts＝tBF+4/5(uGM-tBF)其中tBF是背景与前景之间的门限，uGM为脑白质灰度均值，分别对应高斯曲线的第一和第二个峰值在门限化后，脑组织和颅骨之间还会存在细微的连接，此时可采用形态学中的腐蚀算法，去除所存在的细微连接；通过设定种子点寻求最大连通域可将脑组织从皮肤、骨骼中分离出来；最后针对分割出的脑组织采用形态学中的膨胀算法，恢复腐蚀操作中所损失的脑组织。

2.网格化针对分割出的脑组织，采用四面体单元进行离散化线弹性模型的特点在于每一单元都是常应变单元，本发明方法针对在脑组织边界、脑室边界处相对内部组织变化较大的现象，采用一种具有多分变率的网格划分算法，边界处单元密度大，而内部单元密度小，能精确表达表面，较好的模拟出边界处的变形情况，而对于内部变化相对较小的区域，单元则划分的较大些，减少了计算量所述网格划分算法结合类八叉树和MT(Marching Tetrahedron)算法，采用case表提高算法效率算法组织如下述，1)采用类八叉树算法将三维图像空间划分为由均匀的六面体单元所组成的网格；2)针对每一六面体单元将其划分为五个四面体单元；3)对2)所获得均匀四面体网格，在边界处加以细化，获得多分变率网格；4)采用类MT算法对网格进行切割，去除背景，获得最终的仅包含脑组织的多分变率网格3.物理建模针对网格化后的脑组织，通过对每一单元赋予相应的生物力学属性建立脑组织的物理模型由于脑组织变形过程缓慢而且形变小，应变与应力成线性关系，本发明将其模拟为基于线弹性理论的弹性体模型的数学表达如下▿2ux+11-2v∂∂x(∂ux∂x+∂uy∂y+∂uz∂z)+Fxμ=0]]>▿2uy+11-2v∂∂y(∂ux∂x+∂uy∂y+∂uz∂z)+Fyμ=0]]>▿2uz+11-2v∂∂z(∂ux∂x+∂uy∂y+∂uz∂z)+Fzμ=0]]>其中u为位移矢量，F为力矢量。

v为泊松比，μ＝E/2(1+v)，E为弹性模量4.边界条件物理模型完成后，采用合适的边界条件驱动模型术中最易观察到的是裸露的大脑皮层，所以可通过获得脑皮层变形前和变形后的位置计算出表面位移将该表面位移作为驱动模型的边界条件本方法只需引入三维成像设备LRS，在开颅后，手术前进行扫描后即可移开，操作简单，无需更改现有导航设备且不占用术中空间，所述设备能够提供丰富的几何和纹理信息用于跟踪皮层的运动脑皮层变形前的初始位置可通过刚体配准算法获得，本发明采用一种简单的坐标转换来实现上述刚体配准算法首先借助于术前的点配准(PBR)算法实现图像空间到参考架空间的变换，然后借助于Polaris实现参考架空间到Polaris空间以及Polaris空间到被跟踪器空间的变换，最后借助于校准器，实现被跟踪器空间到LRS空间的变换经过上述一系列变换，最终可以实现图像空间到LRS空间的变换，从而可以获得脑皮层在LRS空间的初始位置变形后脑皮层位置可通过基于活动轮廓模型的非刚体面配准技术实现所述模型所对应的能量函数定义如下述，E(v)=∫∫w10(r,s)|∂v(r,s)∂r|2+w01(r,s)|∂v(r,s)∂s|2+P(v(r,s))drds]]>其中v(r，s)＝(v1(r，s)，v2(r，s)，v3(r，s))为双参数曲面方程，w10(r，s)为弹性系数函数，w01(r，s)为刚性系数函数，P(v(r，s))为势能函数。

引入时间参数t，该能量函数的极值满足欧拉Euler-Lagrange方程 其中v0(s，r)为脑皮层初始位置，可由刚体配准获得F(v)为图像力，定义如下F(v)=▿(P(v))=▿(11+(▿I)2)]]>其中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为脑皮层的扫描结果具体实施方式实施例11.针对256×256×48的三维MRI数据场采用三维自动分割算法，获得脑组织。

门限值分别对应拟合高斯曲线的第一和第二个峰值，腐蚀元素采用半径为5个象素的球状元素，膨胀元素采用半径为6个象素的球状元素2.采用多分辨率网格化算法，将所分割出的脑组织离散为18485个四面体，节点个数为5410边界处最大四面体为7.5×7.5×7.5mm3(用四面体外接六面体大小衡量)，内部最大四面体为15×15×15mm3；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，获得脑。