CN101599181B - 一种代数b样条曲面的实时绘制方法 - Google Patents

Abstract

本发明公开了一种基于NR迭代算法的ABS实时绘制方法，包括七步：1)输入待绘制ABS的相关信息并计算Lipschitz常数，将ABS转换为分片的Bezier曲面片；2)使用MC算法抽取ABS的等值面；3)实时计算ABS的极曲面；4)对等值面进行绘制，得到等值面及其侧影轮廓线的初始值；5)聚集等值面的侧影轮廓线的初始值，根据ABS及其极曲面方程，使用NR迭代算法对初始侧影轮廓线进行计算得到浮点精度的侧影轮廓线，并对其进行校正；6)计算侧影轮廓线的每个像素视线和Bezier曲面的交点；7)通过NR迭代算法利用邻接的Bezier曲面信息重新进行迭代求精，得到所有光线与代数B样条曲面的正确交点，然后再利用交点信息与光源、视点的相对位置以及曲面本身的材质计算光照。本发明方法能实现任意次数的代数曲面的绘制。

Description

一种代数B样条曲面的实时绘制方法

技术领域

本发明涉及计算机图形学实时绘制技术，特别是涉及一种基于Newton-Raphson(简称NR)迭代算法的代数B样条曲面的实时绘制方法。

背景技术

代数曲面的光线投射算法的核心问题为光线跟曲面的求交计算。由于低于五次的一元方程可以解析的求根，Blinn近年来在解析求根的理论基础与鲁棒性上也作了大量的工作，基于此Loop提出了基于GPU的分段代数曲面的解析的鲁棒的实时绘制算法，参见Loop C.，Blinn J.：Real-timeGPU render-ing of piecewise algebraic surfaces.In SIGGRAPH’06，pp.664-670，但是这种方法无法做次数高于4次的代数曲面的绘制。

对于一般代数曲面的求交算法，Kajiya 1982年中使用的算法为Laguerre迭代，该迭代提供了三次的收敛速度，但是每次迭代的代价比较高，参见Kajia J.T.：Ray tracing parametric patches.In SIGGRAPH’82，pp.245-254。Knoll分别在2006年与2007年利用SSE指令集同时进行4条光线的基于区间算法的求交计算，对于超二次曲面达到30fps的性能，随后将其推广到GPU上，参见Knoll A.，Wald I.：Interactive ray tracing of arbitrary implicitfunctions.In Proceedings of the 2nd IEEE/EG Symposium on Interactive RayTracing(2007)，pp.11-18和Knoll A.，Hijazi Y.，Kensler A.，Schott M.，Hansen C.，Hagen H.：Fast ray tracing of arbitrary implicit surfaces withinterval and affine arithmetic.Comput.Graph.Forum 28，1(2009)，26-40，区间算法的缺陷在于需要对于特定的曲面设计特定的区间操作，尤其对于系数很多的曲面区间算法过于保守，导致效率低下。Nishita于1989年使用了Bezier Clipping算法在Bezier曲面的u，v方向上进行裁剪细分进行求交，Nishita T.，Sederberg T.W.，Kakimoto M.：Ray tracing trimmed rationalsurface patches.In SIGGRAPH’90，pp.337-345，该算法需要递归的调用Bezier clipping，不适合并行的GPU运行环境。Reimers等于2008年使用B样条节点插入算法得到实时绘制的效果，Reimers M.，Seland J.：Raycasting algebraic surfaces using the frustum form.Comput.Graph.Forum 27，22008)，361-370，但该方法仅限于单片曲面，而且会引入非常多的复杂的函数复合操作，在求法向与反走样方面表现也不好。

具我们所知，至今没有专门针对高次的分段连续的代数B样条曲面的实时绘制的算法。

发明内容

本发明提出了一种基于NR迭代算法的代数B样条曲面的实时绘制方法，该算法是基于GPU的光线投射算法，适合于绘制分段连续的高次代数曲面。

绘制曲面的方法有两种，一种是绘制粗糙的近似的模型，因为近似的逼近得到的模型相对简单，绘制也方便，但是绘制精度不够；另外一种就是光线投射(或光线跟踪)绘制算法，通过求解视线与曲面的交点达到在屏幕上像素精度的绘制结果，这种方法较为复杂耗时，但结果非常漂亮。因为代数B样条曲面(即代数张量积B样条曲面)具有分段连续性以及一般次数会超过4次(高于4次的方程没有解析解)的性质，所以本发明方法适合分段高次代数曲面的绘制。

由于在曲面空间使用NR迭代算法，在分析了NR迭代算法收敛性要求的基础上提出了分区域绘制的算法，该算法可以保证NR迭代算法的收敛性。

为了达到上述目的，本发明采用以下技术方案如下：

包括以下七个步骤：

(1)输入待绘制的代数B样条曲面的相关信息，计算得到代数B样条曲面的Lipschitz常数，并通过节点插入算法将代数B样条曲面转换为分段的Bezier曲面片；

所述的代数B样条代数曲面为如下定义：

曲面定义于三维空间的一个矩形域R₃＝[a₁，a₂]×[b₁，b₂]×[c₁，c₂]上，令X＝[x₀，x₁，…，x_m+M+1]，Y＝[y₀，y₁，…，y_n+N+1]，Z＝[z₀，z₁，…，z_q+Q+1]表示在x、y和z方向上非递减节点向量，则一个分段的代数张量积B样条曲面表示如(1)式所示：

$F(x,y,z)=\underset{i=0}{\overset{M}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}\underset{k=0}{\overset{Q}{\mathrm{\Σ}}}{w}_{\mathrm{ijk}}{N}_i^m\left(x\right)N_j^n\left(y\right)N_k^q\left(z\right)=0......\left(1\right)$

其中分别表示在节点向量X、Y和Z上的m、n和q次的B样条基函数。标量w_ijk是曲面的权值，与参数曲面的控制顶点有相似的作用。该曲面是张量积的，各个方向都是独立的，曲面的次数为X、Y和Z三个轴方向上的次数和为M+N+Q，M+N+Q通常会大于4，所以代数B样条曲面具有次数高、各个方向上高阶连续性等特点。而代数B样条曲面可以通过节点插入算法得到等价的分段连续的Bernstein基表示的Bezier曲面。

NR迭代算法可以提供二次的收敛速度，也是最常见的求解非线性方程组的方法，但是其算法收敛的要求为良好的初始值，否则算法可能发散。通过分析NR迭代算法的局部收敛特性性质，得到对一般非线性方程组，只需满足如下(2)式，算法即可收敛：

||J(x_*)^-1[J(x)-J(x_*)]||≤||J(x_*)^-1||||[J(x)-J(x_*)]||  ......(2)

≤βγ||x-x_*||

其中该系统的解为x_*，并且有γ，β＞0，使得Jacobi矩阵J(x_*)^-1存在且||J(x_*)^-1||≤β，J∈Lip_γ(N(x_*，r))。

那么我们在屏幕每一个像素上绘制曲面的时候，需要求解通过屏幕像素的视线方程，与曲面联立的非线性方程组，经过推导可得该方程的Lipschitz常数(即βγ)。

所述的屏幕像素的视线方程如(3)式所示：

x_*＝E+dt  ......(3)

其中E为视点坐标，d为视线方向；

所述的Lipschitz常数(即βγ)如(4)式所示：

$\mathrm{\β\γ}=\left|\right|J{\left(x_{*}\right)}^{-1}\left|\right|\frac{\left|\right|[J\left(x\right)-J\left(x_{*}\right)]\left|\right|}{\left|\right|x-x_{*}\left|\right|}=\frac{d^T\·\left|\right|H_F\left|\right|\·d}{\▿F\·d}\≤\frac{1}{\mathrm{\ϵ}}\max \left(\right|\left|H_F\right|\left|\right)......\left(4\right)$

其中

为代数B样条代数曲面函数F(x，y，z)的梯度，H_F为代数B样条代数曲面函数F(x，y，z)的Hessian矩阵，d为绘制的时候相对于每一个像素的视线方向，是d的转置矩阵，

为视线方向与曲面法向的夹角，ε的设定通常会影响线性插值区域的大小，ε值越大，线性插值来逼近的区域也就越大，当视线方向与曲面垂直的时候ε为0，也就是侧影轮廓线的位置，通常选取ε小于0.1就可以得到满意的结果。

在计算的Lipschitz时，涉及到代数B样条曲面的对于各个方向的二阶偏导，且需要计算在其定义域中的Hessian矩阵的最大范数，准确的求解该范数的复杂度非常高，那么我们利用代数B样条曲面的凸包性，可以使用代数B样条曲面每个方向二阶偏导的系数(即Hessian系数)代替函数值来近似的逼近Lipschitz常数，这样计算的结果会放大该常数，导致的结果为MC的分辨率变高，经大量实验验证，取这种方法得到的常数的一般即可满足一样的收敛情况。

(2)使用Marching cubes算法绘制抽取代数B样条曲面的等值面，所述等值面的分辨率为所述的代数B样条曲面的Lipschitz常数；

使用Marching cubes算法(简称MC算法)抽取代数B样条代数曲面(简称ABS)的近似逼近，然后使用该等值面作为NR迭代算法(Newton-Raphson迭代算法)的初始值，并计算等值面上所有点的法向。所述的法向

F_x、F_y、F_z各个分量为曲面表达式在x、y、z三个方向上的偏导。本方法使用的MC方法得到的ABS(代数B样条曲面)的等值面上点及其法向均为准确计算，比线性插值更精确。MC算法采样ABS(代数B样条曲面)的得到的等值面的分辨率使用的为满足NR迭代算法收敛要求的Lipschitz常数。

(3)实时的计算代数B样条曲面的极曲面；

代数B样条曲面的极曲面的表达式如(5)式所示：

$P(x,y,z)=E\·\▿F+F_w=0......\left(5\right)$

其中，F_w为代数B样条曲面所对应的齐次曲面对于齐次坐标w的一阶偏导，式中E代表视点的坐标，

代表法向，各个分量为曲面表达式在x、y、z三个方向上的偏导。

(4)使用DirectX流水线对步骤(2)中的MC算法得到的等值面进行绘制，利用其高效的硬件光栅化算法和Z-Culling技术，得到代数B样条曲面的等值面及其侧影轮廓线的初始值，所述的初始值能进行良好的NR迭代算法。

所述的侧影轮廓线是指在曲面上其法向与视线方向垂直的位置。

(5)聚集步骤(4)中得到的等值面的侧影轮廓线的初始值，根据代数B样条曲面及其极曲面方程，使用NR迭代算法对初始的侧影轮廓线进行计算，得到代数B样条曲面的等值面的浮点精度的侧影轮廓线，并对侧影轮廓线附近的值进行矫正。

浮点精度指的是上述方程的解得到的侧影轮廓线位置是浮点数，该精度远远大于屏幕上绘制的像素精度，适合进行反走样操作。

由于当Lipschitz常数(即βγ)在ε小于某一阈值的时候，Lipschitz常数就会变得很大，这样就要求MC所使用的分辨率很高，这样系统的开销就很大，否则后续的NR迭代算法就不能收敛，所以我们使用曲面的侧影轮廓线来保证绘制的正确性。

曲面的初始的侧影轮廓线为曲面与其极曲面的交线，如(6)式所示：

$\left\{\begin{array}{c}F(x,y,z)=0\\ P(x,y,z)=E\·\▿F+F_w=0\end{array}\right.......\left(6\right)$

其中E为视点的坐标，

为法向

各个分量为曲面表达式在x、y、z三个方向上的偏导，F_w定义为曲面所对应的齐次曲面对于齐次坐标w的一阶偏导。

由于在可见的两个侧影轮廓点之间，曲面是连续可见的，而NR迭代算法在侧影轮廓线附近是很难保证收敛的，所以我们在求取侧影轮廓线的基础上，使用线性插值来替换无法收敛的区域，对侧影轮廓线附近的值进行矫正。

对于侧影轮廓线的NR迭代算法，我们不仅使用了横向、纵向两个方向的轮廓线，而且要针对不同点的信息在不同的方向上使用NR迭代算法。因为侧影轮廓线的法向是垂直于视线方向的，其在屏幕平面上的投影就是其真实的法向，等值面的绘制结果的边界上点很接近侧影轮廓线，其法向也是的侧影轮廓线上法向的近似比较，那么我们就根据在MC算法中求得的法向，将其投影到屏幕空间，比较其在屏幕空间中x，y方向的分量大小来确定其是作为横向还是纵向侧影轮廓线。

(6)将校正后的初始值按照其所属的Bezier曲面片分类，进行聚类后使用NR算法迭代，得到每个像素的视线跟每个Bezier曲面的交点。

但是在某些Bezier定义域之内、与其他Bezier曲面相邻的边界上可能有些交点经过NR算法迭代到了邻接的曲面上。

NR算法的每一步迭代都计算了曲面的法向信息，并保存该信息。

由于ABS(代数B样条曲面)相对于一般的代数曲面具有较多的控制系数，而NR算法的每一次迭代步骤需要将ABS(代数B样条曲面)系数数据以及其极曲面数据载入到流处理器中计算函数值与法向信息，而使用现有的图形卡硬件由于其内存机制的限制，无法对每个线程分别载入数据，所以本方法采用了分片曲面绘制时的数据访问。

针对CUDA现有的计算能力，将并行的线程数设置为64的倍数，这样每一个线程内，访问速度最快的变量——寄存器变量或者shared内存变量的数目就有限制。如果不能使用上述两种变量，可以使用常数内存或者使用取纹理的方法来载入变量，这两种方法都有缓存机制加速访问，最慢的访问是通过直接读取显存的方式。

因为绘制和求解侧影轮廓线时需要使用代数B样条曲面及其极曲面函数的系数，对于3*3*3的曲面而言，需要4*4*4*2个浮点数，相比于常见的曲面，这些数据量还是比较大的。由于使用NR迭代算法，每一步程序主要代价就是函数值以及法向数据的计算，所以显存的优化访问是获取高性能的关键步骤。

显然在迭代开始之前将所有并行线程的相邻曲面以及极曲面(在计算侧影轮廓线时)数据都载入到处理单元中是不可取的，因为无法确定最近的曲面信息，所以每个线程最多可能有3*3*3-1个相邻曲面片，这数据规模是非常大的。如果在程序运行时在根据程序计算结果载入相应的数据，最好也只能使用取纹理的方法来访问，但是由于可能每个线程访问的变量位置不同，会影响程序的并行性。

本发明使用的算法是基于分片的算法，首先根据初始值将绘制数据按照其所属曲面片分类，并根据曲面片信息进行排序聚类，这样属于同一曲面的数据就可以并行处理，需要访问的曲面数据也都是一样的，可以通过最快的shared内存的方式载入。曲面信息的排序聚类我们使用的是cudpp提供的基数排序，这种排序算法使用的是基数排序算法，比较使用于整数的排序。

(7)通过NR迭代算法利用邻接的Bezier曲面信息重新进行迭代求精，得到所有视线与代数B样条曲面的正确交点，然后再利用交点信息与光源、视点的相对位置以及曲面本身的材质计算光照，得到代数B样条曲线。

在某些Bezier曲面片的定义域之内，经过NR算法迭代步骤(6)中得到侧影轮廓线的每个像素视线和每个Bezier曲面的交点到了邻接的Bezier曲面片，因此，需要使用NR迭代算法利用邻接的Bezier曲面片信息重新进行迭代求精，得到最优结果，该次NR迭代算法的内存可以使用取纹理的方法。

当Bezier曲面片在边界位置出现收敛到[0，1]定义域之外时，则在相邻曲面再进行一次使用NR算法利用邻接的Bezier曲面信息重新进行迭代求精，得到最优结果。

由于在NR迭代算法的过程中已经计算了每次迭代结果处的法向，所以只需要简单的将NR迭代算法的最后一步迭代的法向直接进行光照计算，这也免去了以往算法的一些参数变化或使用一元函数求交的相互转化的步骤，效率更高。

走样问题一般发生在曲面与背景，曲面自遮挡的地方，显然这些位置在我们的程序中，侧影轮廓线的位置都是已知的，而且已经具有浮点数的精度，所以可以在亚像素级别上做反走样，侧影轮廓线处的曲面的法向与其在屏幕空间中的法向一致，我们可以根据其关于屏幕像素坐标的梯度来做反走样。

本技术方案中的代数B样条曲面具有可表示任意拓扑的光滑曲面、便于进行外形分析与实体造型等优点，而实时绘制是此类曲面进入实用领域的关键。具有像素精度的光线投射方法可以实现代数曲面的高质量显示。本发明提出了代数B样条曲面的实时GPU光线投射绘制算法。曲面绘制的核心问题为光线与曲面的求交算法，由于高于4次的多项式无法解析求根，所以本发明使用迭代的NR迭代算法，而NR迭代算法的初值由Marching cubes(MC)算法对原始曲面的逼近得到。基于对NR迭代算法的收敛性分析，本发明算法可以得到给定曲面的能保证NR迭代算法收敛的Lipschitz常数，使用在该常数粒度下的MC算法得到的初始值可以保证使用NR迭代算法求根可以达到二次收敛的效果。对于不能满足的收敛性要求的初始值(在侧影轮廓线附近)，可以通过简单的修改来达到保证收敛的效果。与之前的算法相比，该算法能提供在原始的物体空间中二次的收敛速度，收敛更快，而且更适合在现在的图形卡上运行。本发明方法还可以一次处理大量的分段连续高次代数曲面，这是以前的算法没有提到的，并给出了分段多片曲面的一些实验性建议，可以简单的推广到任意次数的代数曲面的绘制。

附图说明

图1是本发明代数B样条曲面实时绘制方法的流程图；

图2是本发明侧影轮廓线附近插值算法示意图。

具体实施方式

下面结合附图对本发明实施例进行详细说明，本发明可以在家用电脑的图形卡硬件中并行处理。

如图1所示，一种基于Newton-Raphson迭代算法的代数B样条曲面的实时绘制方法，包括七个步骤：(1)输入待绘制的代数B样条曲面的相关信息，计算得到代数B样条曲面的Lipschitz常数，并通过节点插入算法将代数B样条曲面转换为分片的Bezier曲面片；(2)使用Marching cubes算法抽取代数B样条曲面的等值面，所述等值面的分辨率为所述的代数B样条曲面的Lipschitz常数；(3)实时的计算代数B样条曲面的极曲面；(4)使用DirectX流水线对步骤(2)中的MC算法得到的等值面进行绘制，利用其高效的硬件光栅化算法和Z-Culling技术，得到代数B样条曲面的等值面的侧影轮廓线的初始值；(5)聚集步骤(4)中得到的等值面的侧影轮廓线的初始值，根据代数B样条曲面的极曲面和代数B样条曲面得到初始的侧影轮廓线，使用NR迭代算法对初始的侧影轮廓线进行计算，得到代数B样条曲面的等值面的浮点精度的侧影轮廓线，并对侧影轮廓线附近的值进行矫正；(6)校正后的侧影轮廓线的初始值按照其所属的Bezier曲面片聚类使用NR迭代算法，得到侧影轮廓线的每个像素视线跟每个Bezier曲面的交点；(7)通过NR迭代算法利用邻接的Bezier曲面信息重新进行迭代求精，得到所有光线与代数B样条曲面的正确交点，然后再利用交点信息与光源、视点的相对位置以及曲面本身的材质计算光照。

具体步骤如下：

(1)输入待绘制的代数B样条曲面的相关信息，计算得到代数B样条曲面的Lipschitz常数，并通过节点插入算法将代数B样条曲面转换为分段的Bezier曲面片；

所述的代数B样条代数曲面为如下定义：

曲面定义于三维空间的一个矩形域R₃＝[a₁，a₂]×[b₁，b₂]×[c₁，c₂]上，令X＝[x₀，x₁，…，x_m+M+1]，Y＝[y₀，y₁，…，y_n+N+1]，Z＝[z₀，z₁，…，z_q+Q+1]表示在x、y和z方向上非递减节点向量，则一个分段的代数张量积B样条曲面表示如(1)式所示：

$F(x,y,z)=\underset{i=0}{\overset{M}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}\underset{k=0}{\overset{Q}{\mathrm{\Σ}}}{w}_{\mathrm{ijk}}{N}_i^m\left(x\right)N_j^n\left(y\right)N_k^q\left(z\right)=0......\left(1\right)$

其中

分别表示在节点向量X、Y和Z上的m、n和q次的B样条基函数。标量w_ijk是曲面的权值，与参数曲面的控制顶点有相似的作用。该曲面是张量积的，各个方向都是独立的，曲面的次数为X、Y和Z三个轴方向上的次数和为M+N+Q，M+N+Q通常会大于4，所以代数B样条曲面具有次数高、各个方向上高阶连续性等特点。而代数B样条曲面可以通过节点插入算法得到等价的分段连续的Bernstein基表示的Bezier曲面。

NR迭代算法可以提供二次的收敛速度，也是最常见的求解非线性方程组的方法，但是其算法收敛的要求为良好的初始值，否则算法可能发散。通过分析NR迭代算法的局部收敛特性性质，得到对一般非线性方程组，只需满足如下(2)式，算法即可收敛：

||J(x_*)^-1[J(x)-J(x_*)]||≤||J(x_*)^-1||||[J(x)-J(x_*)]||  ......(2)

≤βγ||x-x_*||

其中该系统的解为x_*，并且有γ，β＞0，使得Jacobi矩阵J(x_*)^-1存在且||J(x_*)^-1||≤β，J∈Lip_γ(N(x_*，r))。

那么我们在屏幕每一个像素上绘制曲面的时候，需要求解通过屏幕像素的视线方程，与曲面联立的非线性方程组，经过推导可得该方程的Lipschitz常数(即βγ)。

所述的屏幕像素的视线方程如(3)式所示：

x_*＝E+dt ......(3)

其中E为视点坐标，d为视线方向；

所述的Lipschitz常数(即βγ)如(4)式所示：

$\mathrm{\β\γ}=\left|\right|J{\left(x_{*}\right)}^{-1}\left|\right|\frac{\left|\right|[J\left(x\right)-J\left(x_{*}\right)]\left|\right|}{\left|\right|x-x_{*}\left|\right|}=\frac{d^T\·\left|\right|H_F\left|\right|\·d}{\▿F\·d}\≤\frac{1}{\mathrm{\ϵ}}\max \left(\right|\left|H_F\right|\left|\right)......\left(4\right)$

其中

为代数B样条代数曲面函数F(x，y，z)的梯度，H_F为代数B样条代数曲面函数F(x，y，z)的Hessian矩阵，d为绘制的时候相对于每一个像素的视线方向，是d的转置矩阵，

为视线方向与曲面法向的夹角，ε的设定通常会影响线性插值区域的大小，ε值越大，线性插值来逼近的区域也就越大，当视线方向与曲面垂直的时候ε为0，也就是侧影轮廓线的位置，通常选取ε小于0.1就可以得到满意的结果。

在计算的Lipschitz时，涉及到代数B样条曲面的对于各个方向的二阶偏导，且需要计算在其定义域中的Hessian矩阵的最大范数，准确的求解该范数的复杂度非常高，那么我们利用代数B样条曲面的凸包性，可以使用代数B样条曲面每个方向二阶偏导的系数(即Hessian系数)代替函数值来近似的逼近Lipschitz常数，这样计算的结果会放大该常数，导致的结果为MC的分辨率变高，经大量实验验证，取这种方法得到的常数的一般即可满足一样的收敛情况。

(2)使用Marching cubes算法绘制抽取代数B样条曲面的等值面，所述等值面的分辨率为所述的代数B样条曲面的Lipschitz常数；

使用Marching cubes算法(简称MC算法)抽取代数B样条代数曲面(简称ABS)的近似逼近，然后使用该等值面作为NR迭代算法(Newton-Raphson迭代算法)的初始值，并计算等值面上所有点的法向。所述的法向F_x、F_y、F_z各个分量为曲面表达式在x、y、z三个方向上的偏导。本方法使用的MC方法得到的ABS(代数B样条曲面)的等值面上点及其法向均为准确计算，比线性插值更精确。MC算法采样ABS(代数B样条曲面)的得到的等值面的分辨率使用的为满足NR迭代算法收敛要求的Lipschitz常数。

(3)实时的计算代数B样条曲面的极曲面；

代数B样条曲面的极曲面的表达式如(5)式所示：

$P(x,y,z)=E\·\▿F+F_w=0......\left(5\right)$

其中，F_w为代数B样条曲面所对应的齐次曲面对于齐次坐标w的一阶偏导，式中E代表视点的坐标，

代表法向，各个分量为曲面表达式在x、y、z三个方向上的偏导。

(4)使用DirectX流水线对步骤(2)中的MC算法得到的等值面进行绘制，利用其高效的硬件光栅化算法和Z-Culling技术，得到代数B样条曲面的等值面及其侧影轮廓线的初始值，所述的初始值能进行良好的NR迭代算法。

所述的侧影轮廓线是指在曲面上其法向与视线方向垂直的位置。

(5)聚集步骤(4)中得到的等值面的侧影轮廓线的初始值，根据代数B样条曲面及其极曲面方程，使用NR迭代算法对初始的侧影轮廓线进行计算，得到代数B样条曲面的等值面的浮点精度的侧影轮廓线，并对侧影轮廓线附近的值进行矫正。

浮点精度指的是上述方程的解得到的侧影轮廓线位置是浮点数，该精度远远大于屏幕上绘制的像素精度，适合进行反走样操作。

由于当Lipschitz常数(即βγ)在ε小于某一阈值的时候，Lipschitz常数就会变得很大，这样就要求MC所使用的分辨率很高，这样系统的开销就很大，否则后续的NR迭代算法就不能收敛，所以我们使用曲面的侧影轮廓线来保证绘制的正确性。

曲面的初始的侧影轮廓线为曲面与其极曲面的交线，如(6)式所示：

$\left\{\begin{array}{c}F(x,y,z)=0\\ P(x,y,z)=E\·\▿F+F_w=0\end{array}\right.......\left(6\right)$

其中E为视点的坐标，

为法向

各个分量为曲面表达式在x、y、z三个方向上的偏导，F_w定义为曲面所对应的齐次曲面对于齐次坐标w的一阶偏导。

由于在可见的两个侧影轮廓点之间，曲面是连续可见的，而NR迭代算法在侧影轮廓线附近是很难保证收敛的，所以我们在求取侧影轮廓线的基础上，使用线性插值来替换无法收敛的区域，对侧影轮廓线附近的值进行矫正。

如图2所示，该图为算法在某一扫描平面的剖面图，图中曲线F(x，y，z)＝0为所绘ABS的某一剖面，折线F_MC为MC算法得到的等值面的剖面，P_MC为等值面上满足的点，也就是绘制等值面得到的侧影轮廓线的初始值，P_ε为使用P_MC做初始值、通过NR算法得到的、通过屏幕像素(x_ε，y_s)的视线与曲面的真实交点，P_s为使用P_MC做初始值、通过NR算法得到的、处于y_s扫描线上的侧影轮廓线上的点。那么在通过(x_s，y_s)与(x_ε，y_s)之间的像素的视线与曲面交点的初始值就是在“MC等值面”上、P_MC点以后的虚线段表示的，因为无法满足Lipschitz条件，可能会导致NR算法无法收敛，我们使用P_s与P_ε的线性插值(即他们之间的连线)得到对P_s与P_ε之间的曲面的良好的近似，因为不收敛的区域相对较少，这种近似的逼近可以得到较理想的绘制结果。

对于侧影轮廓线的NR迭代算法，我们不仅使用了横向、纵向两个方向的轮廓线，而且要针对不同点的信息在不同的方向上使用NR迭代算法。因为侧影轮廓线的法向是垂直于视线方向的，其在屏幕平面上的投影就是其真实的法向，收敛图边界上点很接近侧影轮廓线，其法向也是的侧影轮廓线上法向的近似比较，那么我们就根据在MC算法中求得的法向，将其投影到屏幕空间，比较其在屏幕空间中x，y方向的分量大小来确定其是作为横向还是纵向侧影轮廓线。

(6)将校正后的侧影轮廓线的初始值按照其所属的Bezier曲面片聚类使用NR算法迭代，得到每个像素视线跟每个Bezier曲面的交点。

NR算法的每一步迭代都计算了曲面的法向信息，并保存该信息。

由于ABS(代数B样条曲面)相对于一般的代数曲面具有较多的控制系数，而NR算法的每一次迭代步骤需要将ABS(代数B样条曲面)系数数据以及其极曲面数据载入到流处理器中计算函数值与法向信息，而使用现有的图形卡硬件由于其内存机制的限制，无法对每个线程分别载入数据，所以本方法采用了分片曲面绘制时的数据访问。

针对CUDA现有的计算能力，将并行的线程数设置为64的倍数，这样每一个线程内，访问速度最快的变量——寄存器变量或者shared内存变量的数目就有限制。如果不能使用上述两种变量，可以使用常数内存或者使用取纹理的方法来载入变量，这两种方法都有缓存机制加速访问，最慢的访问是通过直接读取显存的方式。

因为绘制和求解侧影轮廓线时需要使用代数B样条曲面及其极曲面函数的系数，对于3*3*3的曲面而言，需要4*4*4*2个浮点数，相比于常见的曲面，这些数据量还是比较大的。由于使用NR迭代算法，每一步程序主要代价就是函数值以及法向数据的计算，所以显存的优化访问是获取高性能的关键步骤。

显然在迭代开始之前将所有并行线程的相邻曲面以及极曲面(在计算侧影轮廓线时)数据都载入到处理单元中是不可取的，因为无法确定最近的曲面信息，所以每个线程最多可能有3*3*3-1个相邻曲面片，这数据规模是非常大的。如果在程序运行时在根据程序计算结果载入相应的数据，最好也只能使用取纹理的方法来访问，但是由于可能每个线程访问的变量位置不同，会影响程序的并行性。

本发明使用的算法是基于分片的算法，首先根据初始值将绘制数据按照其所属曲面片分类，并根据曲面片信息进行排序聚类，这样属于同一曲面的数据就可以并行处理，需要访问的曲面数据也都是一样的，可以通过最快的shared内存的方式载入。曲面信息的排序聚类我们使用的是cudpp提供的基数排序，这种排序算法使用的是基数排序算法，比较使用于整数的排序。

(7)通过NR迭代算法利用邻接的Bezier曲面信息重新进行迭代求精，得到所有视线与代数B样条曲面的正确交点，然后再利用交点信息与光源、视点的相对位置以及曲面本身的材质计算光照，得到代数B样条曲线。

在某些Bezier曲面片的定义域之内，经过NR算法迭代步骤(6)中得到侧影轮廓线的每个像素视线和每个Bezier曲面的交点到了邻接的Bezier曲面片，因此，需要使用NR迭代算法利用邻接的Bezier曲面片信息重新进行迭代求精，得到最优结果，该次NR迭代算法的内存可以使用取纹理的方法。

当Bezier曲面片在边界位置出现收敛到[0，1]定义域之外时，则在相邻曲面再进行一次使用NR算法利用邻接的Bezier曲面信息重新进行迭代求精，得到最优结果。

由于在NR迭代算法的过程中已经计算了每次迭代结果处的法向，所以只需要简单的将NR迭代算法的最后一步迭代的法向直接进行光照计算，这也免去了以往算法的一些参数变化或使用一元函数求交的相互转化的步骤，效率更高。

走样问题一般发生在曲面与背景，曲面自遮挡的地方，显然这些位置在我们的程序中，侧影轮廓线的位置都是已知的，而且已经具有浮点数的精度，所以可以在亚像素级别上做反走样，侧影轮廓线处的曲面的法向与其在屏幕空间中的法向一致，我们可以根据其关于屏幕像素坐标的梯度来做反走样。

本技术方案中的代数B样条曲面具有可表示任意拓扑的光滑曲面、便于进行外形分析与实体造型等优点，而实时绘制是此类曲面进入实用领域的关键。具有像素精度的光线投射方法可以实现代数曲面的高质量显示。本发明提出了代数B样条曲面的实时GPU光线投射绘制算法。曲面绘制的核心问题为光线与曲面的求交算法，由于高于4次的多项式无法解析求根，所以本发明使用NR迭代算法，而NR迭代算法的初值由Marchingcubes(MC)算法对原始曲面的逼近得到。基于对NR迭代算法的收敛性分析，本发明算法可以得到给定曲面的能保证NR迭代算法收敛的Lipschitz常数，使用在该常数粒度下的MC算法得到的初始值可以保证使用NR迭代算法求根可以达到二次收敛的效果。对于不能满足的收敛性要求的初始值(在侧影轮廓线附近)，可以通过简单的修改来达到保证收敛的效果。与之前的算法相比，该算法能提供在原始的物体空间中二次的收敛速度，收敛更快，而且更适合在现在的图形卡上运行。我们的方法还可以一次处理大量的分段连续高次代数曲面，这是以前的算法没有提到的，并给出了分段多片曲面的一些实验性建议。该算法还可以简单的推广到任意次数的代数曲面的绘制。

Claims (4)

1.一种代数B样条曲面的实时绘制方法，该方法包括以下七个步骤：

(1)输入待绘制的代数B样条曲面的相关信息，计算得到代数B样条曲面的Lipschitz常数，并通过节点插入算法将代数B样条曲面转换为分片的Bezier曲面片；

(2)使用Marching cubes算法抽取代数B样条曲面的等值面，所述等值面的分辨率为所述的代数B样条曲面的Lipschitz常数；

(3)实时的计算代数B样条曲面的极曲面；

(4)使用DirectX流水线对步骤(2)中的Marching cubes算法得到的等值面进行绘制，利用其高效的硬件光栅化算法和Z-Culling技术，得到代数B样条曲面的等值面及其侧影轮廓线的初始值；所述的侧影轮廓线是指在曲面上其法向与视线方向垂直的位置；

(5)聚集步骤(4)中得到的等值面的侧影轮廓线的初始值，根据代数B样条曲面及其极曲面方程，使用Newton-Raphson迭代算法对初始的侧影轮廓线进行计算，得到代数B样条曲面的等值面的浮点精度的侧影轮廓线，并对侧影轮廓线附近的值进行矫正；所述的初始的侧影轮廓线为曲面与其极曲面的交线；

(6)校正后的初始值按照其所属的Bezier曲面片分类，进行聚类后使用Newton-Raphson算法迭代，得到每个像素的视线跟每个Bezier曲面的交点；

(7)通过Newton-Raphson迭代算法利用邻接的Bezier曲面信息重新进行迭代求精，得到所有光线与代数B样条曲面的正确交点，然后再利用交点信息与光源、视点的相对位置以及曲面本身的材质计算光照。

2.根据权利要求1所述的代数B样条曲面的实时绘制方法，其特征在于：所绘制的代数B样条曲面为如下定义：

曲面定义于三维空间的一个矩形域R₃＝[a₁，a₂]×[b₁，b₂]×[c₁，c₂]上，令X＝[x₀，x₁，…，x_m+M+1]，Y＝[y₀，y₁，…y_n+N+1]，Z＝[z₀，z₁，…，z_q+Q+1]表示在x、y和z方向上非递减节点向量，则一个分段的代数张量积B样条曲面表示如(1)式所示：

$F(x,y,z)=\underset{i=0}{\overset{M}{\mathrm{\Σ}}}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}\underset{k=0}{\overset{Q}{\mathrm{\Σ}}}{w}_{\mathrm{ijk}}{N}_i^m\left(x\right)N_j^n\left(y\right)N_k^q\left(z\right)=0......\left(1\right)$

其中分别表示在节点向量X、Y和Z上的m、n和q次的B样条基函数；标量w_ijk是曲面的权值；该曲面是张量积的，各个方向都是独立的，曲面的次数为X、Y和Z三个轴方向上的次数和为M+N+Q，M+N+Q通常会大于4，所以代数B样条曲面具有次数高、各个方向上高阶连续性的特点；代数B样条曲面可以通过节点插入算法得到等价的分段连续的Bernstein基表示的Bezier曲面。

3.根据权利要求1所述的代数B样条曲面的实时绘制方法，其特征在于：所述的代数B样条曲面的Lipschitz常数如(4)式所示：

$\mathrm{\β\γ}=\left|\right|J{\left(x_{*}\right)}^{-1}\left|\right|\frac{\left|\right|[J\left(x\right)-J\left(x_{*}\right)]\left|\right|}{\left|\right|x-x_{*}\left|\right|}=\frac{d^T\·\left|\right|H_F\left|\right|\·d}{\▿F\·d}\≤\frac{1}{\mathrm{\ϵ}}\max \left(\right|\left|H_F\right|\left|\right)......\left(4\right)$

其中x_*为系统的解，并且有γ，β＞0，使得J(x_*)^-1为Jacobi矩阵，且||J(x_*)^-1||≤β，J∈Lip_γ(N(x_*，r))；

其中

为代数B样条曲面函数F(x，y，z)的梯度，H_F为代数B样条曲面函数F(x，y，z)的Hessian矩阵，d为绘制的时候相对于每一个像素的光线方向，d^T是d的转置，为视线方向与曲面法向的夹角，ε的设定通常会影响线性插值区域的大小，ε值越大，线性插值来逼近的区域也就越大，当视线方向与曲面垂直的时候ε为0，也就是侧影轮廓线的位置。

4.根据权利要求1所述的代数B样条曲面的实时绘制方法，其特征在于：根据代数B样条曲面的极曲面和代数B样条曲面得到的初始的侧影轮廓线表达式如(6)式所示：

$\left\{\begin{array}{c}F(x,y,z)=0\\ P(x,y,z)=E\·\▿F+F_w=0\end{array}\right.......\left(6\right)$

其中E为视点的坐标，

为法向各个分量为曲面表达式在x、y、z三个方向上的偏导，F_w定义为曲面所对应的齐次曲面对于齐次坐标w的一阶偏导。

Images