# Background

## NURBS Surfaces

The purpose of this introduction is to establish notation and serve as a reference for concepts that are sensitive to off-by-one errors or are otherwise subject to definition conflicts/confusion, not to provide a guide for the newcomer to NURBS surfaces (mathematically inclined or otherwise). For more information, please consult the Wikipedia page or one of the many textbooks written on the subject ([Piegl&Tiller,1995] is a classic reference). Be aware that every textbook, web article, or paper on the subject will possibly use different notation. The key to avoiding off-by-one errors is to represent concepts in terms of definitions that everyone agrees on (degree=order-1) and then to search for the analogs of equations (1) and (3) in order to discover whether k or p corresponds to degree or order, whether control points are indexed from P_00 to P_nn or P_11 to P_nn, and so on. It is worth mentioning that Rhino goes beyond the call of duty to use abnormal conventions for the noble purpose of saving 4 floats per surface.

NURBS surfaces are formally defined as a continuous map between a two dimensional rectangle in UV space and a curvy, wavy surface (or at least one that is potentially curvy and wavy) in 3D space.

$S(u,v)=\frac{\sum_{i=0}^n\sum_{j=0}^mw_{ij}P_{ij}N_{ip}(u)N_{jq}(v)}{\sum_{i=0}^n\sum_{j=0}^mw_{ij}N_{ip}(u)N_{jq}(v)}$ (1)

or, if the weights are appended to control points to form 4-dimensional homogeneous coordinates and the projection ("weighing") step is made implicit,

$S(u,v)=\sum_{i=0}^n\sum_{j=0}^mP_{ij}N_{ip}(u)N_{jq}(v)$ (2)

where the basis functions Nip(u) are recursively defined:

$N_{i,0}(u)=\left\{\begin{array}{ll} 1& u_i\leq u(3)

$N_{ip}=\frac{u-u_i}{u_{i+p}-u_i}N_{i,p-1}+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}N_{i+1,p-1}$(4)

Note the three types of parameters are used to control this mapping:

• Control Points $P_{ij}, 0\leq i\leq n, 0\leq j\leq m$
• Weights wij
• Knots $U=\{u_0,u_1,\ldots,u_r\}, V=\{v_0,v_1,\ldots,v_q\}$

Every point on a NURBS surface is an affine combination of control points. Weights control the relative importance of a control point within its patch and knots determine the correspondence between regions in U,V space and control points. A concise and friendly introduction to knots is [Peterson,1990].

In contrast to Bézier surfaces where points on the surface are an affine combination of every point on the control net, points on a NURBS surface are a combination an order*order rectangular subset of the control net. This gives NURBS surfaces the desirable property of local editing: moving one control point does not move the whole surface, only the parts of the surface corresponding to overlapping order*order patches. It also gives NURBS surfaces the ability to force the surface to pass through a point (by stacking (order-1)*(order-1) control points on top of one another) or to embed a completely flat region (a coplanar rectangle of order*order points would guarantee at least a 1x1 planar region in the center). Since a positive affine combination of points lies within the convex hull of those points, a NURBS surface lies within the union of convex hulls of rectangular order*order patches of its control points. For instance, an order 2 NURBS surface is a quadrilateral mesh, since each point on the surface is a combination of the 2*2 neighboring control points.

The fact that every point of the surface is a weighted combination of an order*order patch of control points introduces a complication in the knot vectors where this patch is not well defined, i.e. where it refers to a hypothetical control point beyond the edge of the control net. Considering a NURBS curve with the single coordinate u for simplicity, we find that C(u) is only well defined in [up,urp): before this range of u, the weights don't sum to one, and after this range the definition will refer explicitly to control points that don't exist (assuming that there are #ctrl_pnts+order knots, or equivalently that (r+1)=(n+1)+(p+1)). There are several ways of dealing with this.

Uniform knot vectors $U=\{0,1,2,3,\ldots,r\}$

ignore the problem and simply restrict the domain with the consequence that the curve (surface) edges don't line up with control points at the edge of the control net. The default behavior of blender's "Add NURBS Surface" operator is to instantiate such a NURBS surfaces. I am in favor of changing this. Nevertheless, uniform knot vectors can be useful: if order-1 control points are repeated to form a loop ($P_{n-p}=P_0, P_{n-p+1}=P_1, \ldots, P_n=P_p$) then the curve itself will become cyclic (a possibly deformed cylinder or sphere).

Open/Endpoint/Pinned knot vectors $U=\{\underbrace{0,0,0}_{p+1},1,2,\ldots,\underbrace{r-2p,r-2p,r-2p}_{p+1}\}$

effectively collapse the periphery of the U,V domain (where the surface was not defined) to 0 area (0 length in the case of a curve), which has the effect of ensuring that the edge of the surface lines up with the peripheral points of the control net. In Blender, this corresponds to the "endpoint" option (if "endpoint" is enabled, the curve will be given an open knot vector).

Piecewise Bézier knot vectors $U=\{\underbrace{0,0,0}_{p+1},\underbrace{1,1}_{p},\underbrace{2,2}_{p},\ldots,\underbrace{k,k,k}_{p+1}\}, k=\frac{r+1}{2}-p$

have the effect of ensuring that the NURBS surface can be exactly represented by a piecewise Bézier "quilt." Note that since knots can be inserted without altering the shape of the curve/surface, every NURBS surface can be represented by a quilt of Bézier patches of the same order by repeatedly inserting knots until the knot vector takes the correct form. The Bézier quilt representation could be desirable for artists who want to control surfaces using position and slope at each gridpoint rather than the more abstract NURBS control net (which generally doesn't touch the surface at all). The Bézier quilt is also important in the theoretical and practical sense that it forms the basis of all recursive subdivision algorithms: the 3d location of subdivided 2d parameter points is required to form the subdivided mesh and the first derivative is required to rigorously estimate the error. The quilt is computed implicitly or explicitly through knot insertion, and while this might seem to be a cumbersome approach it turns out that knot insertion is local, fast, and is in fact computationally equivalent to the process of explicitly evaluating S(u,v) and its first derivative at the subdivision points ([Piegl&Tiller,1995])!

## Rendering

Tessellation refers, of course, to the strategy of rendering the surface by creating a polygon mesh to approximate it and then rendering the polygon mesh. This is the approach that is currently within the scope of the GSOC project since 1) the tessellation code would be required regardless of rendering strategy for use by the convert-to-mesh operator and miscellaneous modifiers, and 2) representing the surface by its tessellation ensures consistency between the viewport and the final rendering.

Exact Fragment Shaders for real-time rendering of trimmed NURBS surfaces have been investigated over the last 7-8 years for varying definitions of "exact." Guthe's 2005 thesis provides a very detailed investigation of the case where a tessellated surface is fed to the GPU and normals, trimming, and grazing are computed per-fragment (grazing is when a fragment corresponds to a ray that intersects the tessellated surface but not the exact surface -- these fragments can be discarded for smooth silhouettes). Guthe's work is contained in the OpenSG library and has been thoroughly tested. A possible future direction would be to integrate this into Blender's 3d view, but this will have to wait until code to recursively generate a trimmed Bezier Quilt is completed and bugfixed.

Exact raytracing of NURBS surfaces is not only possible but common practice in commercial CAD packages. This would be another future direction for the NURBS project. Care will need to be taken in order to ensure that tessellation, the exact fragment shader in the 3d view, and the exact raytracing code all produce consistent results without cracks. Implementing tessellation first allows us to temporarily bypass these complex issues.

## Tessellation Terminology

Adaptive tessellation algorithms tend to deposit more polygons in regions of high curvature and fewer polygons in regions of low curvature, often with some guaranteed error bound (as in [Piegl,1998]), while direct tessellation algorithms simply deposit a fixed number of polygons per control-point grid element (as in [Pigel,1995]). The latter strategy also has a guaranteed error bound, but it requires more polygons than an adaptive tessellation to achieve the same error bound.

Trimming is the practice of removing the part of the NURBS surface corresponding to the interior of a curve (Bézier, NURBS, polyline, etc) defined on the coordinate patch U,V. For instance, the above image of an adaptively tessellated surface is trimmed, while the above image of a directly tessellated surface is not. This is coincidental: adaptive/direct has nothing to do with trimmed/untrimmed. Often two types of trimming curves will be defined: one which runs along the perimeter of the surface and defines the region to include and another that defines interior regions to exclude. A different convention for the same notion displays either points with odd total winding numbers or points with even total winding numbers.

Coving (the gerund form of "to cove," not a misspelling of "covering") is the practice of adding triangles between a separately tessellated curve and interior region of a trimmed surface. It is important that there be a general mechanism for creating the coving triangles since tessellation of boundaries is often controlled independently of tessellation of the interior, either by virtue of having a separate error bound (as in [Pigel,1998]) or by virtue of a sewing constraint (as in [Kumar,1997]).

 Two examples of coving.

Sewing is the process of ensuring that the tessellations of two surfaces sharing an edge do not have gaps at the edge. This requires a degree of connectivity information between surfaces. [Kumar,1997] introduces the concept of Super-Surfaces which are defined as sets of NURBS surfaces along with the corresponding connectivity information, potentially inferred from boundary proximity. The concept of a "super-surface" maps very closely onto Blender's notion of a curve object (which can contain multiple actual curves/surfaces). The most intuitive UI would probably be one in which the user could define a "merge tolerance" for the super-surface (blender curve object) below which boundary curves would automatically be sewn together.

 Sewing two tesselated surfaces together

# Tessellation Algorithms

## Uniform Grid Cut

This is the approach used by the Blender Nurbana branch and it is the approach that was used by BRL-CAD until 2012. Since BRL-CAD underwent heavy use by the US defense department during that time, the uniform grid cut tessellation strategy has proven to be satisfactory for visualization.

It is the simplest of the three tessellation approaches, and in outline form the strategy looks like:

• Lay down a uniform grid of vertices in parameter space
• Add a generalized polygon "face" structure for every grid rectangle
• Tessellate the boundary curves into polylines according to sewing and/or error constraints. Keep track of vertices and edges.
• Add vertices at points of intersection between the uniform grid edges and the polyline trim curves.
• Update the face data structures and split them into sub-faces as necessary.
• Triangulate all faces and discard those that were flagged as lying on the wrong side of a trim curve

A formal review of this tessellation strategy can be found in [Pigel,1995].

 Uniform grid with trim polylines Triangulation of the result

## Delaunay Triangulation

Delaunay Triangulation is a very general tessellation strategy that can be applied to almost any surface. Essentially, points are generated along trimming curves ("edge points") and uniformly or adaptively throughout the UV domain ("fill points"). Points that lie within trimming curves are discarded through a scanline (if uniformly deposited) or individual inclusion/exclusion test. The remaining points are then fed to a restricted Delaunay Triangulation algorithm which computes face information, resolving between multiple possible face arrangements by maximizing the minimum angle of each triangulation. This condition leads to high quality meshes that tend to avoid long, skinny triangles. FreeCAD/opencascade use this strategy, although the code is complex and contains error bounding mechanisms that are tightly integrated with the opencascade Constructive-Solid-Geometry (CSG) routines. In other words, the code isn't very portable. There is probably a nice, portable, open-source code out there to do this, but I haven't found it and there isn't a pressing need as the Delaunay Triangulation approach occupies an unfavorable part of the simplicity vs quality diagram compared to Uniform Grid Cut (which wins on simplicity) and Recursive Subdivision (which wins on quality by way of having a good error bound).

 DT of High-Res Circle, Low-Res Square DT with only edge points DT with edge points and fill points

See [Sheng,1992] and [Laurent,2000] for a formal review.

## Recursive Subdivision / Piecewise Bézier

Recursive subdivision is the high-quality, high-complexity algorithm that has been the cornerstone of commercial NURBS codes such as Rhino and Moi. RSD and its descendants have also been the focus of the NURBS tessellation literature for the last decade. The approach goes like this:

• First, convert the surface to a Bézier Quilt
• Insert knots into the U,V knot vectors until they both take the Bézier form. This is computationally equivalent to evaluating positions and derivatives at the corners of each quilt patch.
• Perform a change of basis from the B-Spline basis to the Bernstein polynomial basis
• Compute the derivative of each Bézier patch (this is a fast operation, the derivative is itself a Bézier patch but of lower order)
• If bound provided by the maximum derivative is too large, subdivide the patch (similar to knot insertion)
• Recurse
• Tessellate the trim curve into a polyline, intersect the polyline with the KD-tree composed of Bézier patches
• Add vertices at points of intersection, update faces (may have to split to avoid T junction, see [YingLiang,2002])
• Triangulate all faces

 Uniform Subdivision Adaptive Subdivision Cracks (why explicit face structure must be kept)

Mesa's implementation is relatively straightforward and doesn't contain frills that make it difficult to uproot and transplant. Still, the process will not be trivial. In particular I note that the BRL-2012 implementer wrote 3,000 lines of debug code to visualize an incorrectly behaving subdivision process. I think it would be wise to implement the simpler Uniform Grid Cut algorithm first to easily obtain "known good" results that could be used for debug purposes.

A particularly good review is Guthe's 2005 thesis. Also see [Rockwood,1989], [Abi-Ezzi,1994], [Kumar,1995], [Kumar,1997], [Piegl,1998], [Kahlesz,2002], [Kanai,2007], and [Krishnamurthy,20009].

## Tessellation Algorithms: High-level Comparison

Algorithm Comparison [--Strategy 1: Uniform Grid Cut--] [--Strategy 2: Delaunay Triangulation--] [--Strategy 3: RSD / Piecewise Bézier--]
Complexity 1.0 3.0 6.0
Familiarity with Ref. Code? 90% 20% 80%
Estimated Portability of Ref. Code? Simple Complex Complex (tough to debug)
Literature References? + [Piegl,1995] + [Sheng,1992] [Laurent,2000] [Sinclair,2010] + [Rockwood,1989] [Abi-Ezzi,1994] [Kumar,1995] [Kumar,1997] [Piegl,1998] [Kahlesz,2002] [YingLiang,2002] [Balázs,2004] [Guthe,2005] [Kanai,2007] [Krishnamurthy,2009]
Uses Triangles Efficiently? - + +
Rigorous Error Bound? - Yes but very inefficient [Piegl,1995] - No + Yes
Compatible with Sewing? + Yes + Yes + Yes
Profitable GPU Implementation? - - + [Guthe,2005], [Guthe,2006]

# Notes on Reference Implementations

    Link: http://brlcad.org/d/about
Track Record: Excellent. Comes from US DOD, heavily used and tested for military CAD.

NURBS info site:
New NURBS implementation:
src/librt/primitives/brep
src/libbrep/opennurbs_ext.cpp
Old NURBS implementation (same algo as Nurbana):
Paper: Pigel, 1995
src/librt/primitives/bspline/nurb_tess.c
Trimmed NURBS raytracing:
src/librt/primitives/brep/brep.cpp:676      rt_brep_shot() -> utah_brep_intersect()


## Mesa-GLU

    Link: http://www.mesa3d.org
Track Record: Comes from SGI. Used in ???
glu NRUBS API: http://www.glprogramming.com/red/chapter12.html
Version history: http://mesa3d.org/VERSIONS

gluBeginSurface
gluNurbsSurface()                                 //glu-9.0.0/src/libnurbs/interface/glinterface.cc
NurbsTessellator::nurbssurface()     //glu-9.0.0/src/libnurbs/internals/nurbsinterfac.cc
// First thing it does is make the surface piecewise Bezier
Quilt::toBezier()                           //glu-9.0.0/src/libnurbs/internals/tobezier.cc
gluEndSurface
NurbsTessellator::endsurface()          //glu-9.0.0/src/libnurbs/internals/nurbsinterfac.cc
NurbsTessellator::do_endsurface()     //glu-9.0.0/src/libnurbs/internals/nurbstess.cc
Subdivider::beginTrims() //glu-9.0.0/src/libnurbs/internals/subdivider.cc
foreach trim:
Subdivider::beginLoop() //glu-9.0.0/src/libnurbs/internals/subdivider.cc
Subdivider::endLoop()  //glu-9.0.0/src/libnurbs/internals/subdivider.cc
Subdivider::endTrims()  //glu-9.0.0/src/libnurbs/internals/subdivider.cc
Subdivider::beginQuilts()  //glu-9.0.0/src/libnurbs/internals/subdivider.cc
Subdivider::endQuilts()  //glu-9.0.0/src/libnurbs/internals/subdivider.cc


I intend to base Blender's future recursive subdivision implementation off of Mesa's due to its simplicity and clarity. [Krishnamurthy,2009] notes a bug in Mesa's RSD code which initially led me to consider other implementations:

However, after corresponding with Krishnamurthy and playing with the code he generously provided, I have discovered the source of this bug/misfeature: GLU does not appropriately subdivide trimming polylines, potentially leading to dramatic under-tessellation near the trimmed edges (it subdivides NURBS curves just fine even if they happen to be polylines in disguise). Since this behavior is technically "correct," the bug/misfeature does not affect the suitability of Mesa's code as a reference.

 GLU NURBS Surface Trimmed by Polyline GLU NURBS Surface Trimmed by the same polyline recast in 2D NURBS form

## OpenSG

    Link: http://cg.cs.uni-bonn.de/en/projects/opensg-plus/
Track Record: Excellent. Written by Guthe (heavily published NURBS prof). Lots and lots of demos, papers. http://www.opensg.org/projects/opensg/wiki/Gallery
Readability: Decent, but spread out and contains easily mistaken similar code to what I want.

* Has fragment trimmer AND geometric trimmer
* Geometric code is in opensg/Source/System/NodeCores/Drawables/Nurbs/Internal/OSGNurbsPatchSurface.cpp
CNurbsPatchSurface::ConvertToBezier()
CNurbsPatchSurface::CalculateTrimmingLoops()
* Trimmer in opensg/Source/System/NodeCores/Drawables/Nurbs/Internal/OSGNurbsParSpaceTrimmer.cpp


    Link: http://freecadweb.org   http://www.opencascade.org
Track Record: Excellent. Written and actively maintained by industry consulting firm.
Readability: Eek! (Sprawling. Organized, but there's so much extra stuff to make it difficult to read.)
High Level Concept Intro: https://www.ljll.math.upmc.fr/~perronnet/mit/mit.html

The Delaunay triangulation code is used both for visualization and finite-element mesh generation
GeomAbs_BSplineSurface /  BRepMesh_FastDiscretFace::InternalVertices
opencascade-6.7.1/src/BRepMesh/BRepMesh_FastDiscretFace.cxx:903 (search for  thetype == GeomAbs_BSplineSurface)

Delaunay Triangulation code spends most of its time in: MeshAlgo_CircleInspector::Inspect


## Blender (Nurbana)

    Map Nurbana <-> Trunk
Where Nurbana's Object_NURBS class is mapped to blender data structures: Object_NURBS.h
Where the old blender "Nurb" structure resides: blenderNurbsStruct.h
Length(uv) <-> m_pntsuv[uv]

Where nurbana does tesselation:
nb{Grid,Trim}Tessellator::getVertexes
Grid calls: intern/nurbana/intern/NURBS_Generate.cpp:52         NURBS_Generate::Surface()
Trim calls: intern/nurbana/intern/nbTrimTessellator.cpp:55     nbTrimTessellator::update()

How nurbana does trimming:
Step 0: Convert blender NURBS info (order, num control points, cyclic, etc) into Nurbana-friendly form
Step 1: Tesselate the trim curves into polygons
nbPolygon p; // Accumulates the tesselated trimming curves
foreach (Object_NURBS* trim_curve)
trim_curve->getTess(p)
NURBS_Generate::Curve()
Step 2: Use the polygons to cut a grid
nbTrimPatch  patch(start, last, cyclic);          // nop
nbTessGrid   grid(lengthUV, start, last, step);   // Create a singly linked (->,^) hotpoint,hotrect grid
grid.merge(patch, p); // Add vertices at points of trim-polyline intersection, remove interiors
nbSubDiv.cpp:591 subDiv() // Called for each hotrect, checks for intersection
nbSubDiv.cpp:547 sliceAt() // Called for each intersection
// AT LEAST 1 LINGERING EDGE CASE
// NO SEWING STRATEGY
Step 3: Push the uv polygons forward into xyzw space, project to xyz space


# References

Many of these papers are locked behind the Elsevier iron curtain. Freely available papers are highlighted in green.

[Rockwood,1989] Rockwood, A., Heaton, K., & Davis, T. (1989). Real-time rendering of trimmed surfaces. In ACM SIGGRAPH Computer Graphics (Vol. 23, pp. 107–116). ACM. Retrieved from http://dl.acm.org/citation.cfm?id=74344

[Sheng,1992] Sheng, X., & Hirsch, B. E. (1992). Triangulation of trimmed surfaces in parametric space. Computer-Aided Design, 24(8), 437–444. Retrieved from http://www.sciencedirect.com/science/article/pii/001044859290011X

[Abi-Ezzi,1994] Abi-Ezzi, S. S., & Subramaniam, S. (1994). Fast Dynamic Tessellation of Trimmed NURBS Surfaced1. In Computer Graphics Forum (Vol. 13, pp. 107–126). Wiley Online Library. Retrieved from http://onlinelibrary.wiley.com/doi/10.1111/1467-8659.1330107/abstract

[Kumar,1995] Kumar, S., & Manocha, D. (1995). Efficient Rendering of Trimmed NURBS Surfaces. Computer-Aided Design, 27(7), 509–521. Retrieved from http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.56.2447

[Piegl&Tiller,1995] Piegl Les, Tiller Wayne (1995). The NURBS Book. Monographs in Visual Communication, Springer, 1995.

[Piegl,1995] Piegl, L., & Richard, A. (1995). Tessellating Trimmed NURBS Surfaces. Computer-Aided Design, 27(1), 16–26. Retrieved from http://www.sciencedirect.com/science/article/pii/0000000000000000

[Kumar,1997] Kumar, S., Manocha, D., Zhang, H., & Hoff III, K. E. (1997). Accelerated walkthrough of large spline models. In Proceedings of the 1997 symposium on Interactive 3D graphics (p. 91–ff). ACM. Retrieved from http://dl.acm.org/citation.cfm?id=253313

[Piegl,1998] Piegl, L., & Tiller, W. (1998). Geometry-based Triangulations of Trimmed NURBS Surfaces. Computer-Aided Design, 30(1), 11–18. Retrieved http://www.sciencedirect.com/science/article/pii/S001044859700047X

[Laurent,2000] Laurent, P.-J., & International Conference on Curves and Surfaces. <4, 1999, Saint-Malo> (Eds.). (2000). Curve and surface design: Saint-Malo 99; [Fourth International Conference on Curves and Surfaces]. Nashville, Tenn.: Vanderbilt Univ. Pr. Retrieved http://www.dtic.mil/cgi-bin/GetTRDoc?AD=ADA399461

[Kahlesz,2002] Kahlesz, F., Balázs, Á., & Klein, R. (2002). Multiresolution rendering by sewing trimmed NURBS surfaces. In Proceedings of the seventh ACM symposium on Solid modeling and applications (pp. 281–288). ACM. Retrieved from http://dl.acm.org/citation.cfm?id=566323