基于GPU硬件加速的GraphCuts三维图像分割方法

摘要

本发明涉及一种图像处理技术领域的基于GPU硬件加速的Graph Cuts三维图像分割方法，包括如下步骤：首先，将原始三维图像载入为GPU中的三维纹理，并根据原始图像，调用初始化着色器程序来建立初始图结构，在GPU中生成相应的纹理数据；然后，对纹理数据中的图结构，在GPU着色器程序中进行推进-重标号操作实现最大流最小割的迭代计算，并使用GPU中遮挡查询技术作为迭代结束条件；最后，根据最小割的结果得到分割后的前景区域，并在对应像素上作标记。本发明在速度上有了明显的提升，适用于三维图像，所提出的Graph Cuts中图结构在GPU中的表示方式，也可以在Graph Cuts的其他应用领域中使用。

权利要求

1.一种基于GPU硬件加速的Graph Cuts三维图像分割方法，其特征在于，包括如下步骤：首先，将原始三维图像载入为GPU中的三维纹理，并根据原始图像，调用初始化着色器程序来建立初始图结构，在GPU中生成相应的纹理数据；

然后，对纹理数据中的图结构，在GPU着色器程序中进行推进-重标号操作实现最大流最小割的迭代计算，并使用GPU中遮挡查询技术作为迭代结束条件；

所述使用GPU中遮挡查询技术作为迭代结束条件，具体是指，每次迭代过程中，首先将深度缓存中所有像素的深度值置为0.5，然后对仍需进行迭代计算的像素，写入深度值0，使其通过GPU的深度测试；而对不需进行迭代计算的像素，写入深度值1；再利用GPU遮挡查询获得通过深度测试的像素总数，当总数为零时则结束迭代；

最后，根据最小割的结果得到分割后的前景区域，并在对应像素上作标记。

2.根据权利要求1所述的基于GPU硬件加速的Graph Cuts三维图像分割方法，其特征在于，所述初始图结构，是存储于GPU中的三维纹理，具体是根据原始图像生成Graph Cuts算法所需要的网络图数据结构，其中需要设定的参量包括每个像素与源点和汇点的连接容量，以及每个像素与x轴、y轴、z轴共6个方向上的相邻像素的连接容量，初始盈余，以及初始高度函数值。

3.根据权利要求1所述的基于GPU硬件加速的Graph Cuts三维图像分割方法，其特征在于，所述初始化着色器程序，是GPU中的一段片元着色器程序，它以对应于原始三维图像的纹理为输入，输出表示图结构的三个三维纹理，纹理中的数据即表示初始图结构的参量，它们的组织方式如下表：

4.根据权利要求1所述的基于GPU硬件加速的Graph Cuts三维图像分割方法，其特征在于，所述初始图结构和初始化着色器程序中，需要设定的参量的赋值方式如下：像素p与源点s的连接容量可按如下公式计算：cs，p＝-lnPr(Ip|″bkg″)

与汇点t的连接容量可按如下公式计算：

cp，t＝-lnPr(lp|″obj″)

相邻像素p，q之间的连接容量为：

其中：cp，q为图结构中(p，q)边的容量，Ip为原始图像中像素p的颜色(在单颜色通道图像中即为灰度)，||Ip-Iq||为两像素颜色差的模，||p-q||为两像素间的距离，Pr(Ip|″bkg″)和Pr(Ip|″obj″)为像素颜色Ip分别作为背景或前景时的概率，它们的值是通过对用户交互时预先设定的前景和背景像素统计得出；σ参数的值可取为全图像中相邻像素的颜色差的期望，即：σ2＝E(||Ip-Iq||2)

盈余表示每个像素上入流与出流的差，初始时为0；高度函数表示像素在图中与汇点间的最短距离，初始时为1。

5.根据权利要求1所述的基于GPU硬件加速的Graph Cuts三维图像分割方法，其特征在于，所述推进-重标号操作的GPU实现，具体为利用专用的GPU片元着色器 程序，并发的对每个像素作推进或重标号操作，其中，推进操作为：首先从源点向像素推流，以充满源点和像素间的连接边，并更新像素的盈余，然后若像素盈余不为0，则从当前像素向汇点推进，若像素盈余仍不为0，则将盈余向邻近像素推进；若像素在当前迭代中没有推进操作，则启动重标号操作：将像素的高度数值更新为与其连通的邻域像素的高度函数的最小值。

6.根据权利要求1所述的基于GPU硬件加速的Graph Cuts三维图像分割方法，其特征在于，所述根据最小割的结果得到分割后的前景区域，具体是指，将迭代结束后高度函数值大于设定阈值或盈余大于零的像素归为分割出的前景区域内的像素。 

说明书

技术领域

本发明涉及一种图像处理技术领域的三维图像分割方法，具体是一种基于GPU硬件加速的Graph Cuts三维图像分割方法。 

背景技术

图像分割是数字图像处理领域的重要技术，是计算机视觉技术的基础。以往对二维图像分割的研究成果相对较多，而三维图像在医学、工业制造、地质勘探等领域也均有重要应用，因而对三维图像分割的研究，其意义同样重大。Graph Cuts是近年来备受关注的图像分割框架之一，具有很好的精度和灵活性，在实际应用中倍受青睐。Graph Cuts的缺陷在于，由于它将图像分割转化为流网络的最大流最小割问题来解决，所以计算量较大，对于像素数更多的三维图像尤为明显。 
经对现有技术文献的检索发现，国外有Nandan Dixit等人(Nandan Dixit，Renaud Keriven，Nikos Paragios，“GPU-Cuts：Combinatorial Optimisation，Graphic Processing Units and Adaptive Object Extraction(GPU图割：组合优化，图形处理单元与对象分割)”，法国信息技术与系统教学与研究中心CERTIS研究报告，2005)提出的利用GPU加速的Graph Cuts图像分割方法，但这种方法只能对二维图像进行分割，而对三维图像的分割无能为力。国内的程广斌等人(程广斌，马承华，郝立巍，“基于图形加速器加速的医学图像分割研究”，医疗卫生装备第29卷第2期，2008)利用GPU实现了基于吉布斯随机场(Gibbs random field)的医学图像分割，吉布斯随机场的思想与Graph Cuts类似，但他们的方法仍然局限在二维分割上，无论其数据存储方式还是程序实现，都无法扩展到三维图像分割上。 

发明内容

本发明的目的在于针对现有技术的不足，提供一种利用GPU硬件加速的GraphCuts三维图像分割方法，利用GPU具有高速的实数运算特性和流式并行架构，弥补现有技术在三维图像支持性和计算效率上的缺陷。
本发明是通过以下技术方案实现的： 
本发明基于Graph Cuts分割框架精度高、易扩展的优点，在GPU中实现了三维图像分割。本发明利用GPU硬件的高速并行浮点运算特性，通过GPU中的三维纹理作为数据存储介质，将GPU分割的应用范围从二维扩展到三维。其原理是，对原始图像建立图结构，并将其存储在GPU的三维纹理数据中，然后在GPU中对图结构求取最小割。 
本发明包括如下步骤： 
首先，将原始三维图像载入为GPU中的三维纹理，并根据原始图像，调用初始化着色器程序来建立初始图结构，在GPU中生成相应的纹理数据； 
然后，对纹理数据中的图结构，在GPU着色器程序中进行推进-重标号操作实现最大流最小割的迭代计算，并使用GPU中遮挡查询技术作为迭代结束条件； 
最后，根据最小割的结果得到分割后的前景区域，并在对应像素上作标记。 
所述的初始图结构，是存储于GPU中的三维纹理，具体是根据原始图像生成Graph Cuts算法所需要的网络图数据结构，并对这三个纹理进行赋值。其中需要设定的参量包括每个像素与源点和汇点的连接容量，以及每个像素与x轴、y轴、z轴共6个方向上的相邻像素的连接容量，初始盈余，以及初始高度函数值。其中，每个像素p与源点s的连接容量按下式计算： 
cs，p＝-lnPr(Ip|″bkg″) 
每个像素p与汇点t的连接容量按下式计算： 
cp，t＝-lnPr(Ip|″obj″) 
相邻像素p，q之间的连接容量为： 
$c_{p,q}=\left\{\begin{array}{cc}1,& I_p\le I_q\\ \exp (-\frac{{\left|\right|I_p-I_q\left|\right|}^2}{{2\sigma }^2})\frac{1}{\left|\right|p-q\left|\right|},& I_p>I_q\end{array}\right.$
其中，cp，q为图结构中(p，q)边的容量，Ip为原始图像中像素p的颜色(在单颜 色通道图像中即为灰度)，‖Ip-Iq‖为两像素颜色差的模，‖p-q‖为两像素间的距离。Pr(Ip|″bkg″)和Pr(Ip|″obj″)为像素颜色Ip分别作为背景或前景时的概率，它们的值是通过对用户交互时预先设定的前景和背景像素统计得出；σ参数的值可取为全图像中相邻像素的颜色差的期望，即： 
σ2＝E(‖Ip-Iq‖2) 
盈余表示每个像素上入流与出流的差，初始时为0；高度函数表示像素在图中与汇点间的最短距离，初始时设为1。图结构在GPU中以三个纹理数据的方式进行存储，其中的数据组织方式见附表1。 
所述的初始化着色器程序，是指GPU中的一段片元着色器程序。它以对应于原始三维图像的纹理为输入，输出表示图结构的三个三维纹理，并对这三个纹理进行赋值。并按照初始图结构中提到的赋值方式向这三个纹理写入数据。 
所述的推进-重标号算法的GPU实现，具体为利用GPU片元着色器程序，并发地对每个像素作推进或重标号操作。其中，推进操作为：首先从源点向像素推流，以充满源点和像素间的连接边，并更新像素的盈余；然后若像素盈余不为0，则从当前像素向汇点推进，若像素盈余仍不为0，则将盈余向邻近像素推进。若像素在当前迭代中没有推进操作，则启动重标号操作：将像素的高度函数值更新为： 
$h\left(p\right)=1+\underset{B_{p,q}>0}{\min }\left\{h\left(q\right)\right\}$
所述的使用GPU中遮挡查询技术作为迭代结束条件，具体是指，每次迭代过程中，首先将深度缓存中所有像素的深度值置为0.5，然后对仍需进行迭代计算的像素，写入深度值0，使其通过GPU的深度测试；而对不需进行迭代计算的像素，写入深度值1。再利用GPU遮挡查询获得通过深度测试的像素总数，当总数为零时则结束迭代。 
所述的根据最小割的结果得到分割后的前景区域，具体是指，将迭代结束后高度函数值大于设定阈值或盈余大于零的像素归为分割出的前景区域内的像素。 
本发明基于Graph Cuts分割框架精度高，支持三维分割的优点，利用GPU硬件的高速并行浮点运算特性，弥补了Graph Cuts三维分割时间效率低的问题，同时也是对现有Graph Cuts的GPU实现方法向三维上的扩展和改进，所提出的Graph Cuts中图结构在GPU中的表示方式，也可以在Graph Cuts的其他应用领域中使用。 

附图说明

图1为对一幅图像建立图结构的示意图； 
其中：s为虚拟添加的源点，t为虚拟添加的汇点，灰色圆代表图像像素；源点到像素，像素到汇点，以及相邻像素间有边相连，边的容量分别为cs，p，cp，t以及cp，q；虚线代表图的一个割，它将图像中的像素分为两部分，分别为前景像素和背景像素。 
图2为基于GPU加速的Graph Cuts分割方法的流程图； 
其中，活动像素是指盈余大于0且高度函数值小于一个阈值的像素。 
图3为分割前后的对其中某一片图像的二维显示效果图； 
其中(a)为分割前的图像，(b)为分割图像。 
图4为使用本方法对盆腔中的股骨头进行分割前后的三维绘制效果图； 
其中(a)为分割前的绘制效果，(b)为分割后对前景和背景部分分别进行绘制再融合后的效果。 

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。 
本实施例首先将原始三维图像载入到GPU中的三维纹理texImg，并根据该纹理生成表示图结构的三个纹理texLink1，texLink2和texTerm，然后使用着色器程序进行推进-重标号操作来迭代计算图的最小割，迭代结束条件是图中没有活动像素，最后根据图的最小割得到图像分割结果，并对分割出的像素写入标记。 
本实施例包括如下步骤： 
1.数据初始化。首先将原始三维图像载入为三维纹理，对于图3、4中的医学 图像使用LUMINANCE纹理格式；对于彩色图像，使用RGB纹理格式；然后调用用于初始化操作的GPU片元着色器程序来生成用于存储图结构的纹理数据，该程序的输入为texImg纹理，输出为texLink1、texLink2、texTerm三个纹理数据。三个输出纹理均为RGB格式，三个分量的含义如表1所示。该着色器程序并发地在图像的每个像素上执行，执行内容为设置三个输出纹理的各个分量上的数值。其中，texLink0和texLink1存放的都是朝向相邻像素的边的容量，其值可按如下公式计算： 
$c_{p,q}=\left\{\begin{array}{cc}1,& I_p\le I_q\\ \exp (-\frac{{(I_p-I_q)}^2}{2{\sigma }^2})\frac{1}{\left|\right|p-q\left|\right|},& I_p>I_q\end{array}\right.$
其中，p为当前像素，q为相邻像素，Ip为像素p的颜色，‖p-q‖为像素间的间距，σ2为一个可设置的参数，在本例中，将其设为全图像中相邻像素的颜色差的期望，即： 
σ2＝E(‖Ip-Iq‖2) 
其中E(S)表示S的数学期望。 
将来自源点的流的部分或全部推向汇点，可以将源点边容量与汇点边容量合并为一个量，即它们的差，texTerm中只用一个分量R来表示，其值可按如下公式计算： 
$\mathrm{term}.R=\ln \frac{\mathrm{Pr}\left(I_p\right|{}^{\prime \prime }\mathrm{obj}^{\prime \prime })}{\mathrm{Pr}\left(I_p\right|{}^{\prime \prime }\mathrm{bkg}^{\prime \prime })}$
其中，Pr(Ip|″obj″)表示在对象中像素p的亮度Ip出现的概率，Pr(Ip|″bkg″)为在背景中亮度Ip出现的概率，在本实施例中，它们的值通过对一系列已分割图像进行统计获得。盈余为像素的入流与出流之差，初始时设为0。高度函数表示从当前像素开始通过容量不为0的边，到底汇点所需要经过的边的数目，初始时设为1。 
2.迭代进行推进-重标号操作。首先需要将深度缓存清空为1.0，然后调用着色器程序进行计算。用于推进-重标号操作的片元着色器程序，其输入和输出均为 texLink0，texLink1，texTerm三个纹理数据。该着色器程序并发的在图像的每个像素上执行，其作用是更新三个纹理各个分量上的数值。具体执行步骤为： 
(1)若texTerm.R＞0(texTerm.R表示texTerm纹理在当前像素对应坐标下的R分量，下同)，则从源点推进同等的流至当前像素，并修改盈余： 
texTerm.G+＝texTerm.R.texTerm.R＝0 
(2)处理从相邻像素推向当前像素的入流：检查所有相邻像素，若某相邻像素q满足texTerm(q).G＞0，texTerm(q).B＝texTerm.B+1，且texLink(q，p)＞0(texLink(q，p)为从q到p的边的容量，可根据q与p的相对位置从texLink1或texLink2中的对应分量获取)，则可从像素q向当前像素推进流，对应的更新方式为：texTerm.G+＝min(texTerm(q).G，texLink(q，p))。 
(3)若texTerm.R＜0且texTerm.G＞0，即指向汇点的边的容量未被充满，且盈余大于0，则向汇点推进流： 
flow＝min(-texTerm.R，texTerm.G)， 
texTerm.R+＝flow，texTerm.G-＝flow，texTerm.B＝1 
(4)从当前像素向相邻像素推进流：检查所有相邻像素，若某相邻像素q满足texTerm(q).B＝texTerm.B-1，texLink(p，q)＞0，则可从当前像素向像素q推进流，对应的更新方式为： 
fl＝min(texTerm.G，texLink(p，q)) 
texTerm.G-＝fl，texLink(p，q)-＝fl 
(5)若上述(1)到(4)中，没有任何一项操作对当前像素实施，则进行重标号操作：在相邻像素中，找到一个连接边大于0(texLink(p，q)＞0)且高度函数值最小的像素qmin，然后更新当前像素的高度函数值为： 
texTerm.B＝texTerm(qmin).B+1 
(6)经过(1)到(5)的操作后，若当前像素满足活动像素条件，即 
texTerm.G＞0且texTerm.B＜th1 
则输出当前像素的深度值为0，否则输出为1(上式中th1是一个可设定的阈值参数，在本例中设为三维图像在xyz三方向上的像素数的最大值)。 
在每次迭代完成后，利用GPU的遮挡查询技术获取通过深度测试的像素数N。 若N＞0，则继续迭代，否则完成整个推进-重标号迭代过程，转入下一步。 
3.输出分割标记。分割标记用一个与原始图像同样大小的三维纹理存放，且与原始图像上的像素一一对应，将它表示为texTag。这一步所需的着色器程序，其输入为texTerm，输出为texTag。它并发的在每一个像素上执行，执行内容为更新texTag的数据，具体为：若texTerm.G＞0或texTerm.B＞＝th1，则当前像素为活动像素，可分割为前景区域，于是置texTag为1，否则置为0。更新完成后，texTag可用于后续的体绘制过程。 
图3和图4分别为本例的分割结果的展示，其分割对象为人体盆腔中的股骨头。其中图3为三维图像的其中一片上的分割效果，而从图4中可以看出在分割前后三维体绘制效果的差异。表2中显示了对三组不同规模的数据，分别利用GPU和CPU实现的三维Graph Cuts分割所需耗费的执行时间。 
表1三维Graph Cuts的GPU实现中的纹理数据组织方式 

表2基于本发明的程序的执行时间与传统CPU实现方法的对比 
    数据大小   CPU实现方法的 执行时间(秒)    本发明的  执行时间(秒)    数据组1    512×512×76   30.890    12.219    数据组2    512×512×39   14.797    5.312    数据组3    512×512×68   30.703    11.297