Using NURBS Surfaces in Real-Time Applications

What About Surface Normals?

So now that we have a general idea of a way to tessellate a NURBS surface (or any other parametric surface, for that matter), what else do we need? For one, we need a way to generate surface normals so that we can let the 3D API (Direct3D* in the sample code) do lighting calculations for us. How do we generate these? Well, remember those Calculus classes that we all loved? One of the things we learned is that the derivative of a function is the instantaneous slope of the line tangent to the function at the point where the derivative and function are evaluated. By creating two tangent lines (one in the u and one in the v parameter) we can take a cross product and wind up with a surface normal. Simple enough, you say, but what's the derivative of the function S(u,v)?

Well, there are two partial derivatives: one with respect to u andone with respect to v, and they're ugly! Using the chain-rule:

And, not only is that ugly, we don't really know how to take the derivatives of Bi,p(u) and Bj,q(v). It's possible to take a derivative of Bi,p(u) (and Bj,q(v) ) from it's definition, but there's an easier way. It's possible to come up with a set of equations for calculating the coefficients of the polynomial equation that Bi,p(u) is equivalent to. Then, taking the derivative of BI,p(u) is as simple as multiplying powers by coefficients and reducing the powers by one (if you recall d(Ax^n + Bx^m)/dx = nAx^(n-1) + mBx^(m-1)). You still have to use Equation 5 to compute the derivatives of S(u,v) but it's really not that bad - you're going to be performing the computation of some of the terms any way, and the ones with the derivatives are calculated the same way as the non-derivative terms. We need to be able to calculate the coefficients of the b-spline basis functions when they're represented as follows:

Using a lot of paper and a bit of head scratching, I derived the following formulas to compute the coefficients, Ci,p,k(u).

This seems complex, but unless the knot vector changes, you don't have to re-compute these coefficients after the first time. Also note that Ci,p,k is only dependent on u for the knot span that u is in not on u itself, so we can just evaluate the Ci,p,k for each knot span and store those values. Now we can write the derivative of Bi,p(u) as:

Sample Code

At this point we know what we need to know to talk about the sample code you can download and how to implement this fun stuff. First, everything in the sample code is written in C++ and spread across many files of which mainly two are specific to this article: DRGNURBSSurface.h and DRGNURBSSurface.cpp. Actually, you'll also dive into NURBSSample.cpp if you want to play with the surface control points and knot vectors. DRGNURBSSurface.h contains a class definition for a class called CDRGNURBSSurface (for the curious, C is for "class", DRG is for "Developer Relations Group" which is what the group I'm in at Intel used to be called). The methods of this class of interest to us are Init(), ComputeBasisCoefficients(), ComputeCoefficient(), SetTessellations(), EvaluateBasisFunctions(), TessellateSurface(), and TessellateSurfaceSSE().

Going through these in order, Init() is called to initialize a newly created CDRGNURBSSurface object. The function takes a pointer to a CDRGWrapper class that is part of the framework we wrote for getting at the Direct3D* API. Init() also takes two surface degrees, u and v, and the number of control points in the u and v directions. It takes an array of Point4D structures that contain the weighted control points (x, y, z, and w) stored in u-major order (this means that v values are consecutive in the array). It takes two float arrays that contain the u knots and the v knots. Finally, it takes two optional values that specify the number of tessellations in the u and v directions of the surface. Init() does some calculations to determine how many knots are in the knot vectors and then allocates memory to store some of the information needed to render the surface. Finally, Init() makes a local copy of the incoming data (control points and knots) and then calls ComputeBasisCoefficients().

ComputeBasisCoefficients() calls ComputeBasisCoefficient() which uses the formulas from Equation 6 to compute the coefficients of the polynomials formed from the knot vectors and the degrees of the surface. ComputeBasisCoefficient() calls itself recursively due to the definitions in Equation 6. The coefficients are stored in arrays to be used by EvaluateBasisFunctions(). Because the Ci,p,k(u) are only dependent on the knot span that u belongs in, ComputeBasisCoefficient() takes as an argument this knot span (referred to as an "interval" in the code) rather than the actual value of u.

After Init() has called ComputeBasisCoefficients() to do the one-time calculation of the polynomial coefficients, SetTessellations() is called to set the number of u and v tessellations that will be used for rendering the surface. SetTessellations() can be called at any time after initialization to change the fineness of tessellation of the surface. The sample application calls SetTessellations() whenever the plus key (+) or minus key (-) is pressed to increase or decrease the tessellation of the surface. SetTessellations() allocates memory that's dependent on the number of tessellations used for rendering the surface, sets up some triangle indices for rendering the surface, and then calls EvaluateBasisFunctions().

EvaluateBasisFunctions() uses the coefficients computed in ComputeBasisCoefficients() and a technique called "Horner's method" to evaluate the polynomials that are the expanded form of the basis functions. Horner's method says that f = anx^n+an-1x^(n-1)+…+a1^x + a0 can be evaluated using n multiplications and n additions by rewriting as f = a0+x*(a1+x*(a2+…x*(an-1+x*an)…)). If you think you'll be calling EvaluateBasisFunctions() often because your tessellations will be changing, then other optimizations could be made here (e.g. using a technique called "forward differences" to eliminate the multiplications in the inner loop). Additionally, this method could be optimized using the Streaming SIMD Extensions of the Intel Pentium III processor.

At this point, everything is initialized for tessellating a NURBS surface. Now, at each frame that the sample application renders, the Render() method of the CDRGNURBSSurface object is called and in turns calls TessellateSurface() or TessellateSurfaceSSE() depending on whether or not we've told the object to use the Streaming SIMD Extensions of an Intel Pentium III processor.

TessellateSurface() (or TessellateSurfaceSSE()) uses Equation 4 and Equation 5 to compute the surface points and derivatives at the tessellation steps. A cross-product of the derivatives is used to compute the normal to the surface. We don't check for degenerate normals (see the pitfalls section below) so you'll need to modify these routines if degenerate normals become an issue. During the tessellation, a row of triangle vertices is generated. We alternate between putting the vertices in the odd or the even indices of the vertices buffer. Starting with the second row of generated vertices, we call Direct3D* to render a triangle strip using the strip indices generated in SetTessellations(). We alternate between the sets of indices as well due to the winding order of the triangle strip.