CN105310776A - 一种基于子块的软组织表面变形追踪方法 - Google Patents

Abstract

本发明属医学图像处理及应用领域，涉及一种软组织表面变形追踪算法。本发明方法包括软组织分割算法；对提取的目标组织进行网格化处理，获得初始软组织表面点集；通过三维激光扫描仪或者术中三维成像设备获取变形后软组织表面点集；使用刚性配准方法将对变形后软组织表面和初始软组织表面进行初配准；最后基于子块式能量函数最小非刚性配准算法来获得两个点集中点与点之间的映射关系。本发明的方法实施简单，精度可靠，可集成在现有导航系统中，实现术中软组织变形矫正，从而大幅度提高导航系统精度，有助于临床应用。

Description

一种基于子块的软组织表面变形追踪方法

技术领域

本发明属于医学图像处理及应用领域，涉及一种基于子块的软组织表面变形追踪方法。该方法可以用于手术中软组织表面的变形追踪，继而实现术中软组织变形矫正，大幅度提高手术导航系统精度。

背景技术

研究显示，软组织变形常发生于图像引导手术过程当中。以神经外科导航手术为例，手术过程中，由于开颅后脑脊液的流出、颅内压的改变，以及术中操作和重力等因素的影响会引起脑组织变形，术中脑变形会引起导航定位误差，降低导航精度，导致术后肿瘤残留。生物力学模型利用软组织的生物力学属性约束软组织的运动特征，并借助有限元方程求解组织的变形情况，其特点是操作方便，价格低廉，且效果显著。

采用生物力学模型动态模拟术中软组织变形，关键在于找到合适的边界条件驱动模型，才能准确推导出整个软组织的变形情况。目前生物力学模型通常采用软组织表面变形作为边界条件，软组织表面变形本质上是获取组织有限元模型边界上网格节点变形前后的位移。这样，术中软组织生物力学模型的边界条件的获取就可归结为变形前后的两组点集的非刚性配准问题。

通常，点集配准的关键问题可概括为获取两个点集的相关性，基于相关性得到点集之间的几何变换，从而实现点集的姿势匹配或叠加。ICP(iteratedclosestpoint)算法，是一种刚性配准方法，通过不断迭代在目标曲面上找到与源曲面点距离最近的点。TPS-RPM(thinplatesplinerobustpointmatching)算法利用软分配、确定性退火以及TPS来实现空间映射，利用奇异值丢弃方法求解两个点集之间的相关性和映射参数。该方法有如下缺陷：(1)点集数据量不能过大，即所取软组织表面不能过大，过大则无法收敛；(2)要求软组织表面变形前后的曲面形状一致，尽可能保持对齐；另外还有高斯混合模型(GMM:GaussianMixturemodels)和一致点漂移(CPD：CoherentPointDrift)点集配准算法。所述GMM算法是将模型点集和场景点集分别使用高斯混合模型进行表示，通过最小化两个高斯混合集之间的偏差得出配准结果；所述CPD算法是考虑两个点集的疏密程度，将较密集的点集使用高斯混合模型表示，较稀疏的点集看作密集点集的采样点集，从而将点集配准问题转化为根据稀疏点集的概率分布估计密集点集概率分布的混合密度估计问题。所述两种方法的缺点在于对于软组织表面边界处的点集的匹配效果欠佳。因此，实践中需要一种鲁棒性好的软组织表面变形追踪方法，准确地获取软组织表面变形，驱动软组织生物力学模型，矫正术中软组织变形，提高导航精度。

与本发明有关的参考文献有：

[1]P.J.BeslandN.D.Mckay,“AMethodforRegistrationof3DShapes,”IEEETrans.PatAnalandMachIntel,vol.14,pp.239-56,1992.

[2]H.ChuiandA.Rangarajan,ANewPointMatchingAlgorithmforNonrigidRegistration.ComputVisImageUND,vol.89,pp.114-141,2003.

[3]A.MyronenkoandX.Song,PointSetRegistration:CoherentPointDrift.IEEETransPatAnalandMachIntel,vol.32,pp.2262-75,2010.

[4]B.JianandB.C.Vemuri,RobustPointSetRegistrationUsingGaussianMixtureModels.IEEETransPatAnalandMachIntel,vol.33,pp.1633-45,2011.。

发明内容

本发明的目的在于提供一种软组织表面变形追踪方法，具体涉及一种基于子块的软组织表面变形追踪方法，该方法能准确获取软组织表面变形，可作为边界条件用于驱动软组织生物力学模型，矫正术中软组织变形，提高导航精度。

本发明方法包括软组织分割算法；对提取的目标组织进行网格化处理，获得初始软组织表面点集；通过三维激光扫描仪或者术中三维成像设备获取变形后软组织表面点集；使用刚性配准方法将对变形后软组织表面和初始软组织表面进行初配准；最后基于子块式能量函数最小非刚性配准算法获得两个点集中点与点之间的映射关系。

具体的，本发明的一种基于子块的软组织表面变形追踪算法，其特征在于，其包括步骤：

(1)采用分割算法从术前图像中提取目标软组织；

(2)对提取的目标组织进行网格化处理，并提取表面网格节点作为软组织初始表面；

(3)通过三维激光扫描仪或者术中三维成像设备获取变形后软组织表面点集，作为变形后的软组织表面；

(4)使用刚性配准方法将对变形后软组织表面和初始软组织表面进行初配准；

(5)基于子块式能量函数最小非刚性配准算法获得两个点集中点与点之间的映射关系。

本发明中，所述的目标组织为脑组织或者前列腺组织。

本发明中，追踪软组织表面变形的第一步是从术前图像中分割出目标组织，如脑组织或者前列腺组织。本发明的实施例中采用一种自动分割与手动分割相结合的分割方法，首先采用自动分割，然后用手动分割方法进行完善；本发明方法的步骤(1)中所述的自动分割算法首先通过灰度直方图获取图像的上、下限灰度值以及组织与背景的粗略门限值t。；然后利用门限值t估计图像中组织重心的大体位置，并估计组织的大致尺寸(以球体表示)；之后将组织表面建模成离散三角网格曲面；初始模型为离散三角网格球面，球心位于组织重心，半径为估计组织半径的1/2；最后将初始球面离散三角网格缓慢变形，每次变形一个顶点，当顶点变形至脑组织边界时，需遵循变形力以保证组织表面光滑。

本发明中，针对提取出目标组织进行网格化处理，采用四面体或者六面体网格对目标组织进行离散化，之后提取网格表面的节点作为初始软组织表面。

本发明中，当软组织发生变形后，采用可被光学跟踪设备跟踪的三维激光扫描仪扫描软组织暴露表面，获取变形后软组织表面点集，此方法适用于神经外科开颅手术中；若是在手术中软组织表面未暴露出来，例如在前列腺穿刺手术中，则可以采用术中三维成像设备，例如术中三维B超，获取变形后软组织整个表面点集。

本发明中，获得初始软组织表面和变形后软组织表面后，刚性配准方法将对变形后软组织表面和初始软组织表面进行初配准；对于脑组织表面变形追踪，使用不同空间的坐标系转换实现刚性配准初始软组织表面和变形后软组织表面，包含以下四个空间的坐标变换：(1)将三维激光扫描仪空间变换到跟踪工具空间；(2)将跟踪工具空间变换到跟踪设备空间；(3)将跟踪设备空间变换到参考架空间；(4)将参考架空间变换到图像空间，图像空间为变形前软组织所在空间；其中跟踪工具固定在三维激光扫描仪上，跟踪工具可以被跟踪设备跟踪；参考架固定在患者身体上，参考架也可以被跟踪设备跟踪；上述第一个空间坐标变换通过术前标定过程实现；第二个和第三个空间坐标变换通过跟踪设备的术中跟踪实现；第四个空间坐标变换通过对术前图像和病人空间的标记点刚性配准实现；对于前列腺组织，使用迭代最近点算法(ICP：IterativeClosestPoint)将变形后软组织表面和初始软组织表面进行初配准。

本发明中，完成变形后软组织表面和初始软组织表面进行初配准之后，基于子块式能量函数最小非刚性配准算法获得两个点集中点与点之间的映射关系，其子块式表面跟踪算法包含步骤：

(1)将需要跟踪的软组织初始表面网格节点X＝{x,i＝1,2,……M}划分成若干小块：P＝{p,k＝1,2,……L}，除了小块p_L以外,每个小块p_k包含l个节点，小块p_L包含M-(L-1)*l个节点；

(2)为每个点集小块p_k构建能量最小方程获得点其表面位移；

$\underset{C,f}{\min }E(C,f)=\underset{C,f}{\min }\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{o=1}{\overset{l}{\mathrm{\Σ}}}{c}_{\mathrm{jk}}{\left|\right|y_j-f\left(p_o\right)\left|\right|}^2+\mathrm{\λ}\∫\∫[{\left(\frac{{\∂}^{2}f}{{\∂p}^2}\right)}^2+2{\left(\frac{{\∂}^{2}f}{\∂p\∂y}\right)}^2+{\left(\frac{{\∂}^{2}f}{{\∂y}^2}\right)}^2]$

其中p_k＝{p,o＝1,2,……l}为点集小块，Y＝{y,j＝1,2,……N}为点云的点集，c_jk是定义两个点集之间对应概率的相关性矩阵，λ为权重因子；用薄板样条(Thin-plateSpline)非刚性配准算法计算变形皮层(Y)与未变形皮层(p_k)之间映射函数f，从而获得点集小块表面位移Disp_pk；

(3)集合每个点集小块的位移Disp_pk，获得所有需要跟踪的软组织表面网格节点的位移Disp_x。

本发明的软组织表面变形追踪方法与现有技术相比，具有以下优点：

(1)可以处理大数据点集，即不仅能跟踪局部软组织表面变形，也能够跟踪整个软组织表面变形；

(2)鲁棒性好，不要求软组织表面变形前后的曲面形状一致；

(3)对软组织表面边缘处的跟踪效果更好。

本发明方法能准确获取软组织表面变形，可作为边界条件用于驱动软组织生物力学模型，矫正术中软组织变形，提高导航精度。

附图说明

图1是基于240×240×197MRI脑组织数据场的三维分割算法结果，其中三幅图分别从左到右为横断面，矢状面和冠状面的分割结果。

图2是针对240×240×197MRI脑组织数据场网格化的结果，其中，共划分162650个四面体单元，30150个节点。

图3显示了刚体配准算法的结果，其中点云为变形后软组织，网格为变形前软组织，其中，1为脑组织的结果，通过空间坐标变化得到，2为前列腺组织的结果，通过ICP算法得到。

图4是子块式表面跟踪算法的结果，其中，1为脑组织开颅部位的表面跟踪结果，2为前列腺组织整个表面跟踪的结果。

具体实施方式

实施例1试用于手术中软组织表面的变形追踪

针对240×240×197的三维MRI脑组织数据场采用自动分割和手动分割相结合的算法提取出脑组织，分数灰度阈值设定为0.5；

采用四面体网格将提取出的脑组织离散为162650个网格单元，节点数为30150，四面体最大为5×5×5mm³，提取网格表面节点；

在三维激光扫描仪上安装跟踪工具，通过标定，获取三维激光扫描仪空间到跟踪工具空间坐标变换关系，将参考架固定在患者头部，参考架可以被跟踪设备跟踪到，将三维激光扫描仪安装在一个具备标准自由度的三脚架上，调整三脚架使其采集镜头位于开颅皮层表面法线方向30到60cm处，扫描获取开颅区域变形后软组织表面点云；扫描同时，借助光学跟踪设备获得跟踪工具空间到跟踪设备空间坐标变换关系以及跟踪设备空间到参考架空间的坐标变换关系，通过对术前图像和患者空间的标记点刚性配准得到参考架空间到图像空间坐标变换关系；

获得以上四个空间的坐标关系后，将变形后软组织表面和初始软组织表面配准到统一空间坐标系下，完成初配准；

通过人工勾画，得到开颅区域部位的网格表面节点(变形前软组织表面)和变形后脑组织皮层点云(变形后软组织表面)，其中网格表面节点包含134节点，点云包含18127节点；

将网格表面节点划分为17个小块，前16小块，每个小块包含8个节点，最后一个小块包含6个节点；

为每个点集小块p_k构建能量最小方程获得点其表面位移；

$\underset{C,f}{\min }E(C,f)=\underset{C,f}{\min }\underset{j=1}{\overset{N}{\mathrm{\Σ}}}\underset{o=1}{\overset{l}{\mathrm{\Σ}}}{c}_{\mathrm{jk}}{\left|\right|y_j-f\left(p_o\right)\left|\right|}^2+\mathrm{\λ}\∫\∫[{\left(\frac{{\∂}^{2}f}{{\∂p}^2}\right)}^2+2{\left(\frac{{\∂}^{2}f}{\∂p\∂y}\right)}^2+{\left(\frac{{\∂}^{2}f}{{\∂y}^2}\right)}^2]$

其中p_k＝{p,o＝1,2,……l}为点集小块，其中前16个小块，l＝8,最后一个小块，l＝6，Y＝{y,j＝1,2,……N}，为点云的点集，其中N＝18127，c_jk是定义两个点集之间对应概率的相关性矩阵，λ为权重因子，用薄板样条(Thin-plateSpline)非刚性配准算法计算变形皮层(Y)与未变形皮层(p_k)之间映射函数f，获得点集小块表面位移Disp_pk；

集合每个点集小块的位移Disp_pk，获得所有需要跟踪的软组织表面网格节点的位移Disp_x。

Claims (10)

1.一种基于子块的软组织表面变形追踪方法，其特征在于，其包括步骤：

(1)采用自动分割算法从术前图像中提取目标软组织；

(2)对分割图像进行网格化处理，并提取表面网格节点作为软组织初始表面；

(3)通过三维激光扫描仪或者术中三维成像设备获取变形后软组织表面点集，作为变形后的软组织表面；

(4)使用刚性配准方法将对变形后软组织表面和初始软组织表面进行初配准；

(5)基于子块式能量函数最小非刚性配准算法来获得两个点集中点与点之间的映射关系。

2.根据权利要求1所述方法，其特征在于，所述软组织为脑组织或者前列腺组织。

3.根据权利要求1所述方法，其特征在于，所述步骤(1)中，所述自动分割算法包括步骤：

(1)通过灰度直方图获取图像的上、下限灰度值以及组织与背景的粗略门限值t。；

(2)利用门限值t估计图像中组织重心的大体位置，并估计组织的大致尺寸，以球体表示；

(3)将组织表面建模成离散三角网格曲面，初始模型为离散三角网格球面，球心位于组织重心，半径为估计组织半径的1/2；

(4)将初始球面离散三角网格缓慢变形，每次变形一个顶点，当顶点变形至脑组织边界时，遵循变形力以保证组织表面光滑。

4.根据权利要求1所述方法，其特征在于，所述步骤(2)中的三维网格是四面体网格，或是六面体网格。

5.根据权利要求1所述方法，其特征在于，所述步骤(3)中的术中软组织表面点集通过三维激光扫描仪获取，或通过术中三维成像设备获取。

6.根据权利要求1所述方法，其特征在于，所述步骤(4)中的刚性配准方法通过不同空间的坐标系转换，或通过迭代最近点算法(ICP：IterativeClosestPoint)。

7.根据权利要求1所述方法，其特征在于，所述步骤(5)中，其子块式表面跟踪算法包含步骤：

(1)将需要跟踪的软组织初始表面网格节点X＝{x,i＝1,2,……M}划分成若干小块：P＝{p,k＝1,2,……L}，除了小块p_L以外,每个小块p_k包含l个节点，小块p_L包含M-(L-1)*l个节点；

(2)为每个点集小块p_k构建能量最小方程获得点其表面位移；

其中p_k＝{p,o＝1,2,……l}为点集小块，Y＝{y,j＝1,2,……N}为点云的点集，c_jk是定义两个点集之间对应概率的相关性矩阵，λ为权重因子，用薄板样条(Thin-plateSpline)非刚性配准算法计算变形皮层(Y)与未变形皮层(p_k)之间映射函数f，获得点集小块表面位移Disp_pk；

(3)集合每个点集小块的位移Disp_pk，获得所有需要跟踪的软组织表面网格节点的位移Disp_x。

8.根据权利要求5所述方法，其特征在于，所述三维成设备是术中三维B超，或是术中磁共振或术中三维计算机断层扫描设备。

9.根据权利要求6所述方法，其特征在于，使用不同空间的坐标系转换实现初始软组织表面和变形后软组织表面的初配准包括如下步骤：

(1)将三维激光扫描仪空间变换到跟踪工具空间；

(2)将跟踪工具空间变换到跟踪设备空间；

(3)将跟踪设备空间变换到参考架空间；

(4)将参考架空间变换到图像空间，图像空间为变形前软组织所在空间。

10.根据权利要求10所述方法，其特征在于，所述的跟踪工具固定在三维激光扫描仪上，跟踪工具被跟踪设备跟踪；参考架固定在患者身体上，参考架被跟踪设备跟踪。