一种基于多gpu的人脑变形仿真方法
【技术领域】
[0001] 本发明涉及一种基于多GPU的人脑变形仿真方法，属于计算机虚拟现实领域。
【背景技术】
[0002] 神经外科的导航系统出现于上世纪90年代，通过把现代神经影像（MRI、CT等医 学图像）、计算机三维图像处理、脑立体定位与显微神经外科结合，来提高神经外科手术的 精确性和安全性。但是对于常规导航系统来说脑位移（brainshift)始终是一个难以解决 的技术难题。在实际手术中由于脑组织生物力学属性、重力作用、颅内压变化及手术操作等 影响，常发生脑变形。术中脑变形会导致导航定位精度的下降，超出医学可接受范围后会 严重影响定位系统可信度。因此纠正脑位移误差技术的研究已成为该领域的热点问题。目 前主要途径为术中成像技术（intraoperativeimaging)和模型校正技术（model-updated imaging)。其中模型校正技术仅需要使用现有医学图像和电脑的虚拟模拟技术相结合，因 此已成为国际上脑位移纠正技术的研究热点。
[0003] 研究之初，人们多使用几何模型来预测脑部发生的位移。但是经过数年的实验人 们发现几何模型仅仅能对大脑表面的变形进行较好的拟合，但是无法真实反映大脑内部的 变化。因此，近年来研究人员开始研究使用生物力学模型去预测脑部开颅后的变化。基于生 物力学模型的校正技术在以往的研究中由于硬件性能和算法设计难以实现高度精确的实 时交互的脑漂移预测。为满足医用系统实时性的要求，研究人员往往选择牺牲模型精度，而 使用线弹性模型或固结模型来模拟脑部变形，其结果使得系统的使用受到诸多条件限制， 同时系统的可靠性也有一定程度的损失。
[0004] 相对于线性模型，非线性模型更加符合真实的脑组织特征，通过实验获得合适的 参数后，使用非线性模型计算得到的结果比线性模型误差更小，同时限制条件也相对更加 宽松。但是相对于线性模型，非线性模型的计算复杂度更大，在引入并行计算之前，无法满 足实时交互的要求。到2014年为止，非线性算法的研究实现仅能保证在不到一万个有限元 单元组成的模型下进行实时模拟。计算速度过慢，成为了非线性模型无法投入实际应用的 瓶颈。
【发明内容】
[0005] 本发明要解决的技术问题是：克服现有技术的不足，提供一种基于多GPU的人脑 变形仿真方法，可以在计算机构造的虚拟环境中逼真地模拟脑部的变形过程，满足了人脑 变形仿真中的精度与速度要求。
[0006] 1、本发明技术解决方案：一种基于多GPU的人脑变形仿真方法，其特征在于包括 如下步骤：
[0007] (1)初始化，读入人脑四面体网格模型数据，分析四面体网格单元的拓扑结构关系 得到四面体网格邻接体单元的集合，依据所获得的集合信息和所用GPU核数目，对读入的 四面体网格模型原始数据进行分组、编号，然后基于等参元法计算各四面体单元的形函数 矩阵，用于后面在发生受力或表面位移后内部单元的形变计算；
[0008] (2)初始化完成后，根据手动输入或是力反馈设备捕获的边界受力变化或是扫描 仪或相机捕获的可见表面节点的位移变化，利用Neo-Hookean模型，该模型定义了力与顶 点位移之间的关系，总体拉格朗日方法，计算每一个时间步的组成四面体单元的各顶点的 受力，在算法的实现过程中在多GPU服务器上并行实现，通过这种处理方式来减少多个GPU 之间的数据交互；
[0009] (3)把步骤二中每个GPU核计算出的受力结果进行合并，并根据最终计算结果采 用中心差分法公式
，其中At表示时间步，Μ表示 四面体单元的质量，un+1表示下一个时间步的位移，Rn表示当前时刻的反作用力，Fn表示当 前时刻所施加的外力即步骤二捕获的表面受力，unl表示前一时刻的位移，un表示当前时刻 的位移即步骤二中捕获的表面节点的位移变化，更新四面体单元各顶点的坐标变化，并在 多CPU上实现。
[0010] 步骤（1)中的对读入的四面体网格模型原始数据进行分组、编号，然后基于等参 元法计算各四面体单元的形函数矩阵，包括步骤如下：
[0011] (2. 1)遍历所有四面体网格单元的顶点，为每一个顶点建立所属单元列表，为每个 顶点建立访问标志位并设置为false;
[0012] (2. 2)输入分组起始参考顶点，对四面体网格模型单元进行分组，分组单元数=四 面体总单元数/所用GPU核数，首先将起始点的所属单元加入分组，并将顶点访问标志设为 true，判断当前分组内包含四面体是否达到分组单元数限制，达到则从下一个输入的四面 体单元顶点开始，将四面体放入下一分组，否则读取参考点所属单元列表内访问标志顶点 为false的节点的所属单元并加入当前分组中；
[0013] (2. 3)步骤（2. 2)对四面体单元进行分组完成后，将各组内顶点根据所属单元列 表是否完全在当前分组内进行区分，满足条件的为组内点，否则为边界点；对所有分组内 四面体单元的顶点依照分区次序从低到高，对这些四面体单元组内顶点和边界点进行重编 号；
[0014] (2. 4)根据新的编号，更新所有的单元顶点编号信息；
[0015] (2. 5)根据各单元顶点信息采用等参元法计算得到各单元和标准单元形函数的变 换关系矩阵，等参元法是将各个四面体单元与正四面体的形函数之间建立映射关系，表现 形式为一个3*3矩阵，在每一个时间步中计算顶点坐标发生变化后，通过该变换关系获取 新的形函数矩阵。
[0016] 步骤（2)中采用多GPU加速计算过程，进行数据分组减少多GPU数据交互，包括步 骤如下：
[0017] (3. 1)初始化：初始化完成后，将四面体单元的顶点的坐标信息和单元节点编号 信息存储于同一数据结构中，保证在单元信息传送到GPU端后，已传送单元相关的顶点信 息也已经传送完毕；
[0018] (3. 2)流式传输：为了充分利用多GPU服务器上每个GPU核的计算能力及数据传 输带宽，把输入数据分成若干组，每次传输一组数据后，在开启GPU计算的同时进行下一组 数据的传输；
[0019] (3· 3)在GPU端首先使用由Neo-Hookean模型及总体拉格朗日方法推导出的第二 Piola-Kirchoff应力公5
1其中S表示第二 Piola-Kirchoff应力，μ表示剪切模量，J表示变形梯度，I是三阶单位矩阵，^是右柯西 格林张量第一不变量，C是右柯西格林张量，δ是材料的体积弹性模量，通过该公式计算各 单元的受力情况；
[0020] (3. 4)通过初始化计算时生成的各四面体单元顶点受力求和求整个四面体网格模 型各顶点的总受力，并将结果传送回CPU端。
[0021] 所述步骤（2)中总体拉格朗日方法实现过程为：方法在预测模型变形时通过把整 个连续的变形过程离散为若干个时刻的变形，设计算开始的时间为t，计算变形的下一个时 亥|J和时间t的时间间隔叫做时间步，记做At，首先使用顶点的初始坐标、步骤⑴的形函数 矩阵及捕获的表面位移或受力计算时刻t+△t时每个四面体单元四个顶点的受力情况，然 后通过步骤（3)的计算获得t+At时刻顶点坐标、形函数矩阵及捕获的表面位移或受力计 算时刻t+2* △t时刻t+ △t时每个四面体单元四个顶点的受力情况以此类推知道两个时间 步之间顶点位移不再发生变化。
[0022] 本发明与现有技术相比的优点在于：
[0023] (1)本发明使用总体拉格朗日法、中心差分法及非线性有限元Neo-Hookean模型 来进行脑部变形的显示迭代仿真，通过在GPU端的并行化算法实现加速仿真计算的速度， 在引入多GPU进一步提高并行能力之后，对原始数据采用广度优先搜索方法进行数据分组 和重新编号最小化各结点之间的数据相关性，并通过使用额外的数据结构使得单个节点内 的输入数据及计算过程存在一一对应关系，进而可以使用流式传输实现计算和数据传输的 并行化，最终利用混合CPU/GPU的多核计算架构对整个算法进行加速，满足了人脑变形仿 真中的精度与速度要求。
[0024] (2)本发明充分利用了人脑变形中的细微性以及时空连续性特点，可以在计算机 构造的虚拟环境中逼真地模拟脑部的变形过程。
【附图说明】
[0025] 图1为本发明基于多GPU的人脑变形仿真方法的流程图；
[0026] 图2为本方法中的初始化阶段数据分组，编号更新及准备数据流程图；
[0027] 图3为变形前后模型对比图。
[0028] 具体实施方法
[0029] 如图1、2所示，本发明的具体步骤如下：
[0030] (1)读入人脑模型数据，使用广度搜索分析拓扑结构关系得到邻接体单元集合，并 保存在哈希表中。初始化单元访问位和结点访问位列表为false。
[0031] (2)选取分组起始参考结点，设置为当前访问结点，并将该结点访问标识设置为 true。设置当前分组为第一分组。
[