本发明涉及一种对体数据的分割结果进行人机交互编辑的方法，属于数字化医疗技术领域。
背景技术：
在数字化医疗的过程中，通过ct、mri等诊疗手段获取病人的体数据，需要通过体数据分割技术获取目标区域，有助于进行进一步的分析与诊断。由于医学体数据的复杂性，现有的体数据自动分割技术普适性较差，分割结果通常存在问题，需要使用人机交互技术对分割结果进行进一步的编辑。然而对于三维体数据进行数字化的人机交互编辑是较为困难的，现有的技术不能较好地解决三维人机交互面临的自由度高、交互精度差等问题。为了高效地对体数据分割结果进行编辑，需要设计良好的交互方法。已有的seg3d等软件采用的交互体数据分割方法，该方法中交互方式为，用户在体数据的二维切片上使用笔刷标记出分割目标。为了完成分割任务，需要标记大量二维切片，交互效率较低。或借助软件提供的自动分割方法，得到的分割结果可能存在一定的误差。已有的基于网格编辑的交互体数据分割方法，将分割结果以三角网格模型的形式进行表示，通过用户交互构建三角网格模型完成分割目的，或交互编辑网格形状来对分割结果进行编辑。该方法在网格形状计算过程中没有充分考虑体数据信息，对用户交互的依赖较强，可能需要额外的交互次数来保证得到较好的分割结果。已有的turtleseg等软件采用的基于轮廓编辑的体数据分割方法，该方法对体数据选取一个平面切片，用户在切片上标出分割目标的轮廓线，作为网格生成的约束条件。通过对多条轮廓线标记得到的结果进行曲面重建得到分割结果的网格模型。该方法在精度及效率上存在不足，需要标记较多的轮廓以保证分割结果细节处的精度。相关方法可参见ijirit,yokotah.contour-basedinterfaceforrefiningvolumesegmentation[j].computergraphicsforum,2010,29(7),2153-2160.
现有的基于二维切片交互的体数据交互分割方法需要操作者对大量切片进行标记以得到较好的结果，交互效率低，分割一份数据需要的时间可能长达数小时。基于网格的分割方法中，算法对体数据信息的利用率较低，严重依赖用户交互，由于缺少体数据作为约束，仅考虑网格形状，分割精度较低。基于二维轮廓线的曲面重建方法对于细微结构难以精确分割，同时需要较多的用户交互次数。
技术实现要素：
本发明的目的是提出一种对体数据的分割结果进行人机交互编辑的方法，以克服现有技术中不能结合三维交互与二维交互优势，以及不能充分利用体数据信息的技术问题，进而提升交互效率与分割精度。
本发明提出的对体数据的分割结果进行人机交互编辑的方法，包括以下步骤：
(1)从计算机中读取以二值数据形式存储的体数据的分割结果，使用立方体步进法方法，对分割结果进行等值面提取，得到分割结果的三角网格模型m0，三角网格模型m0表示分割目标的形状和体积特征，立方体步进方法中的提取等值面的目标值为0.5，将上述三角网格模型m0以顶点坐标和三角面片索引的方式进行存储；采用图形处理单元光栅化方法对该三角网格模型m0进行面绘制，并记录三角网格中各面片所在位置的体数据坐标；
(2)从计算机中读取与上述步骤(1)的分割结果相对应的体数据，根据观察者位置，使用光线投射方法，分别从与每个计算机屏幕像素对应的空间位置，对体数据进行体绘制，即沿视线方向，向体数据发出射线，得到投射光线，根据上述步骤(1)的三角网格中各面片所在位置的体数据坐标，对投射到体数据上的光线进行判断，若光线投射到上述步骤(1)中的三角网格面片上，则结束光线投射，若光线未投射到面片上，则投射光线与体数据相交，得到相交线段，在该相交线段上等距获取多个采样点，读取采样点处的体数据灰度值，沿光线投射方向进行体数据灰度值的累加，得到投射光线在体数据中传播过程中积累的亮度值，根据该亮度值得到每个屏幕像素的颜色值，得到体绘制的二维图像；
(3)在同一视图中，采用透明度混合方法，对上述步骤(1)的面绘制二维图像和上述步骤(2)的体绘制二维图像进行混合，得到统一显示图像；
(4)在上述步骤(3)的统一显示图像中，记录用户进行鼠标交互的笔触线，将该笔触线以屏幕像素的坐标点序列形式进行存储，笔触线中每两个坐标点依次相连形成折线段，使用上述步骤(2)中的光线投射方法，沿观察者视线方向，从该折线段上与所有点相对应的屏幕像素，向屏幕内部方向延伸得到多个平行射线，多个平行射线与上述步骤(2)的体数据相交得到多个平行线段，相邻两条平行线段组成一个四边形，将该四边形根据上述笔触线坐标点序列顺序连接，得到一个延伸面，该延伸面以三角网格的形式进行存储，记为截面网格ms，使该截面网格ms与体数据相交，并根据该截面网格ms的网格顶点坐标得到相应的体数据坐标，进而得到体数据的灰度，并在该截面网格ms上显示体数据的灰度；
(5)在步骤(4)中得到的截面网格ms与步骤(1)的三角网格模型m0中各取一个三角面片，对与两个三角面片相对应的两个三角形进行空间相交测试，得到相交的两个三角形在空间相交的交线段的两个端点的三维空间坐标和重心坐标，遍历两个三角网格中所有的三角面片选取方式，得到上述两个三角网格中三角面片相交的所有交线段，将上述所有交线段按照空间相邻关系连接成连续的折线段，记为截面网格ms与三角网格模型m0的相交线c0；
(6)将上述步骤(4)截面网格ms中的多个依次相连的四边形，按照步骤(4)中的笔触线坐标点序列顺序排列在一个二维平面上，得到二维平面上的截面网格ps，采用图形处理单元光栅化方法，对该二维截面网格ps进行面绘制，得到体数据灰度在截面网格ps的显示结果图像；根据上述截面网格ps，设定上述步骤(5)的相交线c0在截面网格ps中的三角面片上的重心坐标不变，得到二维平面上的相交线cp；使上述二维截面网格ps与二维相交线cp共同显示在同一视图中，将该视图记为二维视图；
(7)对上述步骤(6)的二维视图中二维截面网格上的体数据灰度的视觉边界与该二维视图中的二维相交线cp所在的位置进行对比，当相交线偏离上述步骤(1)的分割目标的边界时，使用鼠标移动相交线上的点，将相交线上的点移动到分割目标的边界上，得到相交线cn，并得到相交线cn上的点的重心坐标；
(8)根据上述步骤(7)的相交线cn上的点的重心坐标以及上述步骤(4)中的截面网格ms的网格顶点坐标，得到相交线cn上点的三维空间坐标c；
(9)建立基于体数据梯度的三角网格优化模型，该三角网格优化模型的目标函数e为：
e＝‖l(v′)-t(v′)δ0‖2+ω‖mv′-c‖2-κg(v′)
其中，v′为目标函数的自变量，即优化后的三角网格模型m0的顶点坐标，通过最小化目标函数e，求解得到v′，
‖l(v′)-t(v′)δ0‖2为拉普拉斯形状约束，其中l为拉普拉斯坐标算符，
l(v′)为优化后三角网格模型m0的拉普拉斯坐标，δ0为初始网格m0的拉普拉斯坐标，t(v′)为用顶点坐标线性组合表示的顶点局部变换；
ω‖mv′-c‖2为用户交互约束，ω为权重参数，c为上述步骤(8)中得到的相交线cn上点的三维空间坐标c，m为参数矩阵，即将相交线cn上点的三维空间坐标c表示为三角网格模型m0顶点坐标的线性组合，并将该线性组合表示成矩阵形式，即可得到参数矩阵m；
-κg(v′)为体数据梯度约束，κ为该项的权重参数，g为体数据梯度函数，g(v′)表示顶点v′处的体数据梯度大小，该梯度函数g使用差分进行数值计算；
(10)将上述步骤(9)的优化函数e改写为e′：
e＇＝‖l(v′)-t(v′)δ0‖2+ω‖mv′-c‖2-κg(p)+θ‖v＇-p‖2
其中，p为辅助变量，表示优化过程中三角网格模型m0的顶点坐标，θ为松弛变量，表示‖v′-p‖2的权重系数，在优化过程中逐渐增大；
(11)提取上述步骤(10)中e′的以v′为变量的相关项，得到优化子问题s1，该优化子问题s1为已知辅助变量p，求解顶点坐标v′，其优化函数如下：
上述优化子问题s1为线性方程，利用线性最小二乘方法求解该线性方程，得到v′，并更新v′；
(12)提取上述步骤(10)中e′的以p为变量的相关项，得到优化子问题s2，该优化子问题s2为已知顶点坐标v′，求解辅助变量p，其优化函数如下：
上述优化子问题s2为梯度优化问题，该目标函数关于每个顶点独立，对于三角网格m0的每个顶点，利用邻域搜索的方法，求解该梯度优化问题，遍历三角网格m0的顶点，得到p，并更新p；
(13)设定松弛变量θ的最大值为θmax，初始化时，令松弛变量θ的初始值为θ0，取v′＝v0，其中v0为上述步骤(1)的三角网格m0的顶点坐标，随后开始进行迭代：
(13-1)令p＝v′，遍历三角网格m0的顶点，进行邻域搜索，求解上述优化子问题s2，得到并更新p；
(13-2)根据上述步骤(13-1)的p，使用线性最小二乘方法，求解上述优化子问题s1，得到并更新v′；
增大θ值，重复上述步骤(13-1)和(13-2)，对松弛变量θ进行判断，若θ≥θmax，则停止迭代，进行步骤(14)，若θ＜θmax，则返回(13-1)；
(14)显示与上述步骤(13)的顶点坐标v′对应的三角网格，该三角网格即为对体数据分割结果进行一次人机交互编辑得到的结果。
本发明提出的对体数据的分割结果进行人机交互编辑的方法，其优点是：
本发明方法，通过构建体绘制与面绘制结合且统一显示的交互环境，便于识别出分割目标；使用基于笔触绘制的切片生成方法，对分割结果进行相交线求取，能够适应不同形状的分割目标，该方法有利于更加直观地确认所编辑分割结果中存在的问题；将三维切片平面化生成二维切片，有利于更加准确地进行交互，提高交互编辑精度；由基于体数据梯度的网格优化算法将二维编辑结果传播至三维，节省了用户交互的次数，同时将分割网格边界靠近体数据边界，提升了分割精度。
具体实施方式
本发明提出的对体数据的分割结果进行人机交互编辑的方法，包括以下步骤：
(1)从计算机中读取以二值数据形式存储的体数据的分割结果，使用立方体步进(marchingcubes)方法，对分割结果进行等值面提取，得到分割结果的三角网格模型m0，三角网格模型m0表示分割目标的形状和体积特征，立方体步进方法中的提取等值面的目标值为0.5，该过程中，输入的体数据分割结果以二值数据形式进行存储，在每个位置上，使用0或1表示该位置是否属于分割目标，若为0表示该位置不属于分割目标，若为1表示该位置属于分割目标，所有1位置组成的三维对象即为该分割数据表示的分割结果，因此立方体步进方法中的提取等值面的目标值设定为0.5。将上述三角网格模型m0以顶点坐标和三角面片索引的方式进行存储；采用图形处理单元(gpu)光栅化方法对该三角网格模型m0进行面绘制，此面绘制过程使用gpu设备的深度缓存进行遮挡剔除，得到正确的面绘制结果。另外，在面绘制过程中，记录三角网格中各面片所在位置的体数据坐标；上述三角网格模型m0的面绘制结果是二维图像，该图像可在屏幕上进行显示。
(2)从计算机中读取与上述步骤(1)的分割结果相对应的体数据，根据观察者位置，使用光线投射方法，分别从与每个计算机屏幕像素对应的空间位置，对体数据进行体绘制，即沿视线方向，向体数据发出射线，得到投射光线，根据上述步骤(1)的三角网格中各面片所在位置的体数据坐标，对投射到体数据上的光线进行判断，若光线投射到上述步骤(1)中的三角网格面片上，则结束光线投射，若光线未投射到面片上，则投射光线与体数据相交，得到相交线段，在该相交线段上等距获取多个采样点，读取采样点处的体数据灰度值，沿光线投射方向进行体数据灰度值的累加，得到投射光线在体数据中传播过程中积累的亮度值，根据该亮度值得到每个屏幕像素的颜色值，得到体绘制的二维图像；该二维图像可在屏幕上进行显示。
(3)在同一视图中，采用透明度混合方法，对上述步骤(1)的面绘制二维图像和上述步骤(2)的体绘制二维图像进行混合，得到统一显示图像；
(4)在上述步骤(3)的统一显示图像中，记录用户进行鼠标交互的笔触线，将该笔触线以屏幕像素的坐标点序列形式进行存储，笔触线中每两个坐标点依次相连形成折线段，该折线段表示通过观察上述统一显示的图像，找出分割结果中需要编辑的部位，依据待编辑部位的形状，对三维体数据进行截取时的截面形状。使用上述步骤(2)中的光线投射方法，沿观察者视线方向，从该折线段上与所有点相对应的屏幕像素，向屏幕内部方向延伸得到多个平行射线，多个平行射线与上述步骤(2)的体数据相交得到多个平行线段，相邻两条平行线段组成一个四边形，将该四边形根据上述笔触线坐标点序列顺序连接，得到一个延伸面，该延伸面以三角网格的形式进行存储，记为截面网格ms，使该截面网格ms与体数据相交，并根据该截面网格ms的网格顶点坐标得到相应的体数据坐标，进而得到体数据的灰度，并在该截面网格ms上显示体数据的灰度；使得能够通过观察截面网格了解体数据的灰度情况，并使得通过灰度能够分辨出分割结果对应的分割目标。
(5)在步骤(4)中得到的截面网格ms与步骤(1)的三角网格模型m0中各取一个三角面片，对与两个三角面片相对应的两个三角形进行空间相交测试，得到相交的两个三角形在空间相交的交线段的两个端点的三维空间坐标和重心坐标，遍历两个三角网格中所有的三角面片选取方式，得到上述两个三角网格中三角面片相交的所有交线段，将上述所有交线段按照空间相邻关系连接成连续的折线段，记为截面网格ms与三角网格模型m0的相交线c0；
(6)将上述步骤(4)截面网格ms中的多个依次相连的四边形，按照步骤(4)中的笔触线坐标点序列顺序排列在一个二维平面上，得到二维平面上的截面网格ps，采用图形处理单元光栅化方法，对该二维截面网格ps进行面绘制，得到体数据灰度在截面网格ps的显示结果图像；根据上述截面网格ps，设定上述步骤(5)的相交线c0在截面网格ps中的三角面片上的重心坐标不变，得到二维平面上的相交线cp；使上述二维截面网格ps与二维相交线cp共同显示在同一视图中，将该视图记为二维视图；
(7)对上述步骤(6)的二维视图中二维截面网格上的体数据灰度的视觉边界与该二维视图中的二维相交线cp所在的位置进行对比，当相交线偏离上述步骤(1)的分割目标的边界时，使用鼠标移动相交线上的点，将相交线上的点移动到分割目标的边界上，得到相交线cn，并得到相交线cn上的点的重心坐标；
(8)将上述步骤(7)的二维相交线cp到新相交线cn的形变，转换为三维网格顶点的位置变换。方法为根据上述步骤(7)的相交线cn上的点的重心坐标以及上述步骤(4)中的截面网格ms的网格顶点坐标，得到相交线cn上点的三维空间坐标c；
(9)建立基于体数据梯度的三角网格优化模型，用于将步骤(7)中的二维交互编辑结果传播至三维，使得二维交互编辑操作能够带动三维网格m0变形，从而完成对体数据分割结果的编辑。该三角网格优化模型的目标函数e为：
e＝‖l(v′)-t(v′)δ0‖2+ω‖mv′-c‖2-κg(v′)
其中，v′为目标函数的自变量，即优化后的三角网格模型m0的顶点坐标，通过最小化目标函数e，求解得到v′，
‖l(v′)-t(v′)δ0‖2为拉普拉斯形状约束，其中l为拉普拉斯坐标算符，
l(v′)为优化后三角网格模型m0的拉普拉斯坐标，δ0为初始网格m0的拉普拉斯坐标，t(v′)为用顶点坐标线性组合表示的顶点局部变换；以解决拉普拉斯坐标对旋转缩放变换敏感的问题。
ω‖mv′-c‖2为用户交互约束，ω为权重参数，其取值根据经验进行调整，本发明的一个实施例中取为0.5，c为上述步骤(8)中得到的相交线cn上点的三维空间坐标c，m为参数矩阵，即将相交线cn上点的三维空间坐标c表示为三角网格模型m0顶点坐标的线性组合，并将该线性组合表示成矩阵形式，即可得到参数矩阵m；
-κg(v′)为体数据梯度约束，κ为该项的权重参数，其取值根据经验进行调整，本发明的一个实施例中取为1.0，g为体数据梯度函数，g(v′)表示顶点v′处的体数据梯度大小，该梯度函数g使用差分进行数值计算；
(10)将上述步骤(9)的优化函数e改写为e′：
e′＝‖l(v′)-t(v′)δ0‖2+ω‖mv′-c‖2-κg(p)+θ‖v′-p‖2
其中，p为辅助变量，表示优化过程中三角网格模型m0的顶点坐标，θ为松弛变量，表示‖v′-p‖2的权重系数，在优化过程中逐渐增大；
上式中使用辅助变量p替代了非凸函数g(v′)中的v′，同时增加了‖v′-p‖2项进行约束，使得辅助变量p与被替代的v′尽可能接近。当θ→+∞时，有e′→e；
(11)提取上述步骤(10)中e′的以v′为变量的相关项，得到优化子问题s1，该优化子问题s1为已知辅助变量p，求解顶点坐标v′，其优化函数如下：
上述优化子问题s1为线性方程，利用线性最小二乘方法求解该线性方程，得到v′，并更新v′；
(12)提取上述步骤(10)中e′的以p为变量的相关项，得到优化子问题s2，该优化子问题s2为已知顶点坐标v′，求解辅助变量p，其优化函数如下：
上述优化子问题s2为梯度优化问题，该目标函数关于每个顶点独立，对于三角网格m0的每个顶点，利用邻域搜索的方法，求解该梯度优化问题，遍历三角网格m0的顶点，得到p，并更新p；上述邻域搜索方法为，在每个顶点附近的一定大小的区域内选取一些点，对上述优化子问题s2的目标函数进行计算，得到使目标函数值最小的顶点坐标；
(13)将步骤(11)中的优化子问题s1与步骤(12)中的优化子问题s2进行迭代。
设定松弛变量θ的最大值为θmax，初始化时，令松弛变量θ的初始值为θ0，取v′＝v0，其中v0为上述步骤(1)的三角网格m0的顶点坐标，随后开始进行迭代：
(13-1)令p＝v′，遍历三角网格m0的顶点，进行邻域搜索，求解上述优化子问题s2，得到并更新p；
(13-2)根据上述步骤(13-1)的p，使用线性最小二乘方法，求解上述优化子问题s1，得到并更新v′；
增大θ值，重复上述步骤(13-1)和(13-2)，对松弛变量θ进行判断，若θ≥θmax，则停止迭代，进行步骤(14)，若θ＜θmax，则返回(13-1)；
(14)显示与上述步骤(13)的顶点坐标v′对应的三角网格，该三角网格即为对体数据分割结果进行一次人机交互编辑得到的结果。