CN105389853A - 一种基于多gpu的人脑变形仿真方法 - Google Patents

Abstract

本发明涉及一种基于多GPU的人脑变形仿真方法，使用总体拉格朗日法、中心差分法及非线性有限元Neo-Hookean模型来进行脑部变形的显式迭代仿真，通过在GPU端的并行化算法实现加速仿真计算的速度，在引入多GPU进一步提高并行能力之后，对原始数据采用广度优先搜索方法进行数据分组和重新编号最小化各结点之间的数据相关性，并通过使用额外的数据结构使得单个节点内的输入数据及计算过程存在一一对应关系，进而可以使用流式传输实现计算和数据传输的并行化，最终利用混合CPU/GPU的多核计算架构进行加速，满足了人脑变形仿真中的精度与速度要求。本发明充分利用了人脑变形中的细微性以及时空连续性特点，可以在计算机构造的虚拟环境中逼真地模拟脑部的变形过程。

Description

一种基于多GPU的人脑变形仿真方法

技术领域

本发明涉及一种基于多GPU的人脑变形仿真方法，属于计算机虚拟现实领域。

背景技术

神经外科的导航系统出现于上世纪90年代，通过把现代神经影像(MRI、CT等医学图像)、计算机三维图像处理、脑立体定位与显微神经外科结合，来提高神经外科手术的精确性和安全性。但是对于常规导航系统来说脑位移(brainshift)始终是一个难以解决的技术难题。在实际手术中由于脑组织生物力学属性、重力作用、颅内压变化及手术操作等影响，常发生脑变形。术中脑变形会导致导航定位精度的下降，超出医学可接受范围后会严重影响定位系统可信度。因此纠正脑位移误差技术的研究已成为该领域的热点问题。目前主要途径为术中成像技术(intraoperativeimaging)和模型校正技术(model-updatedimaging)。其中模型校正技术仅需要使用现有医学图像和电脑的虚拟模拟技术相结合，因此已成为国际上脑位移纠正技术的研究热点。

研究之初，人们多使用几何模型来预测脑部发生的位移。但是经过数年的实验人们发现几何模型仅仅能对大脑表面的变形进行较好的拟合，但是无法真实反映大脑内部的变化。因此，近年来研究人员开始研究使用生物力学模型去预测脑部开颅后的变化。基于生物力学模型的校正技术在以往的研究中由于硬件性能和算法设计难以实现高度精确的实时交互的脑漂移预测。为满足医用系统实时性的要求，研究人员往往选择牺牲模型精度，而使用线弹性模型或固结模型来模拟脑部变形，其结果使得系统的使用受到诸多条件限制，同时系统的可靠性也有一定程度的损失。

相对于线性模型，非线性模型更加符合真实的脑组织特征，通过实验获得合适的参数后，使用非线性模型计算得到的结果比线性模型误差更小，同时限制条件也相对更加宽松。但是相对于线性模型，非线性模型的计算复杂度更大，在引入并行计算之前，无法满足实时交互的要求。到2014年为止，非线性算法的研究实现仅能保证在不到一万个有限元单元组成的模型下进行实时模拟。计算速度过慢，成为了非线性模型无法投入实际应用的瓶颈。

发明内容

本发明要解决的技术问题是：克服现有技术的不足，提供一种基于多GPU的人脑变形仿真方法，可以在计算机构造的虚拟环境中逼真地模拟脑部的变形过程，满足了人脑变形仿真中的精度与速度要求。

1、本发明技术解决方案：一种基于多GPU的人脑变形仿真方法，其特征在于包括如下步骤：

(1)初始化，读入人脑四面体网格模型数据，分析四面体网格单元的拓扑结构关系得到四面体网格邻接体单元的集合，依据所获得的集合信息和所用GPU核数目，对读入的四面体网格模型原始数据进行分组、编号，然后基于等参元法计算各四面体单元的形函数矩阵，用于后面在发生受力或表面位移后内部单元的形变计算；

(2)初始化完成后，根据手动输入或是力反馈设备捕获的边界受力变化或是扫描仪或相机捕获的可见表面节点的位移变化，利用Neo-Hookean模型，该模型定义了力与顶点位移之间的关系，总体拉格朗日方法，计算每一个时间步的组成四面体单元的各顶点的受力，在算法的实现过程中在多GPU服务器上并行实现，通过这种处理方式来减少多个GPU之间的数据交互；

(3)把步骤二中每个GPU核计算出的受力结果进行合并，并根据最终计算结果采用中心差分法公式其中Δt表示时间步，M表示四面体单元的质量，u_n+1表示下一个时间步的位移，R_n表示当前时刻的反作用力，F_n表示当前时刻所施加的外力即步骤二捕获的表面受力，u_n-1表示前一时刻的位移，u_n表示当前时刻的位移即步骤二中捕获的表面节点的位移变化，更新四面体单元各顶点的坐标变化，并在多CPU上实现。

步骤(1)中的对读入的四面体网格模型原始数据进行分组、编号，然后基于等参元法计算各四面体单元的形函数矩阵，包括步骤如下：

(2.1)遍历所有四面体网格单元的顶点，为每一个顶点建立所属单元列表，为每个顶点建立访问标志位并设置为false；

(2.2)输入分组起始参考顶点，对四面体网格模型单元进行分组，分组单元数＝四面体总单元数/所用GPU核数，首先将起始点的所属单元加入分组，并将顶点访问标志设为true，判断当前分组内包含四面体是否达到分组单元数限制，达到则从下一个输入的四面体单元顶点开始，将四面体放入下一分组，否则读取参考点所属单元列表内访问标志顶点为false的节点的所属单元并加入当前分组中；

(2.3)步骤(2.2)对四面体单元进行分组完成后，将各组内顶点根据所属单元列表是否完全在当前分组内进行区分，满足条件的为组内点，否则为边界点；对所有分组内四面体单元的顶点依照分区次序从低到高，对这些四面体单元组内顶点和边界点进行重编号；

(2.4)根据新的编号，更新所有的单元顶点编号信息；

(2.5)根据各单元顶点信息采用等参元法计算得到各单元和标准单元形函数的变换关系矩阵，等参元法是将各个四面体单元与正四面体的形函数之间建立映射关系，表现形式为一个3*3矩阵，在每一个时间步中计算顶点坐标发生变化后，通过该变换关系获取新的形函数矩阵。

步骤(2)中采用多GPU加速计算过程，进行数据分组减少多GPU数据交互，包括步骤如下：

(3.1)初始化：初始化完成后，将四面体单元的顶点的坐标信息和单元节点编号信息存储于同一数据结构中，保证在单元信息传送到GPU端后，已传送单元相关的顶点信息也已经传送完毕；

(3.2)流式传输：为了充分利用多GPU服务器上每个GPU核的计算能力及数据传输带宽，把输入数据分成若干组，每次传输一组数据后，在开启GPU计算的同时进行下一组数据的传输；

(3.3)在GPU端首先使用由Neo-Hookean模型及总体拉格朗日方法推导出的第二Piola–Kirchoff应力公式 $S_0^t=\mathrm{\μ}\×J^{-2/3}(I-\frac{1}{3}{I}_1\×{C_0^t}^{-1})+\mathrm{\δ}(J-1)\×J\×{C_0^t}^{-1}$ 其中S表示第二Piola–Kirchoff应力，μ表示剪切模量，J表示变形梯度，I是三阶单位矩阵，I₁是右柯西格林张量第一不变量，C是右柯西格林张量，δ是材料的体积弹性模量，通过该公式计算各单元的受力情况；

(3.4)通过初始化计算时生成的各四面体单元顶点受力求和求整个四面体网格模型各顶点的总受力，并将结果传送回CPU端。

所述步骤(2)中总体拉格朗日方法实现过程为：方法在预测模型变形时通过把整个连续的变形过程离散为若干个时刻的变形，设计算开始的时间为t，计算变形的下一个时刻和时间t的时间间隔叫做时间步，记做Δt，首先使用顶点的初始坐标、步骤(1)的形函数矩阵及捕获的表面位移或受力计算时刻t+Δt时每个四面体单元四个顶点的受力情况，然后通过步骤(3)的计算获得t+Δt时刻顶点坐标、形函数矩阵及捕获的表面位移或受力计算时刻t+2*Δt时刻t+Δt时每个四面体单元四个顶点的受力情况以此类推知道两个时间步之间顶点位移不再发生变化。

本发明与现有技术相比的优点在于：

(1)本发明使用总体拉格朗日法、中心差分法及非线性有限元Neo-Hookean模型来进行脑部变形的显示迭代仿真，通过在GPU端的并行化算法实现加速仿真计算的速度，在引入多GPU进一步提高并行能力之后，对原始数据采用广度优先搜索方法进行数据分组和重新编号最小化各结点之间的数据相关性，并通过使用额外的数据结构使得单个节点内的输入数据及计算过程存在一一对应关系，进而可以使用流式传输实现计算和数据传输的并行化，最终利用混合CPU/GPU的多核计算架构对整个算法进行加速，满足了人脑变形仿真中的精度与速度要求。

(2)本发明充分利用了人脑变形中的细微性以及时空连续性特点，可以在计算机构造的虚拟环境中逼真地模拟脑部的变形过程。

附图说明

图1为本发明基于多GPU的人脑变形仿真方法的流程图；

图2为本方法中的初始化阶段数据分组，编号更新及准备数据流程图；

图3为变形前后模型对比图。

具体实施方法

如图1、2所示，本发明的具体步骤如下：

(1)读入人脑模型数据，使用广度搜索分析拓扑结构关系得到邻接体单元集合，并保存在哈希表中。初始化单元访问位和结点访问位列表为false。

(2)选取分组起始参考结点，设置为当前访问结点，并将该结点访问标识设置为true。设置当前分组为第一分组。

(3)从存储邻接体单元的集合的哈希表中提取当前结点包含的体单元，将所有单元访问位为true的单元压入当前单元分组，并设置这些单元的单元访问位为true，然后将所包含单元其它结点中结点访问位为false的结点依次压入结点队列。

(4)判断所有节点是否完成遍历，满足条件转步骤6，否则判断单元数是否达到分组单元数上限，条件满足则将当前分组设置为下一分组。

(5)从队列中提取第一个结点编号，并将该结点访问标识设置为true，转步骤3。

(6)分组完成后，将各组内结点根据所属单元列表是否完全在当前分组内进行分组，满足条件的为组内点，否则为边界点。对所有结点依照分区及优先组内点，然后编号边界点的顺序进行重编号。

(7)根据新的编号，更新所有的单元结点编号信息。并将结点坐标信息和单元结点编号信息存储在一起。

(8)根据各单元顶点信息采用等参元法计算得到各单元和标准单元形函数的变换关系矩阵，各结点质量矩阵及单元受力索引。其中等参元法是将各个四面体单元与正四面体的形函数之间建立映射关系，表现形式为一个3*3矩阵，在每一个时间步中计算顶点坐标发生变化后，通过该变换关系获取新的形函数矩阵。同时输入模型参数杨氏模量、泊松比及时间步。

(10)GPU端初始化：将除结点坐标信息外的其它信息传送到GPU端。

(11)根据手动输入或是力反馈设备捕获的边界受力变化或是扫描仪或相机捕获的可见表面节点的位移变化，利用Neo-Hookean模型，该模型定义了力与顶点位移之间的关系，总体拉格朗日方法，计算每一个时间步的组成四面体单元的各顶点的受力。在多GPU的并行实现中将四面体单元的顶点坐标信息分成若干段依次传输给GPU，同时在每段数据传输完成后开启GPU端该段数据的计算函数 $S_0^t=\mathrm{\μ}\×J^{-2/3}(I-\frac{1}{3}{I}_1\×{C_0^t}^{-1})+\mathrm{\δ}(J-1)\×J\×{C_0^t}^{-1},$ 其中S表示第二Piola–Kirchoff应力，μ表示剪切模量，J表示变形梯度，I是三阶单位矩阵，I₁是右柯西格林张量第一不变量，C是右柯西格林张量，δ是材料的体积弹性模量，通过该公式计算各单元结点的受力情况。其中分段大小依据GPU性能动态设定。这种分段处理方式可以减少多个GPU之间的数据交互。

(12)各单元结点受力计算完成后，通过预存的单元受力索引计算结点所受合力。

(13)将受力信息传送回CPU端之后，调用碰撞检测和响应程序进行处理后，采用中心差分法公式各节点的坐标变化计算出下一步的结点位移。其中Δt表示时间步，M表示四面体单元的质量，u_n+1表示下一个时间步的位移，R_n表示当前时刻的反作用力，F_n表示当前时刻所施加的外力，u_n-1表示前一时刻的位移，u_n表示当前时刻的位移即步骤二中捕获的表面节点的位移变化，更新四面体单元各顶点的坐标变化。判断模型是否已经收敛稳定。条件满足则中止计算。否则转步骤11继续计算模型变化。

如图3所示，左侧为变形前的模型，右侧为变形后的结果。

本发明中未详细阐述的部分属于本领域技术人员的公知技术。

Claims (4)

1.一种基于多GPU的人脑变形仿真方法，其特征在于包括如下步骤：

(1)初始化，读入人脑四面体网格模型数据，分析四面体网格单元的拓扑结构关系得到四面体网格邻接体单元的集合，依据所获得的集合信息和所用GPU核数目，对读入的四面体网格模型原始数据进行分组、编号，然后基于等参元法计算各四面体单元的形函数矩阵，用于后面在发生受力或表面位移后内部单元的形变计算；

(2)初始化完成后，根据手动输入或是力反馈设备捕获的边界受力变化或是扫描仪或相机捕获的可见表面节点的位移变化，利用Neo-Hookean模型，该模型定义了力与顶点位移之间的关系，总体拉格朗日方法，计算每一个时间步的组成四面体单元的各顶点的受力，在实现过程中采用多GPU加速计算过程，进行数据分组减少多GPU数据交互；

(3)把步骤二中每个GPU核计算出的受力结果进行合并，并根据最终计算结果采用中心差分法公式 $\left(\frac{1}{{\mathrm{\Δt}}^2}M\right)u_{n+1}=R_n-\underset{i}{\Σ}{F}_n^i-\frac{1}{{\mathrm{\Δt}}^2}M(u_{n-1}-2u_n),$ 其中Δt表示时间步，M表示四面体单元的质量，u_n+1表示下一个时间步的位移，R_n表示当前时刻的反作用力，F_n表示当前时刻所施加的外力即步骤二捕获的表面受力，u_n-1表示前一时刻的位移，u_n表示当前时刻的位移即步骤二中捕获的表面节点的位移变化，更新四面体单元各顶点的坐标变化，并在多CPU上实现。

2.根据权利要求1所述的基于多GPU的人脑变形仿真方法，其特征在于：步骤(1)中的对读入的四面体网格模型原始数据进行分组、编号，然后基于等参元法计算各四面体单元的形函数矩阵，包括步骤如下：

(2.1)遍历所有四面体网格单元的顶点，为每一个顶点建立所属单元列表，为每个顶点建立访问标志位并设置为false；

(2.2)输入分组起始参考顶点，对四面体网格模型单元进行分组，分组单元数＝四面体总单元数/所用GPU核数，首先将起始点的所属单元加入分组，并将顶点访问标志设为true，判断当前分组内包含四面体是否达到分组单元数限制，达到则从下一个输入的四面体单元顶点开始，将四面体放入下一分组，否则读取参考点所属单元列表内访问标志顶点为false的节点的所属单元并加入当前分组中；

(2.3)步骤(2.2)对四面体单元进行分组完成后，将各组内顶点根据所属单元列表是否完全在当前分组内进行区分，满足条件的为组内点，否则为边界点；对所有分组内四面体单元的顶点依照分区次序从低到高，对这些四面体单元组内顶点和边界点进行重编号；

(2.4)根据新的编号，更新所有的单元顶点编号信息；

(2.5)根据各单元顶点信息采用等参元法计算得到各单元和标准单元形函数的变换关系矩阵，等参元法是将各个四面体单元与正四面体的形函数之间建立映射关系，表现形式为一个3*3矩阵，在每一个时间步中计算顶点坐标发生变化后，通过该变换关系获取新的形函数矩阵。

3.根据权利要求1所述的多GPU的人脑变形仿真方法，其特征在于：步骤(2)中采用多GPU加速计算过程，进行数据分组减少多GPU数据交互，包括步骤如下：

(3.1)初始化：初始化完成后，将四面体单元的顶点的坐标信息和单元节点编号信息存储于同一数据结构中，保证在单元信息传送到GPU端后，已传送四面体单元相关的顶点信息也已经传送完毕；

(3.2)流式传输：为了充分利用多GPU服务器上每个GPU核的计算能力及数据传输带宽，把输入数据分成若干组，每次传输一组数据后，在开启GPU计算的同时进行下一组数据的传输；

(3.3)在GPU端首先使用由Neo-Hookean模型及总体拉格朗日方法推导出的第二Piola–Kirchoff应力公式 $S_0^t=\mathrm{\μ}\×J^{-2/3}(I-\frac{1}{3}{I}_1\×{C_0^t}^{-1})+\mathrm{\δ}(J-1)\×J\×{C_0^t}^{-1}$ 其中S表示第二Piola–Kirchoff应力，μ表示剪切模量，J表示变形梯度，I是三阶单位矩阵，I₁是右柯西格林张量第一不变量，C是右柯西格林张量，δ是材料的体积弹性模量，通过该公式计算各单元的受力情况；

(3.4)通过初始化计算时生成的各四面体单元顶点受力求和求整个四面体网格模型各顶点的总受力，并将结果传送回CPU端。

4.根据权利要求1所述的多GPU的人脑变形仿真方法，其特征在于：所述步骤(2)中总体拉格朗日方法实现过程为：方法在预测模型变形时通过把整个连续的变形过程离散为若干个时刻的变形，设计算开始的时间为t，计算变形的下一个时刻和时间t的时间间隔叫做时间步，记做Δt，首先使用顶点的初始坐标、步骤(1)的形函数矩阵及捕获的表面位移或受力计算时刻t+Δt时每个四面体单元四个顶点的受力情况，然后通过步骤(3)的计算获得t+Δt时刻顶点坐标、形函数矩阵及捕获的表面位移或受力计算时刻t+2*Δt时刻t+Δt时每个四面体单元四个顶点的受力情况以此类推知道两个时间步之间顶点位移不再发生变化。