NURBS Boolean Evaluation Development Guide

The technical approach used for evaluating Boolean operations on NURBS entities involves calculating geometric intersections, trimming surfaces accordingly, and stitching together a resulting object. Boundary representation NURBS objects are composed of faces, edges, and vertices; and these topologically describe how surfaces, curves, and points are joined together to represent geometry. The geometry is used to find intersections. The topology is used in the application of Boolean logic. This is oversimplified as there are also trimming curves, loops, orientations, and other details involved, but this is nonetheless a useful foundation for understanding the basic approach employed.

BRL-CAD heavily integrates and extends the OpenNURBS Toolkit by Robert McNeel & Associates, developers of the proprietary Rhinoceros freeform NURBS modeling software. BRL-CAD heavily uses OpenNURBS for trimmed NURBS geometry representation (both in-memory and on-disk) and implements functionality not included with OpenNURBS such as ray tracing and Boolean evaluation. ^[1] Additional functionality implemented for BRL-CAD is primarily contained within the boundary representation and ray trace libraries (LIBBREP and LIBRT). ^[2]

The overall NURBS Boolean evaluation algorithm is based on the paper BOOLE: A System to Compute Boolean Combinations of Sculptured Solids. (S. Krishnan et al., 1995). The main implementation file for the Boolean evaluation algorithm is boolean.cpp.

The NURBS surface-surface intersection algorithm is based on the "NURBS Intersection Curve Evaluation" section of the paper Performing Efficient NURBS Modeling Operations on the GPU (A. Krishnamurthy et al., 2008). A detailed outline of the algorithm, as implemented, appears in the main implementation file for the NURBS Surface-Surface intersection algorithm, intersect.cpp.

The ON_Boolean() function performs a single Boolean evaluation on two B-Rep objects. A single execution of the brep command in MGED or Archer may involve passing several successive pairs of B-Rep objects to this function.

	
int
ON_Boolean(
    ON_Brep *evaluated_brep,
    const ON_Brep *brep1,
    const ON_Brep *brep2,
    op_type operation);
	
	

In the nontrivial case where the bounding boxes of brep1 and brep2 intersect, get_evaluated_faces() is called to get the trimmed NURBS faces of the evaluated Boolean result. The faces are then combined into a single B-Rep object returned via the evaluated_brep argument.

	
ON_ClassArray< ON_SimpleArray<Trimmed Face *> >
get_evaluated_faces(
    const ON_Brep *brep1,
    const ON_Brep *brep2,
    op_type operation);
	
	

The intersection curves between the faces of brep1 and brep2 are found by get_face_intersection_curves(). These curves are used to split the original surfaces into pieces, each becoming a new trimmed NURBS face. The categorize_trimmed_faces() function is used to identify which pieces, based on the Boolean operation, are part of the evaluated result. Each TrimmedFace whose m_belong_to_final member is marked TrimmedFace::BELONG is used by ON_Boolean() to create the final evaluated result.

	
ON_ClassArray< ON_SimpleArray<SSICurve> >
get_face_intersection_curves(
    ON_SimpleArray<Subsurface *> &surf_tree1,
    ON_SimpleArray<Subsurface *> &surf_tree2,
    const ON_Brep *brep1,
    const ON_Brep *brep2,
    op_type operation);
	
	

Each pair of brep1 and brep2 surfaces whose bounding boxes intersect are passed to the ON_Intersect() surface-surface intersection routine. The get_subcurves_inside_faces() routine is used to remove irrelevant parts of the surface-surface intersection curves based on the trimming curves of the associated faces.

	
int
ON_Intersect(const ON_Surface *surfA,
             const ON_Surface *surfB,
             ON_ClassArray<ON_SSX_EVENT> &x,
             double isect_tol,
             double overlap_tol,
             double fitting_tol,
             const ON_Interval *surfaceA_udomain,
             const ON_Interval *surfaceA_vdomain,
             const ON_Interval *surfaceB_udomain,
             const ON_Interval *surfaceB_vdomain,
             Subsurface *treeA,
             Subsurface *treeB);
	
	

The first stage of the surface-surface intersection algorithm attempts to identify overlap intersections (areas where the two surfaces are coincident). Our assumption is that the boundary curve of any overlap region must be formed from isocurves of the overlapping surfaces.

Subcurves of isocurves that intersect both surfaces, such that the surfaces are coincident on one side of the curve but not the other, potentially form part of overlap boundaries. These curves are identified using find_overlap_boundary_curves(). To avoid wasted computations, this function also returns intersection points and non-boundary intersection curves which were found during the search for boundary curves.

Then, the split_overlaps_at_intersections() function is run, and curve pieces that share endpoints are stitched together. The stitched boundary curves which close to form loops are recorded as overlap intersection events.

The second stage of the surface-surface intersection algorithm attempts to identify other intersection curves and points. The input surfaces surfA and surfB are subdivided into four subsurfaces, whose bounding boxes are tested in pairs to see which subsurfaces potentially intersect. This subdivision repeats to a fixed depth determined by the constant MAX_SSI_DEPTH (defined in brep_defines.h).

Subsurfaces that lie completely inside an overlap region identified in the first stage are discarded. Each remaining pair of subsurfaces with intersecting bounding boxes is tested for intersection. This is accomplished by approximating each subsurface with two triangles (i.e. a 'split' quad whose corners coincide with those of the actual subsurface patch, which has been split diagonally for a more accurate fit). The triangles are then intersected, and the average of all intersection points is used as the initial guess for a Newton iterative solver, implemented by newton_ssi(), which searches for a point close to the initial guess point which lies on both surfaces.

Solved points that reside inside an overlap region identified in the first stage are discarded. Of the remaining solved intersection points between surfA and surfB, those which are near one another are stitched together into polyline curves. If a line or conic curve can be fit to the polyline curves in 2D, the fit curve replaces the original surfA and/or surfB polyline curve.

OpenNURBS includes two general array classes, ON_ClassArray and ON_SimpleArray, which are similar to C++'s std::vector. Besides having slightly friendlier interfaces, they also feature some higher-level member functions like Reverse() and Quicksort().

The primary difference between the two classes is that ON_SimpleArray doesn't bother constructing and destructing its items. This makes it more efficient than ON_ClassArray, but unsuitable for class objects (though pointers to objects are fine). ON_ClassArray requires items to have correctly implemented copy and assignment functions.

The NURBS Boolean evaluation implementation generally employs a combined array of known size to index elements from two input objects. For example, if brepA has i faces and brepB has j faces, a single array of i + j elements is created.

	
	ON_ClassArray< ON_SimpleArray<SSICurve> > curves_array(face_count1 + face_count2);
	curves_array.SetCount(curves_array.Capacity());
	
	  

Curves and surfaces are nearly always allocated on the heap and referenced by pointers, both in the OpenNURBS library, and in the NURBS Boolean evaluation implementation.

Mostly these allocations are simply done with the new keyword as with any other class. However, a few classes, notably ON_Brep have a New() function that wraps the allocation, which is preferred over using new directly for technical reasons specified in the OpenNURBS opennurbs_brep.h header.

Pointers to objects, curves in particular, are generally "stolen" to avoid having to create a new copy of the object.

	      
ON_SimpleArray<ON_SSX_EVENT> x;
...
for (int i = 0; i < csx_events.Count(); ++i) {
    // copy event
    x.Append(csx_events[i]);

    // clear pointers from original so they aren't deleted by the
    // ON_SSX_EVENT destructor
    csx_events[i].m_curveA = NULL;
    csx_events[i].m_curveB = NULL;
    csx_events[i].m_curve3d = NULL;
}
	      
	    

The OpenNURBS routines make extensive use of the symbol ON_ZERO_TOLERANCE in calculations to test if a result is to be considered equal to zero, or if two values are to be considered equal.

The ON_2dPoint and ON_3dPoint classes intuitively implement operators such as + and * to allow points to be easily summed and scaled.

The operator[] functions are notable because coordinates are not actually stored as arrays in these classes, but rather in the named members x, y, and z. So while accessing coordinates as pt[0], pt[1] is possible, the more readable pt.x, pt.y, is more typically utilized.

The most frequently used member function is DistanceTo(const ON_3dPoint &p), used to check inter-point distances, either as part of an intersection test or to identify closeable gaps or duplicate points.

ON_BoundingBox is returned by the BoundingBox(), GetTightBoundingBox(), and GetBBox() functions, which are implemented by all geometry classes inheriting from ON_Geometry.

The most commonly used members of ON_BoundingBox are Diagonal() (usually in an expression such as bbox.Diagonal().Length() used as a scalar size estimate), and IsPointIn() and MinimumDistanceTo() (used in intersection tests).

ON_Interval is used to represent the domains of parametric curves and surfaces. The domain starts at m_t[0] and ends at m_t[1]. These members can be set directly or via Set(double t0, double t1).

The ParameterAt(double x) function translates a normalized parameter (from a domain starting at 0.0 and ending at 1.0) into a real parameter. Thus, the start of the domain is at domain.ParameterAt(0.0), the midpoint is at domain.ParameterAt(.5), etc.

The most frequently used geometry class is ON_Curve, a generic container for parametric curves. The curve is interrogated by using the PointAt(double t) method to evaluate points at arbitrary values inside the curve's domain, which is specified by the ON_Interval returned by the Domain() method. The start and end points of the curve have dedicated access methods, PointAtStart() and PointAtEnd().

All the PointAt() methods return an ON_3dPoint, though in the common case where ON_Curve objects are representing 2D trim curves, the z coordinate will be 0.0.

It is sometimes necessary to reverse a curve's domain. This is done using the Reverse() method to facilitate stitching curves together. The function has a Boolean int return value that must be checked.

if (curveA->PointAtStart().DistanceTo(curveB->PointAtStart()) < dist_tol) {
  if (curveA->Reverse()) {
      curveA = link_curves(curveA, curveB);
  }
  /* curves that cannot be reversed are degenerate and discarded */
}

	  

A copy of a curve is easily made using the Duplicate() member function, which simply wraps a standard copy procedure:

ON_Curve* Duplicate()
{
  ON_Curve *p = new ON_Curve;
  if (p) *p = *this;
  return p;
}
	

This function is common to all OpenNURBS geometry classes, but curves are by far the most frequently duplicated objects. However, if curves are simply being retained from a working set of container objects, the curve pointers are generally "stolen" rather than copied, with curve members set to NULL so that the curves aren't destructed with the containers.

ON_Line is used to represent an infinite line, defined by two points, from and to.

ON_Line is not a subclass of ON_Curve and should not be confused with ON_LineCurve (which has an ON_Line member), though it does have some of the same methods as an ON_Curve class, including PointAt(double t). However, because the line has an infinite domain, it can be evaluated at any t value, though evaluating at 0.0 returns from and evaluating at 1.0 returns to, as if the line was a parametric curve with a domain between 0.0 and 1.0.

ON_Line has helpful line-specific methods such as ClosestPointTo(const ON_3dPoint &point). Again, because the line is treated as infinite, this function doesn't necessarily return a point in the segment between from and to.

An ON_Surface has a similar interface to an ON_Curve, but adapted to support the surface's two domains, u and v (sometimes called s and t). These also correspond to as the 0 and 1 surface domains (as in the first example in following) or with an x and y parameterization (as shown in the second example).

ON_Interval udom = surface->Domain(0);
ON_Interval vdom = surface->Domain(1);
ON_3dPoint surf_midpt_3d = surface->PointAt(udom.ParameterAt(.5), vdom.ParameterAt(.5));

	  

ON_Interval tdom = trim_curve->Domain();
ON_3dPoint trim_midpt_uv = trim_curve->PointAt(tdom.ParameterAt(.5));
ON_3dPoint trim_midpt_3d = surface->PointAt(trim_midpt_uv.x, trim_midpt_uv.y);

	  

ON_Brep is the top-level OpenNURBS class used to represent the two input objects and the evaluated result of the ON_Boolean() function. The geometry is encoded as a collection of faces, which for our purposes should be topologically connected to enclose solid volumes.

An object's faces are ON_BrepFace objects stored in the ON_Brep face array, m_F[].

Each ON_BrepFace is defined as the subset of an ON_Surface lying inside the face's outerloop (a.k.a. the face boundary) and outside all of its innerloops (a.k.a. trim loops or just trims).

The loops of an ON_BrepFace are listed in its loop array m_li[] as indexes into the associated ON_Brep object's ON_BrepLoop array, m_L[]. The first (and possibly only) loop listed in the face's loop index array is the outerloop, and all following loops are inner trim loops. The type of the loop is also recorded in the loop's m_type member.

brep->m_L[brep->m_F[0]->m_li[0]].m_type;      // ON_BrepLoop::outer
brep->m_L[brep->m_F[0]->m_li[1]].m_type;      // ON_BrepLoop::inner
...
brep->m_L[*brep->m_F[0]->m_li.Last()].m_type; // ON_BrepLoop::inner
	

There are two OpenNURBS classes for representing intersections. ON_X_EVENT is used for curve-curve and curve-surface intersections. ON_SSX_EVENT is used for surface-surface intersections.

The intersection classes enumerate a number of intersection types. Over the course of an evaluation, the m_type of intersection events is repeatedly checked to determine how each event should be processed.

When two curves are coincident with one another over a portion of their domains, m_type will be ON_X_EVENT::ccx_overlap.

Two curves (blue on the left and red on the right) overlap over a portion of their domains (white).

Figure 2. Curve-Curve Overlap Intersection

When two surfaces are coincident over a portion of their domains, m_type will be ON_SSX_EVENT::ssx_overlap.

Two surfaces (blue on the left and red on the right) overlap over a region of their domains (green).

Figure 3. Surface-Surface Overlap Intersection

There are two ways that two surfaces can intersect in a curve. If the normals of the surfaces are parallel over all points of the curve, the intersection m_type is ON_SSX_EVENT::ssx_tangent, and ON_SSX_EVENT::ssx_transverse otherwise.

Two surfaces (red and blue) intersect in a curve (white). The normals of both surfaces are parallel to each other everywhere along the curve.

Figure 4. Surface-Surface Tangent Intersection

Two surfaces (red and blue) intersect in a curve (yellow). The normals of both surfaces are not parallel to each other along the curve.

Figure 5. Surface-Surface Transverse Intersection

Similarly, if two surfaces intersect at a point, the intersection m_type is ON_SSX_EVENT::ssx_tangent_point if the normals of the two surfaces are parallel at that point, and ON_SSX_EVENT::ssx_transverse_point otherwise.

The m_type of an intersection event determines how values in the m_a[], m_b[], m_A[], and m_B[] array members of the event instance are to be interpreted (documented in the OpenNURBS opennurbs_x.h header).

For an ON_X_EVENT representing a curve-curve intersection whose m_type is ON_X_EVENT::ccx_overlap, (m_a[0], m_a[1]) represents the portion of the first curve's domain that overlaps with the second curve, whereas in other cases m_a[1] is simply a duplicate of m_a[0].

As a result, a pattern seen repeatedly in the NURBS Boolean evaluation implementation is a loop over intersection events that gathers intersection points for processing, including overlap endpoints if the event represents an overlap.

            
for (int i = 0; i < x_event.Count(); ++i) {
    x_points.Append(x_event[i].m_a[0]);
    if (x_event[i].m_type == ON_X_EVENT::ccx_overlap) {
        x_points.Append(x_event[i].m_a[1]);
    }
}
            
	

Implicit in working with parametric geometry is that some operations are done in 2D while others are done in 3D and it's very important to know the dimension currently being worked in at all times.

As mentioned in the section above on 2D and 3D points, 3D classes are often used in the implementation to store 2D points, and thus are not a reliable indication that an operation is happening in 3D.

Being that operations in 2D tend to be a lot simpler, 2D is normally the dimension being worked in. However, because parametric curves and surfaces of different objects have different parameterizations, determining where two objects intersect can't be done by comparing 2D parameters; it must happen in 3D.

3.5.1.3. Curve Traversal Directions

It's important to remember that because parameterizations are arbitrary, there is no correspondence whatsoever between a 2D curve in one surface's domain and another surface's domain, even when those 2D curves evaluate to the same 3D curve. In particular, you cannot assume that traversing different curves along their domain from m_t[0] to m_t[1] translates to a consistent traversal direction in 3D, or even that each 2D curve's m_t[0]/m_t[1] corresponds to the same 3D starting point on a closed curve.

The projections (red and blue) of two 2D curves in different domains represent the same 3D curve (within tolerance). However, the start and end points as well as the 3D traversal directions differ between the projections.

Figure 6. Different Traversals of the Same Curve

By the nature of the math involved in representing parametric geometry (e.g. converting between 2D and 3D, and solving intersections between objects with different parameterizations) values that are expected to be identical are generally only equal within a certain tolerance, or error.

Over the course of the evaluation, the same data is interrogated and processed a number of times. If ignored, the error introduced in one stage of the evaluation can grow over subsequent stages, causing an incorrect determination that leaves a curve unclosed, a surface unsplit, and ultimately an incorrect evaluated result.

As a consequence, it's generally a good idea to remove fuzziness when you find it, and avoid algorithms that introduce more error.

The Archer brep command (also implemented in MGED) can be used to get structural information about B-Rep objects and visualize different subcomponents.

Most importantly, brep <obj> info will report summary information, including the number of NURBS surfaces and faces and brep <obj> plot S <index> can be used to plot individual surfaces in 3D.

This is the primary way you can conceptually link a surface or face index to the 3D geometry it represents. So if you notice an error in an object while viewing it in the editor, you can use the brep command to determine the index of the surface with the error, and then inspect the in-memory object in a debugger using that index into the final surface array, tracing that surface object to where it was created, etc.

The dplot command is used to visualize the results of different stages of the NURBS Boolean evaluation algorithm. This makes it easier to isolate the source of a problem in an evaluation.

Unlike the brep command, the dplot command is purely a development tool. Its implementation does not honor library boundaries and does not conform to the typical conventions for editor commands, and for this reason is only available as an Archer command in the NURBS Boolean evaluation development branch (https://sourceforge.net/p/brlcad/code/HEAD/tree/brlcad/branches/brep-debug/).

In the development branch, the NURBS Boolean evaluation source code contains additional calls to DebugPlot functions (implemented in debug_plot.cpp) that create wireframe visualizations of data produced during the evaluations.

For development convenience, these wireframes are not saved as database objects, but rather are written as files in the current directory, with names of the form bool1_*.plot3. An additional bool1.dplot is written which describes the .plot3 files that were written in a format understood by the dplot command.

One set of files is written for each evaluation. Between evaluations, a static counter increments the numeric suffix that's used in the output filenames. So for a combination consisting of three objects, the bool1* files will hold results from the intermediate boolean evaluation between the first two objects in the combination, and the bool2* files will hold results from the final evaluation between the intermediate evaluated object and the remaining leaf of the original comb.

The DebugPlot functions always use the same file names and do not check if written files already exist. It is assumed that you will run an evaluation, inspect the generated files using the dplot command, and then manually remove (or just move) the generated .dplot and .plot3 files before performing another evaluation with the brep command.

It's possible to write out custom curves from any part of the evaluation (i.e. those not covered by dplot) and view them in MGED/Archer.

You can pass a 3D ON_Curve to the DebugPlot::Plot3DCurve() function or a 2D ON_Curve and an associated ON_Surface to the DebugPlot::Plot3DCurve() function.

Both of these functions take an arbitrary filename for a plot3 file the function will write, as well as a color for the curve. The DebugPlot::Plot3DCurve() has an optional vlist parameter which you should omit (see the full definitions in debug_plot.cpp).

            
// somewhere in boolean.cpp
if (face1_curves.Count() > 0 && face1_curves[0] != NULL) {
    static int calls = 0;
    unsigned char mycolor[] = {0, 0, 62};
    std::ostringstream plotname;

    // generate a unique filename
    plotname << "mycurve" << ++calls << ".plot3";

    // plot using method of global DebugPlot instance 'dplot'
    dplot->Plot3DCurveFrom2D(surf1, face1_curves[0],
        plotname.str().c_str(), mycolor);
}
            
	  

After running an evaluation that produces a custom plot3 file, you can draw it using the overlay editor command.

Archer> overlay mycurve1.plot3 1
          

After you notice a problem in the output shown by the dplot command, you need to locate the source code that created the erroneous geometry so you can start debugging. The following sections provide example procedures you can perform in Archer and a debugger to start investigating some common issues.

If the ssx subcommand shows that a surface-surface intersection is missing...

If the isocsx subcommand shows that an isocurve-surface intersection is missing...

If the isocsx subcommand shows that isocurve intersections are incorrect...

If the ssx subcommand shows an incorrect intersection curve...

If the ssx subcommand shows the correct intersections for a given surface pair, but the fcurves subcommand shows those curves are not being correctly clipped by faces...

If the faces subcommand shows that an input face was not split correctly, but the lcurves subcommand shows the relevant intersection was accurate...

What follows is a step-by-step debugging of a real issue affecting the X combination from the BRL-CAD sample database axis.g.

This issue was fixed in revision 65179 in the NURBS Boolean evaluation development branch of the source repository (https://sourceforge.net/p/brlcad/code/HEAD/tree/brlcad/branches/brep-debug/).

If you want to follow along, you can reinstate the error in a checkout of the development branch:

$ svn merge -r 65179:65178 ^/brlcad/branches/brep-debug