[发明专利]一种软组织形变仿真方法

摘要

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

主权项

1.一种软组织形变仿真方法和技术，其特征在于包含以下步骤：步骤1)：采集软组织的数据信息；步骤2)：选择黏弹性模型，构建用于软组织形变仿真的方程组：采用一个弹簧和一个黏壶并联的Kelvin黏弹性模型；首先，构建三维格式的Kelvin黏弹性应力-应变本构方程：依据Kelvin模型，各向同性材料的应力张量σ可分解成它的球形张量和偏斜张量部分，应变张量ε可分离为体积形变和等体积的形状畸变两部分：σαβ=Sαβ+δαβσkk3---(1)]]>ϵαβ=eαβ+δαβϵkk3---(2)]]>其中，α，β＝x，y，z.σ_αβ为应力分量，ε_αβ为应变分量；δ_αβ为Kronecker符号，σ_kk＝σ_xx+σ_yy+σ_zz和ε_kk＝ε_xx+ε_yy+ε_zz分别为体积应力和体积应变；S_αβ和e_αβ分别为偏应力张量和偏应变张量的分量；根据Kelvin模型，偏应力张量和偏应变张量之间、体积应力和体积应变之间的三维黏弹性本构关系可表示为：Sαβ=E·eαβ+η·deαβdt,α,β=x,y,z]]>σkk=E·ϵkk+η·dϵkkdt---(3)]]>上式中，S_αβ和e_αβ分别为偏应力张量和偏应变张量的分量；σ_kk和ε_kk分别为体积应力和体积应变；E为材料的弹性模量，η为黏性系数；分别为偏应变分量和体积应变对时间的导数，即应变率，t为时间；其次，构建应变-位移几何方程：ϵxx=∂u∂x]]>ϵxy=∂u∂y+∂v∂x]]>ϵyy=∂v∂y]]>ϵyz=∂v∂z+∂w∂y---(4)]]>ϵzz=∂w∂z]]>ϵzx=∂w∂x+∂u∂z]]>其中u，v，w为位移在三个坐标方向的分量；ε_xx，ε_yy，ε_zz为正应变分量；ε_xy，ε_yz，ε_zx为剪应变分量；然后，构建用于加速度计算的动量方程，其形式如下：dvdt=1ρ∂σ∂x---(5)]]>其中，v为速度向量，t为时间，ρ为粒子的密度，x为坐标向量；步骤3)：依据步骤1)采集的软组织的数据信息，利用各数据点的位置向量，构建一个没有网格连接的粒子模型，并初始化所述粒子模型中各粒子的位置、质量、速度、加速度、受作用力等信息，构建粒子模型的初始化状态；步骤4)：定义空间栅格，采用链表搜索法搜索粒子的支持域，构建支持域内的光滑核函数；核函数选取三次B-Spline光滑函数：Wij=1πh3×1-32s2+34s30≤s<114(2-s)31≤s<20s≥2---(6)]]>式中，W_ij为由邻近粒子p_j估计粒子p_i运动信息的光滑核函数，为粒子p_i与p_j之间的相对距离，r为粒子p_i与p_j之间的距离，h为光滑长度；对核函数在各方向求导，得到核函数的一阶导数为：∂Wij∂xβ=xiβ-xjβr·1πh4×-3s+94s20≤s<1-34(2-s)21≤s<20s≥2---(7)]]>式中，上标β＝x，y，z表示坐标方向，s，r，h的含义同式(6)，和分别表示粒子p_i与p_j的位置坐标向量在各方向的分量；步骤5)：应用光滑核函数W及其导数对参考粒子的支持域内所有粒子函数的加权平均近似的方法来构建步骤2)中各方程的SPH格式：构建密度方程SPH格式ρi=Σj=1NmjWij---(8)]]>其中，ρ_i为粒子p_i的密度，m_j为p_i的支持域内邻近粒子p_j的质量，N为粒子p_i的支持域内粒子总数，W_ij的含义同式(6)；将动量方程(5)转换为SPH格式DvixDt=Σj=1N[mj(σixxρi2+σjxxρj2)∂Wij∂x+mj(σixyρi2+σjxyρj2)∂Wij∂y+mj(σixzρi2+σjxzρj2)∂Wij∂z]]]>DviyDt=Σj=1N[mj(σiyyρi2+σjyyρj2)∂Wij∂y+mj(σiyxρi2+σjyxρj2)∂Wij∂x+mj(σiyzρi2+σjyzρj2)∂Wij∂z]---(9)]]>DvizDt=Σj=1N[mj(σizzρi2+σjzzρj2)∂Wij∂z+mj(σizxρi2+σjzxρj2)∂Wij∂x+mj(σizyρi2+σjzyρj2)∂Wij∂y]]]>其中，分别为粒子p_i的运动速度在各坐标方向的分量；为粒子p_i应力向量的正应力分量，为粒子p_i应力向量的剪应力分量；为粒子p_j应力向量的正应力分量，为粒子p_j应力向量的剪应力分量；ρ_i，ρ_j分别为粒子p_i和p_j的密度；m_j，W_ij，N的含义同公式(8)；构建几何方程的SPH格式：ϵixx=∂ui∂x=Σj=1Nmjρj·uj·∂Wij∂x]]>ϵiyy=∂vi∂y=Σj=1Nmjρj·vj·∂Wij∂y]]>ϵizz=∂wi∂z=Σj=1Nmjρj·wj·∂Wij∂z---(10)]]>ϵixy=∂ui∂y+∂vi∂x=Σj=1Nmjρj·uj·∂Wij∂y+Σj=1Nmjρj·vj·∂Wij∂x]]>ϵiyz=∂vi∂z+∂wi∂y=Σj=1Nmjρj·vj·∂Wij∂z+Σj=1Nmjρj·wj·∂Wij∂y]]>ϵizx=∂wi∂x+∂ui∂z=Σj=1Nmjρj·wj·∂Wij∂x+Σj=1Nmjρj·uj·∂Wij∂z]]>其中，为粒子p_i应变向量的正应变分量，为粒子p_i应变向量的剪应变分量；u_i，v_i，w_i为粒子p_i的运动位移在各坐标方向的分量，u_j，v_j，w_j为粒子p_j的运动位移在各坐标方向的分量；ρ_j为粒子p_j的密度，m_j，W_ij，N的含义同式(8)；步骤6)：用显示积分法求解步骤5)的常微分方程，计算出粒子的密度、位置、速度等随时间的变化值；步骤7)：循环执行上述步骤3)～步骤6)，计算出各粒子的状态；步骤8)：将粒子模型的当前状态输出之屏幕，经过纹理和光照渲染后，得到软组织的动态形变过程。