一种软组织形变的局部质点弹簧模型的实时动态生成方法

摘要

本发明属于计算机图形学领域，具体地说是一种虚拟手术系统中模软组织形变的局部质点弹簧模型的实时动态生成方法。本发明定义一个用于选择局部区域的圆形拾取工具；根据需要，改变圆形拾取工具的位置和尺寸；将当前场景中的软组织模型投影至圆形拾取工具所在平面，根据投影模型和圆的位置关系，提取出位于圆形拾取工具区域内的顶点和三角形集合；再顶点集合的边界点分离出来；最后实时生成模型，并将边界点固定。该发明很大程度地提高了系统的实时性，为力觉反馈的实时更新奠定了基础。此外本发明的方法简单，速度快，可有效地用于虚拟手术系统的软组织形变仿真中。

说明书

技术领域

本发明属于计算机图形学领域，具体地说是一种虚拟手术系统中模软组织形变的局部质点弹簧模型的实时动态生成方法。 

背景技术

虚拟现实技术涉及计算机图形学、人机交互技术、传感技术、人工智能等领域，它用计算机生成逼真的三维视、听、嗅觉等感觉，使人作为参与者通过适当装置，自然地对虚拟世界进行体验和交互作用。 

虚拟现实技术在医学方面的应用具有十分重要的现实意义。传统手术训练使用人类的尸体、动物或橡胶人体模型作为训练对象，存在着诸多局限性，虚拟手术技术的出现成为解决此类问题的重要手段。这种手术训练系统能够对手术全过程进行训练，这将极大地提高外科手术的学习效率、降低手术训练成本。另外，在远距离遥控外科手术，复杂手术的计划安排，手术过程的信息指导，手术后果预测及改善残疾人生活状况，乃至新药研制等方面，虚拟现实技术都能发挥十分重要的作用。 

在模拟虚拟手术系统中软组织形变过程的时候，通常使用的物理建模策略是质点——弹簧模型（简记为MSS），质点——弹簧模型，顾名思义是基于虚拟质点及弹簧建立起来的一种物理模型。设质点数为n，其相应的弹簧连接数量与质点之间的连接数量一致。其中（1）是系统的形变微分方程： 

$\left(\begin{array}{c}{X}^{\prime }\\ V^{\prime }\end{array}\right)=\left(\begin{array}{c}V\\ F/M\end{array}\right)---\left(1\right)$

其中，X为质点位移；V为质点运动速度；F为虚拟场景中当前受作用的质点x方向、y方向、z方向的合力；作用力F为内力F_int与外力F_ext之和，在模型发生变形时所产生的弹簧力F_i和阻尼力F_d构成内力F_int。因此，作用力F的表达式可以写成如下形式： 

F＝F_int+F_ext＝F_s+F_d+F_ext＝-KX-DX′+F_ext    （2） 

由式(1)和式(2)可知 

MX″+DX′+KX＝F_ext    （3） 

其中，D、M、K分别为阻尼矩阵、质量矩阵及刚度矩阵，均为n阶对角方阵；且K为系数带状矩阵；为1阶列矩阵，元素个数为n。表示质点所受外力的合力。在虚拟场景中的任意一个质点表示为i，应当满足下列方程: 

m_ia＝-f_s-f_d+f_ext    （4） 

其中: 

$f_s=\underset{j\in P\left(i\right)}{\Sigma }{k}_{\mathrm{ij}}\frac{\left(\right|\left|l_{\mathrm{ij}}\right||-l_{\mathrm{ij}}^0)l_{\mathrm{ij}}}{\left|\right|l_{\mathrm{ij}}\left|\right|}---\left(5\right)$

$f_d=\underset{j\in P\left(i\right)}{\Sigma }{c}_{\mathrm{ij}}\left|\right|(v_i-v_j)·\frac{l_{\mathrm{ij}}}{\left|\right|l_{\mathrm{ij}}\left|\right|}\left|\right|\frac{l_{\mathrm{ij}}}{\left|\right|l_{\mathrm{ij}}\left|\right|}=\underset{j\in P\left(i\right)}{\Sigma }{c}_{\mathrm{ij}}\left|\right|(x_i^{\prime }-x_j^{\prime })·\frac{l_{\mathrm{ij}}}{\left|\right|l_{\mathrm{ij}}\left|\right|}\left|\right|\frac{l_{\mathrm{ij}}}{\left|\right|l_{\mathrm{ij}}\left|\right|}---\left(6\right)$

将式（5）和（6）带入式(4)得： 

$m_{i}a+\underset{j\in P\left(i\right)}{\Sigma }{k}_{\mathrm{ij}}\frac{\left(\right|\left|l_{\mathrm{ij}}\right||-l_{\mathrm{ij}}^0)l_{\mathrm{ij}}}{\left|\right|l_{\mathrm{ij}}\left|\right|}+\underset{j\in P\left(i\right)}{\Sigma }{c}_{\mathrm{ij}}\left|\right|(x_i^{\prime }-x_j^{\prime })·\frac{l_{\mathrm{ij}}}{\left|\right|l_{\mathrm{ij}}\left|\right|}\left|\right|\frac{l_{\mathrm{ij}}}{\left|\right|l_{\mathrm{ij}}\left|\right|}=f_{\mathrm{ext}}---\left(7\right)$

整理得: 

$f_{\mathrm{ext}}-m_i{x}_i^{\prime \prime }=\underset{j\in P\left(i\right)}{\Sigma }(k_{\mathrm{ij}}\left(\right|\left|l_{\mathrm{ij}}\right||-l_{\mathrm{ij}}^0)+c_{\mathrm{ij}}||(x_i^{\prime }-x_j^{\prime })·\frac{l_{\mathrm{ij}}}{\left|\right|l_{\mathrm{ij}}\left|\right|}|\left|\right)\frac{l_{\mathrm{ij}}}{\left|\right|l_{\mathrm{ij}}\left|\right|}---\left(8\right)$

这里m_i为单个质点i的质量；k_ij为两相邻质点间弹簧系数，ij表示i质点与j质点（如下ij等同于此意）；c_ij为两相邻质点的阻尼系数；P(i)为与质点i所有相邻质点的集合；l_ij为两相邻质点距离；为两相邻质点初始距离；质点i运动时的加速度用a表示。 

结合上文，质点i的运动位移和速度可用如下微分方程确定： 

$\left(\begin{array}{c}{x}^{\prime }=v\\ v^{\prime }=\frac{f_{\mathrm{ext}}-f_s-f_d}{m}\end{array}\right)---\left(9\right)$

简单的质点——弹簧连接模型如图7示。 

投影变换是计算机图形学中的一个重要应用，投影变换的目的是定义一个可视空间。它决定了一个物体如何映射到屏幕上（即通过透视投影或正投影），并且定义了哪些物体或物体的哪些部分被裁减于最终的图像之外。透视投影最显著的特征就是透视缩短：物体距离照相机越近，它在最终图像上看上去就越小。这是因为透视投影的可视空间是一个金字塔的平截头体（顶部被一个平行于底面的平面所截除）。位于可视空间之内的物体将被投影到金字塔的顶点，也就是照相机或观察点的位置，如图8所示。 

计算机图形学中的几何处理管线擅长于使用视图和投影矩阵以及用于裁剪的视口把顶点世界（物体）坐标变换为窗口（屏幕）坐标。这里，我们还需要逆转这个过程，即用鼠标选择三维空间中的一个位置，鼠标返回一个二维值，也就是鼠标光标所在的屏幕位置。因此，应用程序必须逆转变换过程，确定这个屏幕位置源于三维空间的什么地方。 

发明内容

本发明的目的在于提供一种用于虚拟手术系统中模拟软组织形变的局部质点弹簧模型的实时动态生成方法。 

本发明具体包括如下步骤： 

（1）获取当鼠标左键按下时，鼠标在屏幕坐标下的位置（x，y），选取近侧裁剪平面z＝z_p为投影面，将屏幕坐标（x，y）变换为三维空间坐标（x，y，z_p），然后再以点（x，y，z_p）为中心、长度R_init为半径，进行局部区域的圆形拾取，进行图形渲染，投影面即是透视投影中的近裁剪面z_p； 

（2）定义拾取区域，动态更新拾取区域的位置和尺寸，设此时圆的中心为（x_new，y_new，z_p），半径为R_new； 

（3）将软组织模型的顶点投影至近裁剪平面z_p，以判断投影后的顶点和步骤（2）中以（x_new，y_new，z_p)为中心、R_new为半径的圆的位置关系：1）将软组织模型的全部顶点投影至屏幕坐标下，坐标为（x_i，y_i，z_i）；2）屏幕坐标下的坐标投影至近裁剪面z_p，坐标分别为（xp_i，yp_i，z_p），遍历所有的点（xp_i，yp_i，z_p），根据软组织投影点和圆的位置关系，提取出位于圆内的点集和三角形集合，分别设为P和T； 

（4）根据集合T中三角形的网格拓扑关系，从顶点集合P中分离出边界点集，设为集合B： 

1）建立一个临时点数组Temp，并置为空； 

2）对当前点的所有邻接三角网格上异于当前点的点进行判断：如果临时点数组Temp包含此点，则从临时点数组Temp中删除此点，否则，把此点加入到临时点数组Temp中； 

3）检查临时点数组Temp，如果临时点数组Temp为空，则当前点为内点，否则为边界点，并且临时点数组Temp中的点是和当前点邻接的边界点； 

（5）对位于拾取区域内的三角形集合T，根据顶点之间的连接关系，实时生成局部质点弹簧模型，并将边界点集合B的顶点固定： 

1）遍历顶点集合P，把P中的所有顶点添加至质点弹簧系统中； 

2）遍历三角形集合T中的所有边，若当前边i尚未生成弹簧模型，则对该边添加弹簧，并设置当前边是否已加弹簧的标志位为true； 

3）遍历边界点集合B，并将B中的所有顶点固定。 

本发明与现有技术相比的优势在于：现有的虚拟手术中，软组织形变仿真都是对整个模型进行物理建模，系统实时性差；或者在程序初始化阶段，针对某一固定区域，根据三角形 的拓扑关系，指定该区域的大小，来生成局部质点——弹簧模型，该方法一旦生成局部质点弹簧模型，在程序运行过程中就无法对其进行更改，而且初始化时的更新操作又特别繁琐，对非技术人员来说，这无疑是一个弊病。该发明在程序初始化阶段无须进行任何操作，在程序运行过程中可以实现局部区域可视化选择、动态生成局部质点——弹簧模型，具有较强的灵活性和实用价值，更加符合使用者的操作习惯。该发明很大程度地提高了系统的实时性，为力觉反馈的实时更新奠定了基础。此外本发明的算法简单，速度快，可有效地用于虚拟手术系统的软组织形变仿真中。 

附图说明

图1是本发明方法的基本流程图； 

图2是用于选择局部区域的圆形拾取工具示意图； 

图3是同一尺寸的圆形拾取工具位于不同位置； 

图4是同一位置的圆形拾取工具拥有不同的尺寸； 

图5是对图3中的模型动态生成的局部质点——弹簧模型； 

图6是对图4中的模型动态生成的局部质点——弹簧模型； 

图7是本发明的系统仿真图； 

图8为本发明的系统测试图； 

图9为本发明的系统点阵模拟图。 

通过附图可以看出，本发明具有实时、动态生成局部质点——弹簧模型的特点，其灵活性是传统方法所无法比拟的。 

具体实施方式

传统的虚拟手术在模拟软组织形变过程的时候，都是对整个模型进行物理建模，由于质点弹簧模型在每个时间步长内进行数值计算的复杂程度与质点——弹簧模型的规模呈正比例关系。这就迫切需要提出一种新的局部质点——弹簧模型的动态生成方法。 

（1）首先定义一个用于选择局部区域的圆形拾取工具； 

（2）根据需要，改变（1）中圆形拾取工具的位置和尺寸； 

（3）将当前场景中的软组织模型投影至（2）中圆形拾取工具所在平面，根据投影模型和圆的位置关系，提取出位于圆形拾取工具区域内的顶点和三角形集合； 

（4）再将（3）中顶点集合的边界点分离出来； 

（5）最后根据位于圆内的三角形之间的拓扑关系，实时生成局部质点——弹簧模型，并将（4）中边界点固定。 

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

（1）首先利用OpenGL中的函数gluMouseFunc(HandleMouseButton)和自定义函数void HandleMouseButton（int button，int state，int x,int y)获取当鼠标左键按下时，鼠标在屏幕 坐标下的位置（x，y），选取近侧裁剪平面z＝z_p为投影面，利用三维计算机图形投影变换函数gluUnProtect()将屏幕坐标（x，y）变换为三维空间坐标（x，y，z_p），然后再以点（x，y，z_p）为中心、长度R_init为半径，定义一个用于选取局部区域的圆形拾取工具，如图2所示； 

gluUnProtect()函数的说明： 

int gluUnProtect(GLdouble winx,GLdouble winy,GLdouble winz，const 

GLdouble modelMatrix[16],const GLdouble projMatrix[16]，const GLint 

viewport[4]，GLdouble*objx，GLdouble*objy，GLdouble*objz)； 

这个函数使用由模型视图矩阵（modelMatrix）、投影矩阵（projMatrix）和视口（viewport）所定义的变换把指定的窗口坐标（winx，winy，winz）映射到物体坐标。它所产生的物体坐标在objx，objy，objz中返回。如果程序中使用的是glDepthRange()的默认值，winz为0.0表示经过变换这个的点的全部坐标位于近裁剪面上，如果winz是1.0表示它的坐标位于远裁剪面上。本发明中使用的是glDepthRange()的默认值，并且winz为0.0，故投影面即是透视投影中的近裁剪面，即objz=z_p。 

（2）自定义拾取区域环节，在鼠标按下的状态，通过移动鼠标来控制（1）中定义的用于选择局部区域的圆形拾取工具的位置，并通过键盘按键来增大或减小其半径，以动态更新拾取区域的位置和尺寸，设此时圆的中心为（x_new，y_new，z_p），半径为R_new，如图3和图4所示； 

（3）当鼠标抬起时，将软组织模型的顶点投影至近裁剪平面（即投影面）z_p，以判断投影后的顶点和（2）中以（x_new，y_new，z_p）为中心、R_new为半径的圆的位置关系。这一过程主要分为以下两步：①将软组织模型的全部顶点使用函数gluProtect()函投影至屏幕坐标下，设坐标分别为（x_i，y_i，z_i））；②再将上一步得到的屏幕坐标使用函数gluUnprotect()投影至近裁剪面z_p，设坐标分别为（xp_i，yp_i，z_p）。这样便实现了用于拾取局部区域的圆环和软组织投影点的坐标统一。然后遍历所有的点（xp_i，yp_i，z_p），并根据软组织投影点和圆的位置关系，提取出位于圆内的点集和三角形集合，分别设为P和T； 

gluprotect()函数的说明： 

int gluProtect(GLdouble objx，GLdouble objy，GLdouble objz，const GLdouble 

modelMatrix[16],const GLdouble projMatrix[16]，const Glint viewport[4], 

GLdouble*winx,GLdouble*winy，GLdouble*winz)； 

这个函数使用由模型视图矩阵（modelMatrix）、投影矩阵（projMatrix）和视口（viewport）所定义的变换，把指定的物体坐标（objx，objy，objz）变换为窗口坐标，并在winx、winy、winz中返回。 

其中投影点和圆的位置关系可以用简单的条件表达式来描述：(xp_i-x_new)²+(yp_i-y_new)²+(z_p-z_p)²≤0?圆内：圆外。 

（4）根据集合T中三角形的网格拓扑关系，从顶点集合P中分离出边界点集，设为集合B； 

具体分离办法如下： 

本发明中为了便于检索顶点、三角型网格彼此之间的拓扑关系，采用了改进的三角网格数据结构，下面将对其进行详细的说明。 

typedefvector<ZJGVertex*>vvector； 

typedefvector<ZJGEdge*>evector； 

typedefvector<ZJGFace*>fvector； 

①自定义的顶点结构 

②自定义的边结构 

③自定义的面结构 

判断连续三角网格中的某一点P是否为边界点，可通过点P的邻接点集合V_p来判断。若V_p中的点都能通过三角网格的边相连组成闭合的曲线，则点P为内点，否则为边界点。如图9所示，V_H={A，I，J，E，F，G}，因这些点通过三角形的边AI，IJ，JE，EF，FG，GA组成闭合的曲线，所以点H为内点，而不是边界点；而V_A=（B，I，H，G），因B和G不相邻，所以点A为边界点。 

在实际应用中，由于空间曲面上的点分布不规则，对点进行排序比较困难，因此用点的 邻接三角形网格集合来对点进行判断，基本步骤如下： 

建立一个临时点数组Temp，并置为空。 

对当前点的所有邻接三角网格上异于当前点的点作如下操作：如果Temp包含此点，则从Temp中删除此点，否则，把此点加入到Temp中。 

检查Temp，如果Temp为空，则当前点为内点，否则为边界点，并且Temp中的点是和当前点的邻接边界点。 

（5）最后对位于拾取区域内的三角形集合T，根据顶点之间的连接关系，实时生成局部质点——弹簧模型，并将边界点集合B的顶点固定，如图5和图6所示。 

具体生成办法如下： 

①遍历顶点集合P，把P中的所有顶点添加至质点弹簧系统中。 

②遍历三角形集合T中的所有边，若当前边i尚未生成弹簧模型，则对该边添加弹簧，并设置当前边是否已加弹簧的标志位为true。 

③遍历边界点集合B，并将B中的所有顶点固定。 

循环（2）—（5）步操作，便可以实现虚拟手术系统中局部质点——弹簧模型的实时动态生成。