| WO/2012/009560 | DIGITAL X-RAY FIELD AND LIGHT FIELD ALIGNMENT |
| WO/2003/093869 | X-RAY EXAMINATION APPARATUS INCLUDING A DOSIMETER |
| JP10234716 | X-RAY DEVICE |
宋志坚 (中国上海市医学院路138号251信箱, Shanghai 2, 200032, CN)
WANG, Manning (P. O. Box 251, No. 138 Yixueyuan Road, Shanghai 2, 200032, CN)
王满宁 (中国上海市医学院路138号251信箱, Shanghai 2, 200032, CN)
ZHANG, Chenxi (P. O. Box 251, No. 138 Yixueyuan Road, Shanghai 2, 200032, CN)
复旦大学 (中国上海市邯郸路220号, Shanghai 3, 200433, CN)
SONG, Zhijian (P. O. Box 251, No. 138 Yixueyuan Road, Shanghai 2, 200032, CN)
宋志坚 (中国上海市医学院路138号251信箱, Shanghai 2, 200032, CN)
WANG, Manning (P. O. Box 251, No. 138 Yixueyuan Road, Shanghai 2, 200032, CN)
王满宁 (中国上海市医学院路138号251信箱, Shanghai 2, 200032, CN)
| 权利要求 1、 一种模拟脑组织图像变形及矫正的方法， 基于线弹性理论的物理模型， 利用带纹理成像的三维激光扫描仪获取到开颅后被暴露脑组织皮层表面点云数 据， 用于驱动已建立的模型进行计算， 其特征包括以下步骤： (1)采用基于 Level set 的自动分割方法， 把脑组织从颅骨中分离出来； (2)将分割后的结果进行多分辨率网格化处理， 划分出的四面体网格靠边界 处密度较大， 而内部则密度小； (3)对每一网格单元赋予相应的生物力学属性， 建立脑组织的线弹性模型； (4)将带纹理成像的三维激光扫描仪对准开颅区域进行扫描， 获取变形后皮 层表面点云数据，通过表面跟踪算法计算表面节点位移， 作为有限元计算的边界 条件； (5)根据已获取的边界条件， 用有限元的方法计算整个脑组织的变形场； (6)采用反向插值法来更新术前图像， 并用表面投影算法进行三维可视化。 2、 根据权利要求 1所述的一种模拟脑组织图像变形及矫正的方法， 其特征 是， 所述步骤 （1 ) 中基于 Level set的自动分割方法的主要步骤是： 先进行粗 略的分割， 设定门限值后， 再进行腐蚀、 区域生长、 膨胀、 再次区域生长四步操 作； 然后用 Level set算法进行细化处理。 3、 根据权利要求 1所述的一种模拟脑组织图像变形及矫正的方法， 其特征 是， 所述步骤 （2 ) 中多分辨率网格算法是先使用八叉树算法把三维图像四面体 网格化， 再把靠近边界的四面体细化， 最后采用移动四面体算法（Marching Tetrahedron)把网格从背景中切割出来， 并采用 case 表来提高网格划分的效 率。 4、 根据权利要求 1所述的一种模拟脑组织图像变形及矫正的方法， 其特征 是， 所述步骤 （3 ) 中表示线弹性模型的偏微分方程为， 1 d 3μχ + duy + dik、 t Fy l - 2v dy dx dy dz 其中， "为位移矢量， F为力 矢量， 为泊松比， =^ ？ Γ7 Α ^为弹性模量。 5、 根据权利要求 1所述的一种模拟脑组织图像变形及矫正的方法， 其特征 是， 所述步骤 （4 ) 通过构建能量最小方程实现开颅区域脑皮层的表面跟踪： min E(C, /) = min∑∑ c, || - f (x, )|f + 其中 = (x, i=l, 2, ····· -M)为表面网格的点集， = (y, j=l， 2， ····· -N)为 点云的点集， 是定义两个点集之间对应概率的相关性矩阵， J 为权重因子， 用 TPS (Thin-plate Spl ine)非刚性配准算法计算变形皮层（ 与未变形皮层 U) 之间映射函数 /， 从而获得表面位移《。。 6、 根据权利要求 1所述的一种模拟脑组织图像变形及矫正的方法， 其特征 是， 所述步骤 （5 ) 使用有限元的方法将模型离散化， 求解出这一系列偏微分方 程， 根据已获得的边界条件， 计算整个脑组织的变形场， 通过使用基于 PETSc2. 3. 2 (Portable, Extens ible Toolkit for Scient ific Computat ion)构 建的线性求解器来实现多处理器并行计算， 用以加快运算速度。 7、 根据权利要求 1所述的一种模拟脑组织图像变形及矫正的方法， 其特征 是， 所述步骤 （6 ) 中反相插值算法为， 从变形后的网格单元出发， 寻找出单元 内的整数坐标点，利用形函数获得该点在未变形前的位置， 再利用三线性插值获 得该点的灰度值， 然后更新术前图像的三维数据场。 |
[ 1] M. Knauth, C. R. Wirtza, V. M. Tronniera, N. Arasa, S. Kunzea, and
K. Sartora, ^Intraoperative MR imaging increases the extent of tumor resection in patients with high-grade gl iomas,〃 American journal of neuroradiology vol. 20， pp. 1642 - 1646 1999.
[2] F. A. Jolesz, ^Image-guided procedures and the operating room of the future Radiology, vol. 204， pp. 601-602， 1997.
[3] 杜固宏， 周良辅， and 毛颖， "神经导航术中脑移位的研究 " 中国微侵袭 神经外科杂志 vol. 7， pp. 65-68， 2007
[4] R. D. Bucholz, D. D. Yeh, J. Trobaugh, L. L. McDurmont, C. D. Sturm,
C. Baumann, J. M. Henderson, A. Levy, and P. Kessman, 〃The correction of stereotactic inaccuracy caused by brain shift using an Lecture Notes In Computer Science vol. 1225， pp. 459 - 466， 1997.
[5] A. Nabavi, P. M. Black, D. T. Gering, C. _F. Westin, V. Mehta, R. S.
Pergol izzi, M. Ferrant, S. K. Warf ield, N. Hata, R. B. Schwartz, W.
M. W. I I I， R. Kikinis, and F. A. Jolesz, ^Serial intraoperative magnetic resonance imaging of brain shift.，〃 Neurosurgery vol. 48， pp. 787-797， 2001.
[6] 0. Skrinjar, A. Nabavi, and J. Duncan, ^Model-driven brain shift compensation, " Medical image analysis vol. 6， pp. 361 - 373， 2002.
[7] M. I. Miga, K. D. Paulsen, F. E. Kennedy, J. Hoopes, A. Hartov, and
D. W. Roberts, "Initial In-Vivo Analysis of 3d Heterogeneous Brain
Computations for Model-Updated Image-Guided Neurosurgery, /r Lecture
Notes In Computer Science, vol. 1496， pp. 743 - 752， 1998.
[8] K. D. Paulsen, M. I. Miga, F. E. Kennedy, P. J. Hoopes, A. Hartov, and D. W. Roberts, 〃A computational model for tracking subsurface tissue deformation during stereotactic neurosurgery,〃 IEEE transactions on bio-medical engineering vol. 46， pp. 213-225， 1999.
[9] D. Terzopoulos, J. Piatt, A. Barr, and K. Fleischer, 〃Elastical ly deformable models, /r Computer Graphics, vol. 21， pp. 206-214， 1987.
[ 10] 0. M. Skrinjar, C. Studholme, A. Nabavi, and J. S. Duncan, 〃Steps toward a stereo-camera-guided biomechanical model for brain shift compensation, " Lecture Notes In Computer Science, vol. 2082， 2001.
[11] T. K. Sinha, B. M. Dawant, V. Duay, D. M. Cash, R. J. Wei l, R. C.
Thompson, K. D. Weaver, and M. I. Miga, 〃An atlas-based method to compensate for brain shift： Prel iminary results 〃 Medical Image Analysis, vol. 11， pp. 128-145， 2007.
[12] M. I. Miga, K. D. Paulsen, J. M. Lemery, S. D. Eisner, A. Hartov, F. E. Kennedy, and D. W. Roberts, "Model-updated image guidance： initial cl inical experiences with gravity-induced brain deformation,〃 IEEE transactions on medical imaging vol. 18， pp. 866-874， 1999.
[ 13] D. L. G. Hi l l, C. R. Maurer, M. Y. Wang, R. J. Maciunas, J. A. Barwise : and J. M. Fitzpatrick, ^Estimation of intraoperative brain surface movement 〃 Lecture Notes In Computer Science, vol. 1205， pp. 449 - 458 1997.
[ 14] M. A. Audette, K. Siddiqi, F. P. Ferrie, and T. M. Peters, 〃An integrated range-sensing, segmentation and registration framework for the characterization of intra-surgical brain deformations in image-guided surgery,〃 Computer Vision and Image Understanding, vol. 89， pp. 226 - 251， 2002.
发明内容 本发明的目的是提供一种模拟脑组织图像变形 及矫正的方法，是一种实施简 单、 实时性较好、精度较高的神经外科手术导航系 统中脑变形校正的方法， 使手 术导航在临床应用中更精确、 实用和方便。
本发明采用基于线弹性理论的物理模型，利用 带纹理成像的三维激光扫描仪 ( Laser Range Scanner, LRS ) 获取到开颅后被暴露脑组织皮层表面点云数据 ， 用于驱动已建立的模型进行计算， 最后将更新后的图像传给 IGNS开始实施术中 导航。
本发明的技术方案是： 首先采用基于 Level set的自动分割方法， 把脑组织 从颅骨中分离出来，将分割后的结果用进行多 分辨率网格化， 建立相应的线弹性 模型。借助三维激光扫描仪获取的术中点云数 据， 结合表面跟踪算法获得开颅后 被暴露脑皮层的位移， 并以此作为线弹性模型有限元方程的边界条件 ， 求解整个 脑组织的变形场。最后采用反向插值法来更新 术前图像， 并用光线投射算法进行 三维显示， 最后将更新后的图像传给 IGNS开始实施术中导航。
本发明方法通过下属步骤实现：
1、 采用基于 Level set 的自动分割方法， 把脑组织从 MRI颅部图像中分离 出来， 并建立相应的网格模型； 采用基于 Level set的自动分割方法来完成自动 分割， 其步骤是： 首先进行粗略的分割， 设定门限值后， 再进行腐蚀、区域生长、 膨胀、 再次区域生长四步操作； 然后用 Level set算法进行细化处理。
2、 将分割后的结果进行多分辨率网格化处理， 划分出的四面体网格靠边界 处密度较大， 而内部则密度小； 其中， 网格化是将连续三维空间离散化， 划分成 若干个通过节点相连的单元， 由于脑组织边界的变形比内部大， 为了精确模拟脑 组织的边界变形情况,我们对其进行多分辨率 格化处理， 划分出的四面体网格 靠边界处密度较大， 而内部则密度小。多分辨率网格化结合了八叉 树和移动四面 体两种算法， 并采用 case 表来提高计算效率。 具体算法如下：
(1) 采用类八叉树算法将三维图像空间划分为由均 匀的六面体单元所组成的 网格；
(2) 针对每一六面体单元将其划分为五个四面体单 元；
(3) 对于（2 ) 中所获得均匀四面体网格， 在边界处加以细化， 获得多分辨率 网格； (4) 采用移动四面体算法 (Marching Tetrahedron) 对网格进行切割，去除背 景，把网格从背景中切割出来， 并采用 case 表来提高网格划分的效率， 获得最终 的仅包含脑组织的多分辨率网格。
3、 对每一网格单元赋予相应的生物力学属性， 建立脑组织的线弹性模型， 所述的线弹性模型基于把大脑看成均匀的、各 向同性的弹性体， 用一系列偏微分 方程来来表示， 其中， 所述线弹性模型的偏微分方程表达如下：
其中， "为位移矢量， 为力矢量， ^为泊松比， H =E/2 (1+ v ), ^为弹性模 模型建立后， 需要用表面跟踪来计算对应表面节点的位移， 作为有限元计算 的边界条件来代表术中大脑变形过程中所受的 驱动力。
在表面跟踪前， 必须将术中点云和表面网格进行配准， 这首先需要完成以下 四个空间的坐标转换： （l) LRS空间与跟踪工具空间的坐标转换； （2) Polaris空 间与跟踪工具空间的坐标转换； （3) Polaris空间与参考架空间的坐标转换； （4) 参考架空间与图像空间的转换 (参考架固定在猪头部）。在实验过程中， 第一个坐 标转换在 BDS中的校准模块中完成； 第二、 三、 四个坐标转换关系在 IGNS中生 成，并通过通讯模块传给 BDS。获取四个坐标转换关系后，将初始皮层表 面（image place)和变形后皮层表面（LRS-space)配准到同一 标空间 （如图 4) ， 进而再 用 RPM ( Robust Point Matching) 算法对脑皮层的变化进行追踪。
4、 将带纹理成像的三维激光扫描仪对准开颅区域 进行扫描， 获取变形后皮 层表面点云数据，通过表面跟踪算法计算表面 节点位移， 作为有限元计算的边界 条件；
在上述方案基础上， 所述步骤 （4) 通过构建能量最小方程实现开颅区域脑 皮层的表面跟踪：
E(C,/) = min∑∑ C , ||) ; f
dx 其中 = (x, i=l，2， ····· ·Μ}为表面网格的点集， r= (y, j=l, 2, ····· -Nj为 点云的点集， ¾ 是定义两个点集之间对应概率的相关性矩 阵， J 为权重因子， 用 TPS (Thin-plate Spl ine)非刚性配准算法计算变形皮层（ 与未变形皮层 0) 之间映射函数 /， 从而获得表面位移 α。。
5、 根据已获取的边界条件， 用有限元的方法计算整个脑组织的变形场， 使 用有限元的方法将模型离散化， 求解出这一系列偏微分方程， 根据已获得的边界 条件， 计算整个脑组织的变形场， 离散化之后的有限元方程为：
Κα = Ρ ί为刚性矩阵， α为结点位移， Ρ为结点载荷。 Κ由杨氏弹性模量 Ε和泊松 比 两个参数决定。
通过使用基于 PETSc2. 3. 2 (Portable, Extens ible Toolkit for Scient ific Computat ion)构建的线性求解器来实现多处理器并行计 ， 用以加快运算速 度。
6、 经过有限元计算， 采用反向插值法来更新术前图像， 并用光线投射算法 进行显示， 即表面投影算法进行三维可视化， 所述的反相插值算法为， 从变形后 的网格单元出发， 寻找出单元内的整数坐标点， 利用形函数获得该点在未变形前 的位置，再利用三线性插值获得该点的灰度值 ，然后更新术前图像的三维数据场。
本发明具有下述优点：
( 1)采用了基于 Level set的自动分割方法， 使颅部 MRI图像的脑组织分割 出的结果更加准确、 精细。
(2)提出了一种基于 TPS的表面跟踪算法， 确保了表面不光滑的扫描点云与 术前图像之间特征点的一一对应的双向约束， 具有较强的有效性和鲁棒性。
(3)在有限元计算中使用了多处理器并行计算， 具有更高的计算效率。 有利 于图像的快速更新， 增强了导航的实时性。
附图说明 图 1为脑变形校正的流程图；
图 2为基于 256 X 256 X 48MRI数据场的基于 Level set的自动分割算法结果； 图 3为针对 256 X 256 X 48MRI数据场的多分辨率网格化结果（a.断层图 , b. 断层网格, c.三维网格）；
图 4 为初始皮层表面与变形后皮层表面配准到同一 坐标空间；
图 5为表面跟踪的结果。 其中左图为初始表面低于变形后的实际脑皮层 表 面， 右图为是变形前跟踪结束后网格初始表面和点 云表面重叠到一起；
图 6为脑变形模型预测结果的叠加显示（灰色组 为未变形的脑组织， 红色 网格为模型预测的变形结果） ；
图 7为位移场的三维可视化结果， 箭头表示位移方向， 颜色表示位移大小， 颜色由红到蓝表示位移由小到大。
具体实施方式 实施例 1
1. 针对 256 X 256 X 48MRI数据场采用基于 Level set 的自动分割算法， 将脑组 织从 MRI 颅部图像中分离出来；
2. 采用多分辨率网格化的方法， 建立术前的网格模型。 将所分割出的脑组织离 散为 18485个四面体， 节点个数为 5410。 边界处最大四面体 7. 5 X 7. 5 X
7. 5 mm 3 ;
3. 在线弹性模型的每一个单元设置脑组织生物力 学属性参数。杨氏模量 =3Kpa 泊松比 =0. 45
4. 采用坐标变化实现刚体配准， 获得脑皮层在 LRS空间的初始位置。 启动 Polaris 对参考架和跟踪工具进行跟踪， 获取 Polaris 空间与跟踪工具空 间以及参考架空间之间的坐标转换关系。
5. 在开颅后， 手术前利用 LRS 进行一次扫描， LRS距离病人开颅位置的距离控 制在 450 左右， 扫描时间 5 7s， 分辨率 512 (水平） X 500 (垂直） 根 据设备提供的变形后的三维表面信息以及由步 骤 4获得的变形前三维表面信 息， 应用表面跟踪算法， 获得表面节点位移， 作为有限元求解的边界条件。
6. 采用有限元方法将线弹性模型方程化为矩阵形 式： Au =b (A 为总刚矩阵， u 为待求位移矢量， b为节点力矢量） 采用共轭梯度法求解含有 5410 (节 点数） X 3 (维度） =16230个方程的线性方程组。 引入边界条件， 消除矩阵 A的奇异性， 求出位移矢量 u， 获得脑组织内部任意节点处的位移， 再结合 有限元中的形函数就可以插值出脑组织内部任 意位置的变形。
由变形后的数据场， 根据形函数计算出所有对显示有贡献的坐标点 （整数坐 标点）变形前的位置， 利用三线性插值计算出该点的灰度值。 最后对所有这 些整数坐标点组成的三维数据场采用光线透射 法加以可视化， 用于导航手 术。
在奔腾 IV2. 6G、 1G内存、 Win2000操作系统下运行， 需要 2分钟左右。