一种软组织形变仿真方法技术

本发明专利技术涉及一种基于光滑粒子流体动力学的软组织形变仿真方法，属于图形处理技术领域，该方法选取光滑粒子流体动力学法，以黏弹性力学模型来反映软组织的生物力学特性，包含以下步骤：依据黏弹性模型，构建软组织形变仿真计算相关的一系列方程；选择合适的支持域搜索策略和光滑核函数，采用粒子近似法对方程组的各相关项进行近似计算，通过显示积分法计算各粒子的密度、位置、速度等随时间的变化值；动态将粒子模型每个时间步长的状态输出到屏幕上，并进行纹理光照的渲染，显示软组织器官受力情况下的实时形变过程。该方法无需繁琐的网格计算，可提高软组织形变仿真的准确性和实时性。

【技术实现步骤摘要】

本专利技术涉及虚拟手术
，具体地说是一种基于光滑粒子流体动力学的软组 织形变仿真方法，属于图形处理
的物体形变实时模拟
，其IPC专利分类 号为 G06T17/00。
技术介绍
虚拟手术是虚拟现实技术在现代医学中的应用，是一个多学科交叉的研究领域。 它主要由医学数据可视化与建模、人体软组织器官受力形变仿真两部分构成，在视觉上与 触觉感官上为使用者提供了手术场景的真实再现。可用于外科医生的培训、手术结果的预 测、手术计划的辅助制定、手术导航等，具有重大的理论意义和广阔的应用前景。软组织形变仿真是虚拟手术中的最重要技术之一，能否逼真、实时的模拟软组织 器官在外力作用下的形变是整个系统的关键。常见的软组织形变计算模型分为两大类基于几何的形变模型和基于物理的形 变模型。其中，第一类模型仅仅考虑了几何形态的变化，而忽略了软组织的实际力学本构 方程以及形变过程中物体质量、力或其他物理现象的作用，因此不能真实的反映软组织的 形变过程，该模型目前已较少使用。而物理模型则基于软组织的力学本构方程，通过相应 的计算模型得出组织受力时的形变，能够更加真实的反映组织的形变，因此目前使用较多。 质点-弹簧模型(Mass-Spring Model)是较早出现的物理模型，它由无质量的弹簧和阻尼 器网格连接的质点集合组成，该模型建模简单、易于实现，且计算速度相对较快；但很难通 过实验来获取模型中成千上万的弹簧、质点和阻尼器参数以反映软组织内部真实的力学特 性，因此，使用该模型实现形变仿真的逼真度不高；另外，前期需要耗费很多精力来构造弹 簧阻尼器网格，仿真计算对网格结构依赖程度大，并且不能准确描述大变形问题。有限元模 型(Finite Element Model，简写为FEM)是另一种常用的有网格物理模型，它将软组织离 散为若干子单元，通过这些子单元边界上的节点相互结成为组合体，并利用变分原理建立 求解应变量的代数方程组，从而计算形变。用有限元模型实现软组织形变仿真的方法已在 近几年得到广泛的研究和应用，该方法具有坚实的弹性力学基础，较之其他有网格模型计 算精度高，但它的计算复杂度大，很难实现软组织形变的实时模拟。边界元模型(Boundary Element Model，简写为BEM)可以看作对有限元模型的改进，它只对模型的边界进行离散， 以降低问题的维数，简化计算。该方法无需考虑内部节点位移，较有限元法计算简单，但这 种方法只能解决具有均质各向同性的线性问题。形变仿真中软组织所采用的力学模型是形变仿真的基础与关键，决定形变仿真的 准确性。当前常见的仿真系统为了能够实现组织的实时形变仿真，大多采用较为简化的力 学模型，如弹性力学模型、线弹性力学模型等，这些类型的模型大多直接采用经典材料力 学中的力学模型，与软组织的实际力学特性差异很大，因此尽管可以实现能够满足一定速 度需求的仿真，但精度很低，形变效果不如人意。如前所述，现在常用的有网格物理模型需要耗费很大的精力来构造网格模型，后续的计算过程大都紧密的依赖于这种网络结构，并且计算复杂度往往较大，因此无法满足 实时性、健壮性、精确性的要求，更难以准确描述软组织大变形问题；另外，经典的、简化的 材料力学模型不能真实反映软组织的生物力学特性，依据这种力学模型对软组织的形变进 行仿真计算，逼真度也很低。软组织形变仿真是虚拟手术系统中的关键环节，为此，如何设 计一个性能更加良好的方法来实现软组织的形变仿真，以尽可能的满足系统关于实时性、 健壮性、精确性的要求，成为虚拟手术系统所面临的首要问题。
技术实现思路
本专利技术的目的是解决虚拟手术系统中软组织形变仿真的实时性、逼真性和大变形 仿真问题，克服了现有技术的不足，以较能真实反映软组织力学特性的黏弹性模型为力学 模型，提供一种基于无网格法的软组织形变仿真方法，并实现其中最关键的基础部分，使其 满足虚拟手术系统的需要。为完成本专利技术的目的，本专利技术采取的技术方案是以黏弹性力学模型反映软组织 的生物力学特性，采用无网格连接的粒子模型来表示软组织形变计算模型，并选取光滑粒 子流体动力学法(Smoothed Particle Hydrodynamics，简写为SPH)作为无网格模型的计算 方法，完成软组织形变的仿真计算。其包含以下几个步骤步骤1)采集软组织的数据信息；步骤幻选择黏弹性模型，构建用于软组织形变仿真的方程组应力-应变本构方 程(1)、O)、(3)、应变-位移几何方程(4)和用于加速度计算的动量方程(5)；步骤3)依据步骤1)采集的软组织的数据信息，利用各数据点的位置向量，构建 一个没有网格连接的粒子模型，并初始化所述粒子模型中各粒子的位置、质量、速度、加速 度、受作用力等信息，构建粒子模型的初始化状态；步骤4)定义空间栅格，采用链表搜索法搜索粒子的支持域，构建支持域内的光 滑核函数W;核函数W选取三次B-Spline光滑函数(6)，对核函数在各方向求导，得到核函 数的一阶导数⑵；步骤幻应用光滑核函数W及其导数对参考粒子的支持域内所有粒子函数的加权 平均近似的方法来构建步骤2、中各方程的SPH格式构建密度方程SPH格式(8)；将动量方程(5)转换为SPH格式(9)；构建几何方程的SPH格式(10)；步骤6)用显示积分法求解步骤幻的常微分方程，计算出粒子的密度、位置、速度 等随时间的变化值；步骤7)循环执行上述步骤幻 步骤6)，计算出各粒子的状态；步骤8)将粒子模型的当前状态输出之屏幕，经过纹理和光照渲染后，得到软组 织的动态形变过程。进一步地，步骤6)的求解过程如下1)对模型的某一质点ρ施加外力f。utp ；2)对模型中的各粒子Pi (i = 1,2,..，M)循环完成下述3) 7)步骤的运算处理；3)对当前粒子Pi，以h为光滑半径，搜索支持域内的相邻粒子Pj(j = 1,2,.. ,N)；4)使用公式(6) (7)，计算当前粒子Pi与支持域内各邻近粒子Pj (j = 1,2, ..,N) 间的光滑核函数Wu及其导数；5)使用方程(8)，采用密度求和法计算粒子的密度P i ；6)计算粒子的加速度^，加速度的计算采用如下方法权利要求1. 和技术，其特征在于包含以下步骤 步骤1)采集软组织的数据信息；步骤幻选择黏弹性模型，构建用于软组织形变仿真的方程组 采用一个弹簧和一个黏壶并联的Kelvin黏弹性模型； 首先，构建三维格式的Kelvin黏弹性应力-应变本构方程依据Kelvin模型，各向同性材料的应力张量ο可分解成它的球形张量和偏斜张量部 分，应变张量ε可分离为体积形变和等体积的形状畸变两部分2.根据权利要求1的软组织形变仿真方法，其特征在于步骤6)的求解过程如下1)对模型的某一质点P施加外力f。utp；2)对模型中的各粒子Pi(i = 1，2，. .，M)循环完成下述3) 7)步骤的运算处理；3)对当前粒子Pi，以h为光滑半径，搜索支持域内的相邻粒子=1,2,.. ,N)；4)使用公式(6)(7)，计算当前粒子Pi与支持域内各邻近粒子= 1，2，..，N)间的 光滑核函数Wu及其导数；5)使用方程(8)，采用密度求和法计算粒子的密度Pi ；6)计算粒子的加速度％，加速度的计算采用如下方法全文摘要本专利技术涉及一种基于光滑粒子流体动力本文档来自技高网...

【技术保护点】
一种软组织形变仿真方法和技术，其特征在于包含以下步骤：步骤１）：采集软组织的数据信息；步骤２）：选择黏弹性模型，构建用于软组织形变仿真的方程组：采用一个弹簧和一个黏壶并联的Ｋｅｌｖｉｎ黏弹性模型；首先，构建三维格式的Ｋｅｌｖｉｎ黏弹性应力－应变本构方程：依据Ｋｅｌｖｉｎ模型，各向同性材料的应力张量σ可分解成它的球形张量和偏斜张量部分，应变张量ε可分离为体积形变和等体积的形状畸变两部分：σ↓［αβ］＝Ｓ↓［αβ］＋δ↓［αβ］σ↓［ｋｋ］／３（１）ε↓［αβ］＝ｅ↓［αβ］＋δ↓［αβ］ε↓［ｋｋ］／３（２）其中，α，β＝ｘ，ｙ，ｚ．σ↓［αβ］为应力分量，ε↓［αβ］为应变分量；δ↓［αβ］为Ｋｒｏｎｅｃｋｅｒ符号，σ↓［ｋｋ］＝σ↓［ｘｘ］＋σ↓［ｙｙ］＋σ↓［ｚｚ］和ε↓［ｋｋ］＝ε↓［ｘｘ］＋ε↓［ｙｙ］＋ε↓［ｚｚ］分别为体积应力和体积应变；Ｓ↓［αβ］和ｅ↓［αβ］分别为偏应力张量和偏应变张量的分量；根据Ｋｅｌｖｉｎ模型，偏应力张量和偏应变张量之间、体积应力和体积应变之间的三维黏弹性本构关系可表示为：Ｓ↓［αβ］＝Ｅ.ｅ↓［αβ］＋η.ｄｅ↓［αβ］／ｄｔ，α，β＝ｘ，ｙ，ｚσ↓［ｋｋ］＝Ｅ.ε↓［ｋｋ］＋η.ｄε↓［ｋｋ］／ｄｔ（３）上式中，Ｓ↓［αβ］和ｅ↓［αβ］分别为偏应力张量和偏应变张量的分量；σ↓［ｋｋ］和ε↓［ｋｋ］分别为体积应力和体积应变；Ｅ为材料的弹性模量，η为黏性系数；ｄｅ↓［αβ］／ｄｔ，ｄε↓［ｋｋ］／ｄｔ分别为偏应变分量和体积应变对时间的导数，即应变率，ｔ为时间；其次，构建应变－位移几何方程：ε↓［ｘｘ］＝＊ｕ／＊ｘε↓［ｘｙ］＝＊ｕ／＊ｙ＋＊ｖ／＊ｘε↓［ｙｙ］＝＊ｖ／＊ｙε↓［ｙｚ］＝＊ｖ／＊ｚ＋＊ｗ／＊ｙε↓［ｚｚ］＝＊ｗ／＊ｚε↓［ｚｘ］＝＊ｗ／＊ｘ＋＊ｕ／＊ｚ（４）其中ｕ，ｖ，ｗ为位移在三个坐标方向的分量；ε↓［ｘｘ］，ε↓［ｙｙ］，ε↓［ｚｚ］为正应变分量；ε↓［ｘｙ］，ε↓［ｙｚ］，ε↓［ｚｘ］为剪应变分量；然后，构建用于加速度计算的动量方程，其形式如下：ｄｖ／ｄｔ＝１／ρ＊σ／＊ｘ（５）其中，ｖ为速度向量，ｔ为时间，ρ为粒子的密度，ｘ为坐标向量；步骤３）：依据步骤１）采集的软组织的数据信息，利用各数据点的位置向量，构建一个没有网格连接的粒子模型，并初始化所述粒子模型中各粒子的位置、质量、速度、加速度、受作用力等信息，构建粒子模型的初始化状态；步骤４）：定义空间栅格...

【技术特征摘要】

【专利技术属性】
技术研发人员：刘雪梅，皇甫中民，向明森，闫雒恒，赵振国，赵晶，闫新庆，吴慧欣，郝爱民，刘明堂，孙新娟，杨礼波，石秋华，刘欢，
申请(专利权)人：华北水利水电学院，
类型：发明
国别省市：41[中国|河南]