**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.