一种基于GPU加速的体数据非真实感绘制方法

摘要

本发明一种基于GPU加速的体数据非真实感绘制方法，该方法包括以下步骤：在多个视点下，对医学体数据进行直接体绘制，获取多个渲染结果，对二维的渲染结果提取特征线，定义全局能量最小化规则，计算三维特征线，利用卷积曲面方法，定义一标量场，对三维特征线进行聚类，并利用B样条曲线拟合，完成风格化绘制。本方法借助GPU强大的计算能力，大大提高了渲染速度和效率。

说明书

技术领域

本发明属于非真实渲染及体绘制技术领域，具体涉及基于GPU加速的体数据非真实感 绘制方法。

背景技术

非真实感绘制作为真实感绘制的有效补充，其目的不在于图形的真实性，而主要在于产 生一种与艺术家作品相近的艺术效果，表现出图形的艺术特质，它已在教育、娱乐、艺术等 领域得到了广泛应用。然而由于三维数据规模不断增大，处理方法日益复杂，非真实感体绘 制也面临着与传统体绘制同样的难题：图像质量与绘制速度的矛盾。而随着可编程图形硬件 的出现，基于GPU(Graphics Process Unit，图形处理单元)的非真实感绘制成为该领域的一个 研究热点。

在90年代，非真实感绘制研究的重点集中在对自然材质的模拟，如钢笔画效果、铅笔 画效果、油画效果、水彩画效果等。Winkenbach等基于三维空间，引入“笔划纹理”的概念， 用笔划纹理填充三维模型的表面以模拟钢笔画。通过省略细节提高线条画的可理解性;同时， 笔划参数的不同可以适应不同分辨率曲面的表现。Curtis等提出了一种自动生成水彩画的方 法。该方法利用Kubelka-Munk色彩混合方式获取相当真实的水彩效果，但是计算量较大。 Hertzmann等提出了一种从照片创建手绘效果的模型，通过使用不同大小和形状的笔画绘制 图像，以模拟油画效果。绘制是基于图层的，每一层选择不同的笔触，每一层在上一层的基 础上绘制。

从艺术的角度看，轮廓线也占据很重要的位置，针对轮廓线的研究一直是NPR领域的 热点。在制作卡通画中，为了丰富的表达模型的形状，常用轮廓线来勾勒模型的大致轮廓， 并用粗线来表现模型的重要特征，达到吸引观众注意力的目的。目前，很多研究致力于这方 面的工作，其中最基本的是轮廓线的计算。轮廓线计算一般是基于图形空间和图像空间的相 关算法。基于图像空间的算法一般有：基于Depth Map、Normal Map探测图像的轮廓线，背 面线框绘制。

基于Depth Map的算法的主要思想是对于属于不同物体的点，像素深度值变化大:而对 于邻近的点，属于同一个物体上的点，其像素深度值变化小。Decaudin等提出利用用向量映 像来探测图像的轮廓线。该方法假设物体在轮廓线附近的面法向发生变化，其中向量映像是 特殊绘制的图像，图像中的每一个像素存储的是对应物体表面的法向量。

基于边的遍历检测算法直接遍历所有的边，判断视线与面法向量乘积不小于0，该面为 朝前面，反之，该面为朝后面，该方法简单但效率低下，Gooch提出了通过高斯球来加速该 算法。基于面的遍历检测算法是Buchanan提出的利用Edge Buffer数据结构遍历获取轮廓线 的方法，该方法只需要知道每个面的组成边，在没有硬件加速的特殊情况下可以极大提高检 测效率。

获取轮廓线后，研究者考虑如何能使用更丰富的线条展示模型的形状，而不是单纯的利 用简单的直线。研究者提出了轮廓线的风格化绘制(Stylized Silhouette Rendering，SSR)，它 是非真实感绘制中一种重要的表现形式。通过SSR可以模拟各种材质和绘制手法。

Northrup等首先提出对轮廓线进行风格化的实现方法。首先，利用随机检测的算法找出 轮廓线，然后，利用图形图像相结合的方法筛选出可见的轮廓线，对于可见的轮廓线，作者 利用纹理映射的方法对特征线进行风格化。其主要方法是将特征边沿相邻两边的角平分线方 向将特征线扩展成四边形，将代表不同材质的二维纹理映射到四边形中。Isenberg也提出了 一个风格化轮廓线的算法。利用渲染时的z-buffer信息来判断轮廓线的可见性，其主要原理 为轮廓线出现在z-depth不连续的地方，达到良好的视觉效果。

体绘制技术是一种重要的科学可视化技术，采用该技术，不经过复杂的特征算法，就能 够交互地浏览数据内部结构特征。该项技术广泛应用于医学、气象和地质等领域。目前直接 体绘制算法主要有Ray Casting算法、Shear-Warp算法、足迹算法、Texture Mapping算法等。 LevoyM等将体数据递归均匀子分成一个八叉树结构，在光线投射算法中遍历八叉树来跳过 包含当前采样点的最大空区域，以此达到加速的效果。

所谓特征点的跟踪与拟合是指在体数据中得到相应的特征点后，需要根据这些特征点的 一些特性，比如结构张量，曲率等等，把相应的特征点拟合成能够表达体数据内相关结构信 息的特征线，这样的特征线能够清晰的表达了体数据的内部重要结构和特征，是用户关心的 重要区域。Lawson等提出的基于三角形区域重构的Bezier方法。针对曲线曲面相交问题普 遍存在计算机辅助设计与加工制造中，国内外很多学者研究了这个问题。Anderssion等从曲 线自相交情况的条件入手研究，给出了Bezier曲线避免自相交的条件。Tatiana等研究了二 维曲线造型中出现自相交情况如何去除的问题。

本发明属于非真实渲染及体绘制技术领域，具体涉及基于GPU加速的体数据非真实感 绘制方法。

发明内容

本发明要解决的技术问题是：提供一种对体数据进行非真实绘制的方法，并利用GPU 硬件的计算能力，提高了计算与绘制效率，该方法主要利用具有视觉意义的二维特征线，并 根据全局能量最小化规则，计算对应的三维特征线并完成矢量化及风格化绘制效果。

本发明解决上述技术问题的技术方案为：基于GPU加速的体数据非真实感绘制方法， 包括如下步骤：

1）二维特征线的提取

步骤（1）、从多个视点出发进行渲染，具体为在XYZ笛卡尔坐标系中在3个轴的正负 方向六个视点，对场景进行Ray Casting（光线投射）渲染，该方法是对每个像素从视点出发 投射一条光线，该光线在体数据内投射过程中，按照一定的步长进行多次采样，并计算采样 点相应的颜色和不透明度，直到光线穿出体数据，该过程可以在GPU上实现，利用CUDA （计算统一设备架构）的硬件加速，大大提高渲染速度。并利用渲染到纹理技术，把绘制结 果渲染到多个纹理中，得到多个二维渲染结果。

步骤（2）、对步骤（1）的二维效果图利用图像处理相关技术，主要是Canny算子和Contour  Detect算法对二维图进行特征线的提取。Canny算子也是一阶算子，实质是用一个准高斯函 数做平滑运算，然后用带方向的一阶微分算子来定位导数最大值，其贡献在于给出了一种寻 找最优算子的思路，即确立了最优准则。Contour Detect主要是对由Canny Detect的计算结 果来计算点的连续性，这样就能在二维上得到特征线。

2）二维特征线向三维线的映射

步骤（3）、二维特征线在二维上是具有明显视觉意义的线，我们认为这些线是三维的线 在二维的投影效果，所以需要找到每条二维特征线对应三维的特征线。然而一个二维点是很 多个三维点的投影，所以问题转化为一个全局最小化问题。一个三维点有位置信息，不透明 度信息，颜色信息，这些是在最小化问题需要考虑的因子。

步骤（4）、根据步骤（3）定义的最小化问题，问题可以转化为类似于多阶段多节点的 最短路径问题，利用动态规划算法可以求解该问题，并且借助CUDA架构，每条二维线分 配一个线程负责计算对应的三维线，这样大大加快了计算效率。

3）三维特征线的聚类与拟合

步骤（5）、利用卷积曲面的方法，定义一标量场，该标量场每个采样点的标量值是衡量 其空间范围内的曲线密度，即曲线的多和少，并根据该标量场，保留标量值大的采样点，对 这些采样点进行细化操作，进一步化简。

步骤（6）、对于步骤（5）细化后的采样点，类似于模型的中心线，但是这些点并没有 连接信息，我们采用最小生成树算法来对这些点进行连接，首先为这些点构建一个图G，图 的每个节点代表每个采样点，而图的边表示采样点的连接关系，边的权重表示两个点的距离， 为该图可以利用Prim算法构建最小生成树T，并根据T找到叶节点和分叉点，从而能生成 曲线。

步骤（7）、对于步骤（6）得到的曲线，由于是把各个点进行连接得到，这样的曲线是 不够圆滑的，我们采用B样条曲线对其进行拟合，原来曲线上的点相当于控制点，这样能够 生成更加光滑的曲线，使得绘制效果更加逼真。

总之，与之前的方法相比，本发明能够快速地绘制体数据非真实感效果并且能得到较好 的效果。本发明主要有两点贡献：第一，给出了一种出图像空间计算三维特征线的方法，该 方法是从直接体绘制的渲染效果，并提取二维特征线出发，具有快速并较准确的优点。第二， 利用卷积曲面方法对三维曲线进行聚类并利用B样条进行拟合，使得绘制效果更加良好和逼 真。

附图说明

图1为算法整体流程图；

图2为光线投射直接体绘制原理图；

图3为渲染二维特征线效果图；

图4为计算三维特征线示意图；

图5为三维特征线效果图；

图6为采样点的曲线密度场效果图；

图7为对场进行细化效果图；

图8为对细化后的点连接成曲线效果图；

图9为Bezier曲线拟合效果图。

具体实施方式

下面结合附图以及本发明的具体实施方式进一步说明本发明。

对于由体素组成的给定体数据（图2的voxel grid），本方法针对该体数据的非真实感绘 制包括以下步骤:

步骤（1）我们设置六个视点（图2中的eye点，其他五个视点未示出）在X，Y，Z坐 标轴正负的六个方向，对场景进行光线投射（Ray Casting）绘制。该方法是基于图像空间的 算法，从图像的每一个像素，沿着视线方向发射一条光线，光线穿越整个场景，并在这个过 程中，对场景进行采样获取颜色信息，同时依据光线吸收模型将颜色值和不透明度进行累加， 直到光线穿出场景。其中颜色和不透明度的累加公式如公式（1）（从背面到前面进行采样合 成）：

$C_i^{\Delta }=(1-A_{i-1}^{\Delta })C_i+C_{i-1}^{\Delta }$

$A_i^{\Delta }=(1-A_{i-1}^{\Delta })A_i+A_{i-1}^{\Delta }---\left(1\right)$

其中C，A分部表示颜色和不透明度，光线在体数据穿射过程中，是按照一定步长进行采样 的，存在公式（2）：

t＝t_start+d*delta   (2)

其中t_start是光线与场景的交点初始位置，d是光线方向矢量，delta则为步长。该方法完全 基于GPU加速，大大加快了渲染速度，得到多个二维渲染结果，如图3（这里只示意了一个）。

步骤（2）、对步骤（1）的二维效果图利用图像处理相关技术，主要是Canny算子和Contour  Detect算法对二维图进行特征线的提取。Canny算子也是一阶算子，实质是用一个准高斯函 数做平滑运算，然后用带方向的一阶微分算子来定位导数最大值。具体求解步骤如下：

A)用高斯滤波器平滑图像，可用公式（3）：

$h(x,y,\sigma )=\frac{1}{{2\mathrm{\pi \sigma }}^2}{\theta }^{-\frac{x+y}{{2\sigma }^2}}$

g(x,y)=h(x,y,σ)*f(x,y)   (3)

其中h(x,y,σ)表示对原图像f(x,y)的平滑，g(x,y)为平滑后的图像，“*”代表卷积。

B)用一阶偏导有限差分计算梯度幅值和方向

C)对梯度幅值进行非极大值抑制

D)用双阈值算法检测和连接边缘

在对二维效果图利用Canny边缘检测后，再利用轮廓线检测（Contour Detect)算法对这 些特征点进一步计算，可用提取出特征线，效果图如图4，图5，图5中同一颜色的线代表 为一特征线。

步骤（3）二维特征线在二维上是具有明显视觉意义的线，我们认为这些线是三维的线 在二维的投影效果，所以需要找到每条二维特征线对应三维的特征线。如图4（a），其中在 X-视点下，计算得到一条二维的特征线（黄色线所示），我们需要在三维中计算出该线对应 的曲线（红色虚线所示）。然而我们注意到，一个二维点是很多个三维点的投影，也就是说 一个点对应着三维的很多点，如图4（a）所示的绿色点。

该问题可以转化为一个全局最小化问题，即使得所有阶段的全局代价最小（路径最短）， 该问题可以转化为最短路径问题，如图4（b）所示，红色的路径即我们需要的路径。如果单 纯使用暴力求解，则需要O(n!)复杂度。我们动态规划算法可以求解该问题，有如下方程：

Dist^[k][x]=min{Dist[y]+Weight(x,y)},y∈k+1阶段节点   (4)

其中Dist[x]表示节点x到终节点的最短路径长度，Weight(x,y)表示节点x，y之间的路径， 由于一个三维点有位置信息，不透明度信息，颜色信息，这些是在最小化问题需要考虑的因 子。所以可以有：

Weight(x,y)=w1*|x.col-y.col|+w2*|x.pos-y.pos|   (5)

其中col和pos分别表示颜色RGB和不透明度及空间位置，w1和w2代表权重。我们借助 基于GPU的CUDA架构，每条二维线分配一个线程负责计算对应的三维线，这样大大加快 了计算效率，效果如图5。

步骤（4）、可以看到，初步得到的三维特征线，能够较好的反应出体数据的结构，但是 这些线还比较多和杂乱，因此需要对这些线进一步聚类。我们采用卷积曲面的方法，定义一 标量场，该标量场衡量其空间一定空间范围内的曲线密度，即曲线的多和少。卷积曲面定义 如下：设P∈R³,f:R³→R为一个势函数，g:R³→R为表示某一个几何骨架的三变量 函数，

$g\left(P\right)=\left(\begin{array}{cc}1,& P\in \Omega \\ 0,& \mathrm{otherwise}\end{array}\right)---\left(6\right)$

Ω为骨架，点Q∈Ω，由骨架Ω产生的势函数定义为：

F(P)＝∫_Ωg(Q)f(P-Q)dΩ   (7)

写成卷积形式为

$F\left(\mathrm{PP}\right)=(f\oplus g)\left(P\right)---\left(8\right)$

我们采样的势函数为高斯核函数，如公式（9）所示：

$f\left(x\right)=\frac{1}{\sigma \sqrt{2\pi }}{\theta }^{-\frac{{(x-\mu )}^2}{{2\sigma }^2}}---\left(9\right)$

构建该标量场的时候，为每个采样点P建立半径为R的采样球，当有曲线穿越过该球时，会 对该曲线进行积分操作，如公式（10），|p-q|表示两点的欧式距离。

$f\left(P\right)={\int }_{\Omega }\frac{1}{\sigma \sqrt{2\pi }}{\theta }^{-\frac{{|F-Q|}^2}{{2\sigma }^2}}\mathrm{d\Omega }---\left(10\right)$

并根据该标量场，保留标量值大的采样点，留作进一步化简。见效果图6，其中颜色越白的 点表示该点曲线密度值越高，也就是越多曲线在该点附近区域。

步骤（5）、对于步骤（4）得到场内采样点，仍然较多，我们采用细化算法对这些点进 一步化简，从中提取出能够反应体数据结构的点。细化算法类似于“剥洋葱”的过程，经过一 层层的剥离，从原来的图像中去掉一些点，但仍要保持原来的形状，直到得到图像的骨架， 它可用于二维图像和三维体数据，对本发明而言，一个采样点是否应该被删除,取决于以下条 件:

a)该点是否是表面的点，该判断是根据该点的空间26-邻居属性。

b)该点不能是一条曲线的起点或者终点。

c)删除该点不会改变欧拉示性数（Euler characteristic)，例如，删除该点不会导致原数 据产生空洞。

d)删除该点不会改变联通联结数，联结数是指在该点的3*3*3领域中，和该点连接的 图形分量的个数。

如果满足以上条件，则我们可以删除该采样点，直到所有的点都不能删除为止，那么剩下的 点就是我们需要的中心线上的点。见效果图7。

步骤（6）对应步骤（5）求得的细化后的点，能够得到该场的中心线，但是这些点并没 有连接信息，我们采用最小生成树算法来对这些点进行连接，步骤如下：

a)首先为这些点构建一个图G(V,E)，V代表采样点，E代表采样点的连接关系。边的权 重表示两个点的欧式距离，如果两个点距离过大，则认为是不连接的。

b)利用Prim算法为图G构建最小生成树T，在本方法中，该图是无向非连通图，所以 会有多颗树。

c)根据这些最小生成树，找到叶节点和分叉点，并遍历树，计算出每个曲线经过的节点， 从而生成曲线。

假设其中某颗最小生成树如图8（a）所示，则其中红色的点则为分叉点和叶节点，则可 以分段构建曲线路径：

曲线0：包含节点1-2-3；曲线1：包含节点3-4-5；曲线2：包含节点3-6-7； 曲线3：包含节点1-8；曲线4：包含节点8-9-10；曲线5：包含节点8-13-14； 曲线6：包含节点8-11-12。效果图见图8（b）。

步骤（7）对于步骤（6）得到的曲线，由于是把各个点进行连接得到，这样的曲线是不 够圆滑的，我们采用B样条曲线对其进行拟合，我们采用的是Bezier曲线。

假设给出n+1个控制点位置：P_k=(x_k，y_k，z_k)，这里k可以取0到n。这些点将混合产 生下列位置P(u)，用来描述P₀和P_n间逼近Bezier多项式函数的路径：

$P\left(u\right)={\Sigma }_{k=0}^n{P}_k{\mathrm{BEZ}}_{k,n}\left(u\right),0\le u\le 1---\left(11\right)$

Bezier混合函数BEZ_k，n(u)是Bernstein多项式：

BEZ_k，n(u)＝C(n，k)u^k(1-u)^n-k   (12)

这里，参数C(n,k)是二项式系数：

$C(n,k)=\frac{n|}{k\left|(n-k)\right|}---\left(13\right)$

在本发明，原来曲线上的点相当于Bezier曲线的控制点，然后利用公式（11）（12）（13） 插值计算出更多的点，这样能够生成更加光滑的曲线，使得绘制效果更加逼真，效果图见图 9。

其中，步骤（1）（3）（4）（7）计算均是用CUDA实现的，可以利用GPU强大的计算 能力，大大提高了效率。

本文算法的实现使用的软件平台为Microsoft visual studio2010与OpenGL，使用了 CUDA来加速并行算法的计算效率。硬件平台为3.4GHz Inter(R)Core(TM)i7-2600CPU、4GB 内存以及NVIDIA GeForce GTX570GPU。方法效果图如图9所示。输出屏幕分辨率为 512*512，体数据三维大小为256*256*128，其中利用光线投射算法进行直接体绘制，并利用 Canny算子和轮廓线检测算法计算二维特征线的时间大概是0.70s，利用动态规划求解二维线 对应的三维特征线时间大概是1.22s，利用卷积曲面方法构建密度场时间大概为5.26s，其中 对场进行细化，并最后拟合的时间大概为1.60s。