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

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

【技术实现步骤摘要】

本专利技术涉及一种基于多GPU的人脑变形仿真方法，属于计算机虚拟现实领域。
技术介绍
神经外科的导航系统出现于上世纪90年代，通过把现代神经影像(MRI、CT等医学图像)、计算机三维图像处理、脑立体定位与显微神经外科结合，来提高神经外科手术的精确性和安全性。但是对于常规导航系统来说脑位移(brainshift)始终是一个难以解决的技术难题。在实际手术中由于脑组织生物力学属性、重力作用、颅内压变化及手术操作等影响，常发生脑变形。术中脑变形会导致导航定位精度的下降，超出医学可接受范围后会严重影响定位系统可信度。因此纠正脑位移误差技术的研究已成为该领域的热点问题。目前主要途径为术中成像技术(intraoperativeimaging)和模型校正技术(model-updatedimaging)。其中模型校正技术仅需要使用现有医学图像和电脑的虚拟模拟技术相结合，因此已成为国际上脑位移纠正技术的研究热点。研究之初，人们多使用几何模型来预测脑部发生的位移。但是经过数年的实验人们发现几何模型仅仅能对大脑表面的变形进行较好的拟合，但是无法真实反映大脑内部的变化。因此，近年来研究人员开始研究使用生物力学模型去预测脑部开颅后的变化。基于生物力学模型的校正技术在以往的研究中由于硬件性能和算法设计难以实现高度精确的实时交互的脑漂移预测。为满足医用系统实时性的要求，研究人员往往选择牺牲模型精度，而使用线弹性模型或固结模型来模拟脑部变形，其结果使得系统的使用受到诸多条件限制，同时系统的可靠性也有一定程度的损失。相对于线性模型，非线性模型更加符合真实的脑组织特征，通过实验获得合适的参数后，使用非线性模型计算得到的结果比线性模型误差更小，同时限制条件也相对更加宽松。但是相对于线性模型，非线性模型的计算复杂度更大，在引入并行计算之前，无法满足实时交互的要求。到2014年为止，非线性算法的研究实现仅能保证在不到一万个有限元单元组成的模型下进行实时模拟。计算速度过慢，成为了非线性模型无法投入实际应用的瓶颈。
技术实现思路
本专利技术要解决的技术问题是：克服现有技术的不足，提供一种基于多GPU的人脑变形仿真方法，可以在计算机构造的虚拟环境中逼真地模拟脑部的变形过程，满足了人脑变形仿真中的精度与速度要求。1、本专利技术技术解决方案：一种基于多GPU的人脑变形仿真方法，其特征在于包括如下步骤：(1)初始化，读入人脑四面体网格模型数据，分析四面体网格单元的拓扑结构关系得到四面体网格邻接体单元的集合，依据所获得的集合信息和所用GPU核数目，对读入的四面体网格模型原始数据进行分组、编号，然后基于等参元法计算各四面体单元的形函数矩阵，用于后面在发生受力或表面位移后内部单元的形变计算；(2)初始化完成后，根据手动输入或是力反馈设备捕获的边界受力变化或是扫描仪或相机捕获的可见表面节点的位移变化，利用Neo-Hookean模型，该模型定义了力与顶点位移之间的关系，总体拉格朗日方法，计算每一个时间步的组成四面体单元的各顶点的受力，在算法的实现过程中在多GPU服务器上并行实现，通过这种处理方式来减少多个GPU之间的数据交互；(3)把步骤二中每个GPU核计算出的受力结果进行合并，并根据最终计算结果采用中心差分法公式其中Δt表示时间步，M表示四面体单元的质量，un+1表示下一个时间步的位移，Rn表示当前时刻的反作用力，Fn表示当前时刻所施加的外力即步骤二捕获的表面受力，un-1表示前一时刻的位移，un表示当前时刻的位移即步骤二中捕获的表面节点的位移变化，更新四面体单元各顶点的坐标变化，并在多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应力公式S0t=μ×J-2/3(I-13I1×C0t-1)+δ(J-1)×J×C0t-1]]>其中S表示第二Piola–Kirchoff应力，μ表示剪切模量，J表示变形梯度，I是三阶单位矩阵，I1是右柯西格林张量第一不变量，C是右柯西格林张量，δ是材料的体积弹性模量，通过该公式计算各单元的受力情况；(3.4)通过初始化计算时生成的各四面体单元顶点受力求和求整个四面体网格模型各顶点的总受力，并将结果传送回CPU端。所述步骤(2)中总体拉格朗日方法实现过程为：方法在预测模型变形时通过把整个连续的变形过程离散为若干个时刻的变形，设计算开始的时间为t，计算变形的下一个时刻和时间t的时间间隔叫做时间步，记做Δt，首先使用顶点的初始坐标、步骤(本文档来自技高网...

【技术保护点】
一种基于多GPU的人脑变形仿真方法，其特征在于包括如下步骤：(1)初始化，读入人脑四面体网格模型数据，分析四面体网格单元的拓扑结构关系得到四面体网格邻接体单元的集合，依据所获得的集合信息和所用GPU核数目，对读入的四面体网格模型原始数据进行分组、编号，然后基于等参元法计算各四面体单元的形函数矩阵，用于后面在发生受力或表面位移后内部单元的形变计算；(2)初始化完成后，根据手动输入或是力反馈设备捕获的边界受力变化或是扫描仪或相机捕获的可见表面节点的位移变化，利用Neo‑Hookean模型，该模型定义了力与顶点位移之间的关系，总体拉格朗日方法，计算每一个时间步的组成四面体单元的各顶点的受力，在实现过程中采用多GPU加速计算过程，进行数据分组减少多GPU数据交互；(3)把步骤二中每个GPU核计算出的受力结果进行合并，并根据最终计算结果采用中心差分法公式(1Δt2M)un+1=Rn-ΣiFni-1Δt2M(un-1-2un),]]>其中Δt表示时间步，M表示四面体单元的质量，un+1表示下一个时间步的位移，Rn表示当前时刻的反作用力，Fn表示当前时刻所施加的外力即步骤二捕获的表面受力，un‑1表示前一时刻的位移，un表示当前时刻的位移即步骤二中捕获的表面节点的位移变化，更新四面体单元各顶点的坐标变化，并在多CPU上实现。...

【技术特征摘要】
1.一种基于多GPU的人脑变形仿真方法，其特征在于包括如下步骤：
(1)初始化，读入人脑四面体网格模型数据，分析四面体网格单元的拓扑结构关系得
到四面体网格邻接体单元的集合，依据所获得的集合信息和所用GPU核数目，对读入的四
面体网格模型原始数据进行分组、编号，然后基于等参元法计算各四面体单元的形函数矩
阵，用于后面在发生受力或表面位移后内部单元的形变计算；
(2)初始化完成后，根据手动输入或是力反馈设备捕获的边界受力变化或是扫描仪或
相机捕获的可见表面节点的位移变化，利用Neo-Hookean模型，该模型定义了力与顶点位
移之间的关系，总体拉格朗日方法，计算每一个时间步的组成四面体单元的各顶点的受力，
在实现过程中采用多GPU加速计算过程，进行数据分组减少多GPU数据交互；
(3)把步骤二中每个GPU核计算出的受力结果进行合并，并根据最终计算结果采用中
心差分法公式(1Δt2M)un+1=Rn-ΣiFni-1Δt2M(un-1-2un),]]>其中Δt表示时间步，M表示四面
体单元的质量，un+1表示下一个时间步的位移，Rn表示当前时刻的反作用力，Fn表示当前时
刻所施加的外力即步骤二捕获的表面受力，un-1表示前一时刻的位移，un表示当前时刻的位
移即步骤二中捕获的表面节点的位移变化，更新四面体单元各顶点的坐标变化，并在多
CPU上实现。
2.根据权利要求1所述的基于多GPU的人脑变形仿真方法，其特征在于：步骤(1)中的
对读入的四面体网格模型原始数据进行分组、编号，然后基于等参元法计算各四面体单元
的形函数矩阵，包括步骤如下：
(2.1)遍历所有四面体网格单元的顶点，为每一个顶点建立所属单元列表，为每个顶点
建立访问标志位并设置为false；
(2.2)输入分组起始参考顶点，对四面体网格模型单元进行分组，分组单元数＝四面体总
单元数/所用GPU核数，首先将起始点的所属单元加入分组，并将顶点访问标志设为true，
判断当前分组内包含四面体是否达到分组单元数限制，达到则从下一个输入的四面体单元
顶点开始，将四面体放入下一分组，否则读取参考点所属单元列表内访问标志顶点为false
的节点的所属单元并加入当前分组中；
(2.3)步骤(2.2)对四面体单元进行分组完成后，将各组内顶点根据所属单元列表是否完

\t全在当前分组内进行区分，满足条件的为组内点，否则为边界点；对所有分组内四面体单
元的顶点依照分区次序从低到高，...

【专利技术属性】
技术研发人员：胡勇，田野，沈旭昆，
申请(专利权)人：北京航空航天大学，
类型：发明
国别省市：北京;11