Cull BREP tracing work and report telemetry

Reuse one de Casteljau triangle for one-axis restrictions touching a
parameter endpoint, and add exact high-order enclosure and performance
gates.

include/brep.h

@@ -44,6 +44,7 @@ __BEGIN_DECLS
 #include "brep/quality.h"
 #include "brep/intersect.h"
 #include "brep/boolean.h"
+#include "brep/cobb.h"
 #include "brep/cdt.h"
 #include "brep/csg.h"
 

include/brep/CMakeLists.txt

@@ -4,6 +4,7 @@ set(
   bbnode.h
   assembly.h
   boolean.h
+  cobb.h
   brnode.h
   cdt.h
   csg.h

include/brep/bbnode.h

@@ -29,6 +29,13 @@
 
 #include "common.h"
 
+#ifdef __cplusplus
+extern "C++" {
+#  include <cstdint>
+#  include <list>
+}
+#endif
+
 #include "brep/defines.h"
 #include "brep/curvetree.h"
 #include "brep/ray.h"
@@ -41,6 +48,7 @@
 extern "C++" {
 namespace brlcad {
     class BBNode;
+    class BBNodeArena;
 }
 }
 #endif
@@ -65,7 +73,7 @@ namespace brlcad {
     /**
      * Bounding Box Hierarchy Node
      */
-    class BREP_EXPORT BBNode : public PooledObject<BBNode> {
+    class BREP_EXPORT BBNode {
     public:
 	explicit BBNode(const ON_BoundingBox &node, const CurveTree *ct = NULL);
 	BBNode(const CurveTree *ct,
@@ -77,6 +85,8 @@ namespace brlcad {
 	~BBNode();
 
 	BBNode(Deserializer &deserieralizer, const CurveTree &ctree);
+	static BBNode *deserializeTree(Deserializer &deserializer,
+		const CurveTree &ctree);
 	void serialize(Serializer &serializer) const;
 
 	const ON_BrepFace &get_face() const;
@@ -115,18 +125,36 @@ namespace brlcad {
 	/** Test whether a ray intersects the 3D bounding volume of the
 	 * node - if so, and node is not a leaf node, query children.  If
 	 * leaf node, and intersects, add to list.
+	 *
+	 * The fixed-array overload appends while count is below capacity,
+	 * increments count for every intersected leaf, and sets overflow if
+	 * any leaf cannot be stored.  The caller initializes count and
+	 * overflow and may use the total count to select a fallback.
 	 */
 	bool intersectsHierarchy(const ON_Ray &ray, std::list<const BBNode *> &results) const;
+	bool intersectsHierarchy(const ON_Ray &ray, const BBNode **results,
+		std::size_t capacity, std::size_t &count,
+		bool &overflow) const;
 
 	ON_2dPoint getClosestPointEstimate(const ON_3dPoint &pt) const;
 	ON_2dPoint getClosestPointEstimate(const ON_3dPoint &pt, ON_Interval &u, ON_Interval &v) const;
 	int getLeavesBoundingPoint(const ON_3dPoint &pt, std::list<const BBNode *> &out) const;
-	bool isTrimmed(const ON_2dPoint &uv, const BRNode **closest, double &closesttrim, double within_distance_tol) const;
+	/** Classify against the prepared trim leaves without allocating.
+	 * isTrimmedAllocating retains the former filtered-list path for
+	 * transition-time equivalence diagnostics only. */
+	bool isTrimmed(const ON_2dPoint &uv, const BRNode **closest,
+		double &closesttrim, double within_distance_tol,
+		std::size_t *candidate_count = NULL) const;
+	bool isTrimmedAllocating(const ON_2dPoint &uv, const BRNode **closest,
+		double &closesttrim, double within_distance_tol,
+		std::size_t *candidate_count = NULL) const;
 
 	void BuildBBox();
-	bool prepTrims();
+	bool prepTrims(double within_distance_tol =
+		BREP_EDGE_MISS_TOLERANCE);
 
-	const std::vector<BBNode *> &get_children() const;
+	std::size_t childCount() const;
+	BBNode *child(std::size_t index) const;
 
 	/** Bounding Box */
 	ON_BoundingBox m_node;
@@ -138,6 +166,8 @@ namespace brlcad {
 	/** Trimming Flags */
 	bool m_checkTrim;
 	bool m_trimmed;
+	/* Arena construction nodes are destroyed shallowly by their arena. */
+	bool m_arena_node;
 
 	/** Point used for closeness testing - usually based on evaluation
 	 * of the curve/surface at the center of the parametric domain
@@ -152,10 +182,15 @@ namespace brlcad {
 
 
     private:
+	friend class BBNodeArena;
+
 	BBNode(const BBNode &source);
 	BBNode &operator=(const BBNode &source);
+	BBNode(Deserializer &deserializer, const CurveTree &ctree,
+		BBNodeArena *arena);
+	void adopt(BBNode &source);
+	void setArenaOwner(BBNodeArena *owner);
 
-	void removeChild(BBNode *child);
 	bool intersectedBy(const ON_Ray &ray, double *tnear = NULL, double *tfar = NULL) const;
 
 	/** Report if a given uv point is within the uv boundaries defined
@@ -163,39 +198,41 @@ namespace brlcad {
 	 */
 	bool containsUV(const ON_2dPoint &uv) const;
 
-	void getTrimsAbove(const ON_2dPoint &uv, std::list<const BRNode *> &out_leaves) const;
 	const BBNode *closer(const ON_3dPoint &pt, const BBNode *left, const BBNode *right) const;
 
-	struct Stl : public PooledObject<Stl> {
-	Stl() : m_children(), m_trims_above() {}
-
-	    std::vector<BBNode *> m_children;
-	    std::list<const BRNode *> m_trims_above;
-	} * const m_stl;
+	BBNode **m_children;
+	uint32_t m_child_count;
+	uint32_t m_child_capacity;
     };
 
-    inline void
-    BBNode::addChild(BBNode *child)
-    {
-	if (LIKELY(child != NULL)) {
-	    m_stl->m_children.push_back(child);
-	}
-    }
+    /**
+     * Per-SurfaceTree construction storage.  The first node uses ordinary
+     * transferable storage and subsequent nodes are bulk allocated.  The
+     * final logical root is adopted into that first node, whose lifetime owns
+     * the arena.
+     */
+    class BBNodeArena {
+    public:
+	BBNodeArena();
+	~BBNodeArena();
 
-    inline void
-    BBNode::removeChild(BBNode *child)
-    {
-	std::vector<BBNode *>::iterator i;
-	for (i = m_stl->m_children.begin(); i != m_stl->m_children.end();) {
-	    if (*i == child) {
-		delete *i;
-		i = m_stl->m_children.erase(i);
-	    } else {
-		++i;
-	    }
-	}
-    }
+	BBNode *make(const ON_BoundingBox &node,
+		const CurveTree *ct = NULL);
+	BBNode *make(const CurveTree *ct, const ON_BoundingBox &node,
+		const ON_Interval &u, const ON_Interval &v,
+		bool checkTrim, bool trimmed);
+	BBNode *finish(BBNode *root);
 
+    private:
+	friend class BBNode;
+	struct Block;
+	void *allocate();
+	BBNode *make(Deserializer &deserializer, const CurveTree &ctree);
+
+	Block *m_blocks;
+	BBNode *m_anchor;
+	bool m_started;
+    };
 
     inline const ON_BrepFace &
     BBNode::get_face() const
@@ -203,16 +240,22 @@ namespace brlcad {
 	return *m_ctree->m_face;
     }
 
-    inline const std::vector<BBNode *> &
-    BBNode::get_children() const
+    inline std::size_t
+    BBNode::childCount() const
     {
-	return m_stl->m_children;
+	return m_child_count;
+    }
+
+    inline BBNode *
+    BBNode::child(std::size_t index) const
+    {
+	return index < m_child_count ? m_children[index] : NULL;
     }
 
     inline bool
     BBNode::isLeaf() const
     {
-	return m_stl->m_children.empty();
+	return m_child_count == 0;
     }
 
     inline void

include/brep/brnode.h

@@ -33,8 +33,7 @@
 #ifdef __cplusplus
 extern "C++" {
 /* @cond */
-#  include <vector>
-#  include <list>
+#  include <algorithm>
 /* @endcond */
 }
 #endif
@@ -49,34 +48,19 @@ extern "C++" {
 namespace brlcad {
 
     /**
-     * Bounding Rectangle Hierarchy
+     * Compact prepared trim-curve leaf.
      */
-    class BREP_EXPORT BRNode : public PooledObject<BRNode> {
+    class BREP_EXPORT BRNode {
     public:
-	explicit BRNode(const ON_BoundingBox &node);
 	BRNode(const ON_Curve *curve,
 	       int trim_index,
 	       int adj_face_index,
 	       const ON_BoundingBox &node,
-	       const ON_BrepFace *face,
 	       const ON_Interval &t,
-	       bool innerTrim,
-	       bool checkTrim,
-	       bool trimmed);
-	~BRNode();
-
+	       bool innerTrim);
 	BRNode(Deserializer &deserializer, const ON_Brep &brep);
 	void serialize(Serializer &serializer) const;
 
-	/** Node management functions */
-	void addChild(BRNode *child);
-
-	/** Return a list of all nodes below this node that are leaf nodes */
-	void getLeaves(std::list<const BRNode *> &out_leaves) const;
-
-	/** Report the depth of this node in the hierarchy */
-	int depth() const;
-
 	/**
 	 * Get 2 points defining bounding box:
 	 *
@@ -91,19 +75,42 @@ namespace brlcad {
 	 */
 	void GetBBox(fastf_t *min, fastf_t *max) const;
 
-	bool isTrimmed(const ON_2dPoint &uv, double &trimdist) const;
+	/** Test a vertical +V parameter-space ray against this monotone
+	 * trim leaf.  Endpoint fuzz bridges small imported loop gaps while
+	 * the exact test uses a half-open interval to avoid double counting
+	 * shared vertices. */
+	bool crossesAbove(const ON_2dPoint &uv, double &distance,
+		double endpoint_tolerance = 0.0) const;
+	/** Test a horizontal +U parameter-space ray using the same half-open
+	 * endpoint ownership as crossesAbove(). */
+	bool crossesRight(const ON_2dPoint &uv, double &distance,
+		double endpoint_tolerance = 0.0) const;
+	/** Return a fixed-workspace estimate of the shortest UV distance from
+	 * a point to this trim leaf. */
+	double curveDistance(const ON_2dPoint &uv) const;
+	/** Lower bound from the leaf's sampled UV bounding box. */
+	double curveBBoxDistance(const ON_2dPoint &uv) const;
+	int trimIndex() const { return m_trim_index; }
+	ON_Interval parameterInterval() const { return m_t; }
+	ON_Interval curveDomain() const {
+	    return m_trim ? m_trim->Domain() : ON_Interval::EmptyInterval;
+	}
+	ON_3dPoint startPoint() const {
+	    return ON_3dPoint(m_start[0], m_start[1], 0.0);
+	}
+	ON_3dPoint endPoint() const {
+	    return ON_3dPoint(m_end[0], m_end[1], 0.0);
+	}
+	ON_2dPoint bboxMinimum() const {
+	    return ON_2dPoint(m_bbox_min[0], m_bbox_min[1]);
+	}
+	ON_2dPoint bboxMaximum() const {
+	    return ON_2dPoint(m_bbox_max[0], m_bbox_max[1]);
+	}
 
-	ON_2dPoint getClosestPointEstimate(const ON_3dPoint &pt) const;
-	ON_2dPoint getClosestPointEstimate(const ON_3dPoint &pt, ON_Interval &u, ON_Interval &v) const;
 	fastf_t getCurveEstimateOfV(fastf_t u, fastf_t tol) const;
 	fastf_t getCurveEstimateOfU(fastf_t v, fastf_t tol) const;
 
-	/** Bounding Box */
-	ON_BoundingBox m_node;
-
-	/** Surface Information */
-	ON_Interval m_v;
-
 	/** Trim Curve Information */
 	int m_adj_face_index;
 
@@ -114,85 +121,23 @@ namespace brlcad {
 	bool m_innerTrim;
 
     private:
-	BRNode(const BRNode &source);
-	BRNode &operator=(const BRNode &source);
-
-	void removeChild(BRNode *child);
-
-	/** Test if this node is a leaf node (i.e. m_children is empty) */
-	bool isLeaf() const;
-
-	fastf_t getLinearEstimateOfV(fastf_t u) const;
-	const BRNode *closer(const ON_3dPoint &pt, const BRNode *left, const BRNode *right) const;
-
-	struct Stl : public PooledObject<Stl> {
-	Stl() : m_children() {}
-
-	    std::vector<const BRNode *> m_children;
-	} * const m_stl;
-
-	const ON_BrepFace *m_face;
-	ON_Interval m_u;
 	const ON_Curve *m_trim;
 	int m_trim_index;
 	ON_Interval m_t;
-	bool m_checkTrim;
-	bool m_trimmed;
-	ON_3dPoint m_estimate;
-	fastf_t m_slope;
-	fastf_t m_bb_diag;
-	ON_3dPoint m_start;
-	ON_3dPoint m_end;
+	double m_bbox_min[2];
+	double m_bbox_max[2];
+	double m_start[2];
+	double m_end[2];
     };
 
-    inline void
-    BRNode::addChild(BRNode *child)
-    {
-	if (LIKELY(child != NULL)) {
-	    m_stl->m_children.push_back(child);
-	}
-    }
-
-    inline void
-    BRNode::removeChild(BRNode *child)
-    {
-	std::vector<const BRNode *>::iterator i;
-	for (i = m_stl->m_children.begin(); i != m_stl->m_children.end();) {
-	    if (*i == child) {
-		delete *i;
-		i = m_stl->m_children.erase(i);
-	    } else {
-		++i;
-	    }
-	}
-    }
-
-    inline bool
-    BRNode::isLeaf() const
-    {
-	if (m_stl->m_children.empty()) {
-	    return true;
-	}
-	return false;
-    }
-
     inline void
     BRNode::GetBBox(fastf_t *min, fastf_t *max) const
     {
-	VSETALL(min, INFINITY);
-	VSETALL(max, -INFINITY);
-	if (m_start != ON_3dPoint::UnsetPoint) {
-	    VMINMAX(min, max, m_start);
-	}
-	if (m_end != ON_3dPoint::UnsetPoint) {
-	    VMINMAX(min, max, m_end);
-	}
+	VSET(min, std::min(m_start[X], m_end[X]),
+	    std::min(m_start[Y], m_end[Y]), 0.0);
+	VSET(max, std::max(m_start[X], m_end[X]),
+	    std::max(m_start[Y], m_end[Y]), 0.0);
     }
-
-    extern bool sortX(const BRNode *first, const BRNode *second);
-    extern bool sortY(const BRNode *first, const BRNode *second);
-
-
 } /* namespace brlcad */
 } /* extern C++ */
 

include/brep/cdt.h

@@ -159,6 +159,9 @@ extern BREP_EXPORT int cdt_bmesh_create(struct cdt_bmesh **m);
 extern BREP_EXPORT void cdt_bmesh_destroy(struct cdt_bmesh *m);
 extern BREP_EXPORT int cdt_bmesh_deserialize(const char *fname, struct cdt_bmesh *m);
 extern BREP_EXPORT int cdt_bmesh_repair(struct cdt_bmesh *m);
+extern BREP_EXPORT int cdt_test_boundary_start(void);
+extern BREP_EXPORT int cdt_test_boundary_steiner_filter(void);
+extern BREP_EXPORT int cdt_test_spurious_components(void);
 
 __END_DECLS
 

include/brep/cobb.h

@@ -0,0 +1,65 @@
+/*                         C O B B . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2026 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ */
+/** @file brep/cobb.h
+ *
+ * Reusable Cobb rational Bezier sphere fixtures.
+ */
+#ifndef BREP_COBB_H
+#define BREP_COBB_H
+
+#include "common.h"
+
+#include "brep/defines.h"
+
+#ifdef __cplusplus
+extern "C++" {
+
+class ON_3dPoint;
+class ON_BezierSurface;
+class ON_Brep;
+
+/** Return one unit Cobb sphere patch after rotations in degrees. */
+BREP_EXPORT extern ON_BezierSurface *
+ON_Brep_CobbSphereFace(double rotation_x, double rotation_z);
+
+/**
+ * Return the legacy six-patch representation with independent face topology.
+ * This is an open plate-mode fixture, not a solid.
+ */
+BREP_EXPORT extern ON_Brep *
+ON_Brep_CobbSphereUnsewn(double radius, const ON_3dPoint &origin);
+
+/**
+ * Return the same six support surfaces with eight shared vertices and twelve
+ * shared edges.  The result is a closed, oriented solid.
+ */
+BREP_EXPORT extern ON_Brep *
+ON_Brep_CobbSphereSewn(double radius, const ON_3dPoint &origin);
+
+/** Measure maximum radial error on a uniform per-surface sample grid. */
+BREP_EXPORT extern double
+ON_Brep_CobbSphereMaxRadialError(const ON_Brep *brep, double radius,
+	const ON_3dPoint &origin, int samples_per_direction = 33);
+
+}
+#endif
+
+#endif  /* BREP_COBB_H */
+
+/*
+ * Local Variables:
+ * mode: C++
+ * tab-width: 8
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */

include/brep/curvetree.h

@@ -32,11 +32,11 @@
 #include "common.h"
 
 #ifdef __cplusplus
-// @cond SKIP_C++_INCLUDE
 extern "C++" {
-#  include <map>
+#  include <cstdint>
+#  include <list>
+#  include <vector>
 }
-// @endcond
 #endif
 
 #include "brep/defines.h"
@@ -54,24 +54,39 @@ namespace brlcad {
     /**
      * CurveTree declaration
      */
-    class BREP_EXPORT CurveTree : public PooledObject<CurveTree> {
+    class BREP_EXPORT CurveTree {
     public:
 	explicit CurveTree(const ON_BrepFace *face);
 	~CurveTree();
 
 	CurveTree(Deserializer &deserializer, const ON_BrepFace &face);
 	void serialize(Serializer &serializer) const;
-	std::vector<std::size_t> serialize_get_leaves_keys(const std::list<const BRNode *> &leaves) const;
-	std::list<const BRNode *> serialize_get_leaves(const std::size_t *keys, std::size_t num_keys) const;
-	void serialize_cleanup() const;
 
 	/**
 	 * Return just the leaves of the surface tree
 	 */
 	void getLeaves(std::list<const BRNode *> &out_leaves) const;
-	void getLeavesAbove(std::list<const BRNode *> &out_leaves, const ON_Interval &u, const ON_Interval &v) const;
+	void getLeavesAbove(std::list<const BRNode *> &out_leaves,
+		const ON_Interval &u, const ON_Interval &v,
+		fastf_t tol = BREP_UV_DIST_FUZZ) const;
 	void getLeavesAbove(std::list<const BRNode *> &out_leaves, const ON_2dPoint &pt, fastf_t tol) const;
 	void getLeavesRight(std::list<const BRNode *> &out_leaves, const ON_2dPoint &pt, fastf_t tol) const;
+	std::size_t preparedLeafCount() const;
+	std::size_t preparedLeafImageCount() const;
+	int preparedMaximumDepth() const;
+
+	/** Classify a face parameter against all prepared trim leaves without
+	 * allocating a filtered list. */
+	bool isTrimmed(const ON_2dPoint &uv, const BRNode **closest,
+		double &closesttrim, double within_distance_tol,
+		std::size_t *candidate_count = NULL,
+		bool heal_join_gaps = true) const;
+
+	/** Certify a UV cell whose material state cannot change without
+	 * crossing a trim or a structural trim endpoint.  Returns false when
+	 * the cell must retain pointwise trim checks. */
+	bool classifyCell(const ON_Interval &u, const ON_Interval &v,
+		double within_distance_tol, bool &trimmed) const;
 
     private:
 	friend class BBNode;
@@ -79,37 +94,50 @@ namespace brlcad {
 	CurveTree(const CurveTree &source);
 	CurveTree &operator=(const CurveTree &source);
 
-	const BRNode *getRootNode() const;
-
-	/**
-	 * Calculate, using the surface bounding volume hierarchy, a uv
-	 * estimate for the closest point on the surface to the point in
-	 * 3-space.
-	 */
-	ON_2dPoint getClosestPointEstimate(const ON_3dPoint &pt) const;
-	ON_2dPoint getClosestPointEstimate(const ON_3dPoint &pt, ON_Interval &u, ON_Interval &v) const;
-
 	void getLeavesRight(std::list<const BRNode *> &out_leaves, const ON_Interval &u, const ON_Interval &v) const;
 
-	int depth() const;
-
 	bool getHVTangents(const ON_Curve *curve, const ON_Interval &t, std::list<fastf_t> &list) const;
 	bool isLinear(const ON_Curve *curve, double min, double max) const;
-	BRNode *subdivideCurve(const ON_Curve *curve, int trim_index, int adj_face_index, double min, double max, bool innerTrim, int depth) const;
-	BRNode *curveBBox(const ON_Curve *curve, int trim_index, int adj_face_index, const ON_Interval &t, bool isLeaf, bool innerTrim, const ON_BoundingBox &bb) const;
-	static ON_BoundingBox initialLoopBBox(const ON_BrepFace &face);
+	void prepareLoopLeaves();
+	void subdivideCurve(const ON_Curve *curve, int trim_index,
+		int adj_face_index, double min, double max, bool innerTrim,
+		int depth);
 
 	const ON_BrepFace * const m_face;
-	BRNode *m_root;
 
+	struct Stl {
+	    struct Loop {
+		struct Leaf {
+		    Leaf() : node_index(0), image{0, 0},
+			extend_minimum{false, false} {}
+		    Leaf(uint32_t leaf_index, int u_image, int v_image,
+			bool extend_u_minimum, bool extend_v_minimum) :
+			node_index(leaf_index), image{u_image, v_image},
+			extend_minimum{extend_u_minimum,
+			    extend_v_minimum} {}
 
-	struct Stl : public PooledObject<Stl> {
-	Stl() : m_sortedX() {}
+		    uint32_t node_index;
+		    int image[2];
+		    bool extend_minimum[2];
+		};
 
-	    std::vector<const BRNode *> m_sortedX;
+		Loop() : index(-1), type(ON_BrepLoop::unknown), minimum(0.0, 0.0),
+		    maximum(0.0, 0.0), have_bbox(false), sorted() {}
+
+		int index;
+		ON_BrepLoop::TYPE type;
+		ON_2dPoint minimum;
+		ON_2dPoint maximum;
+		bool have_bbox;
+		std::vector<Leaf> sorted;
+	    };
+
+	Stl() : m_leaves(), m_loops(), m_maximum_depth(0) {}
+
+	    std::vector<BRNode> m_leaves;
+	    std::vector<Loop> m_loops;
+	    int m_maximum_depth;
 	} * const m_stl;
-
-	mutable std::map<const BRNode *, std::size_t> *m_sortedX_indices;
     };
 
 } /* namespace brlcad */

include/brep/surfacetree.h

@@ -57,6 +57,14 @@ namespace brlcad {
      */
     class BREP_EXPORT SurfaceTree {
     public:
+	enum FailureReason {
+	    FAILURE_NONE = 0,
+	    FAILURE_NULL_SURFACE,
+	    FAILURE_BOUNDING_BOX,
+	    FAILURE_NORMAL_EVALUATION,
+	    FAILURE_SUBDIVISION
+	};
+
 	explicit SurfaceTree(const ON_BrepFace *face, bool removeTrimmed = true, int depthLimit = BREP_MAX_FT_DEPTH, double within_distance_tol = BREP_EDGE_MISS_TOLERANCE);
 	~SurfaceTree();
 
@@ -66,9 +74,29 @@ namespace brlcad {
 		return false;
 	    return true;
 	};
+	/* The first deterministic construction failure and its local cell. */
+	FailureReason Failure() const { return m_failure; };
+	ON_Interval FailureDomain(int direction) const {
+	    return direction == 0 ? m_failure_u : m_failure_v;
+	};
+	int FailureDepth() const { return m_failure_depth; };
 
 	BBNode *getRootNode() const;
 
+	/**
+	 * Transfer ownership of the prepared surface hierarchy to a caller.
+	 * The returned node must be deleted by the caller.  After this call the
+	 * SurfaceTree destructor will not delete the hierarchy.
+	 */
+	BBNode *releaseRootNode();
+
+	/**
+	 * Transfer ownership of the curve hierarchy referenced by the prepared
+	 * BBNodes to a caller.  The curve tree must outlive the BBNode hierarchy.
+	 * After this call the SurfaceTree destructor will not delete it.
+	 */
+	const CurveTree *releaseCurveTree();
+
 	/**
 	 * Calculate, using the surface bounding volume hierarchy, a uv
 	 * estimate for the closest point on the surface to the point in
@@ -100,14 +128,21 @@ namespace brlcad {
 	bool isFlatU(const ON_Plane frames[9]) const;
 	bool isFlatV(const ON_Plane frames[9]) const;
 	bool hasSplit(const ON_Surface *surf, int dir, const ON_Interval &interval, double &split) const;
+	BBNode *initialBBox(const CurveTree *ctree, const ON_Surface *surf, const ON_Interval &u, const ON_Interval &v, int depth) const;
 	BBNode *subdivideSurface(const ON_Surface *localsurf, const ON_Interval &u, const ON_Interval &v, ON_Plane frames[9], int depth, int depthLimit, int prev_knot, double within_distance_tol) const;
-	BBNode *surfaceBBox(const ON_Surface *localsurf, bool leaf, const ON_Plane frames[9], const ON_Interval &u, const ON_Interval &v, double within_distance_tol) const;
+	BBNode *surfaceBBox(const ON_Surface *localsurf, bool leaf, const ON_Plane frames[9], const ON_Interval &u, const ON_Interval &v, int depth, double within_distance_tol) const;
+	void recordFailure(FailureReason reason, const ON_Interval &u, const ON_Interval &v, int depth) const;
 
 	const bool m_removeTrimmed;
 	const ON_BrepFace * const m_face;
 	BBNode *m_root;
+	BBNodeArena *m_node_arena;
 	std::queue<ON_Plane *> * const m_f_queue;
 	SurfaceSplitCache * const m_split_cache;
+	mutable FailureReason m_failure;
+	mutable ON_Interval m_failure_u;
+	mutable ON_Interval m_failure_v;
+	mutable int m_failure_depth;
     };
 
 } /* namespace brlcad */

include/rt/rt_instance.h

@@ -219,8 +219,8 @@ RT_EXPORT extern int rt_gettrees(struct rt_i *rtip,
  * observing nodes or soltabs that the finishing pass will subsequently free.
  *
  * Semaphores used for critical sections in parallel mode:
- * RT_SEM_TREE ====> protects rtip->i->rti_solidheads[] lists, d_uses(solids)
- * RT_SEM_RESULTS => protects HeadRegion, mdl_min/max, d_uses(reg), nregions
+ * RT_SEM_TREE ====> protects rti_solidheads[], directory d_uses, prep publish
+ * RT_SEM_RESULTS => protects HeadRegion, mdl_min/max, nregions
  * RT_SEM_WORKER ==> (db_walk_dispatcher, from db_walk_tree)
  * RT_SEM_STATS ===> nsolids
  *

include/rt/soltab.h

@@ -60,6 +60,7 @@ struct soltab {
     struct rt_i *               st_rtip;        /**< @brief "up" pointer to rt_i */
     long                        st_uses;        /**< @brief Usage count, for instanced solids */
     int                         st_id;          /**< @brief Solid ident */
+    int                         st_prep_complete; /**< @brief prep publication flag */
     point_t                     st_center;      /**< @brief Centroid of solid */
     fastf_t                     st_aradius;     /**< @brief Radius of APPROXIMATING sphere */
     fastf_t                     st_bradius;     /**< @brief Radius of BOUNDING sphere */
@@ -77,7 +78,7 @@ struct soltab {
     struct bound_rpp *          st_piece_rpps;  /**< @brief bounding RPP of each piece of this solid */
 };
 
-#define RT_SOLTAB_INIT_ZERO {BU_LIST_INIT_ZERO, BU_LIST_INIT_ZERO, NULL, NULL, 0, 0, VINIT_ZERO, 0.0, 0.0, NULL, NULL, VINIT_ZERO, VINIT_ZERO, 0, BU_PTBL_INIT_ZERO, NULL, DB_FULL_PATH_INIT_ZERO, 0, 0, NULL}
+#define RT_SOLTAB_INIT_ZERO {BU_LIST_INIT_ZERO, BU_LIST_INIT_ZERO, NULL, NULL, 0, 0, 0, VINIT_ZERO, 0.0, 0.0, NULL, NULL, VINIT_ZERO, VINIT_ZERO, 0, BU_PTBL_INIT_ZERO, NULL, DB_FULL_PATH_INIT_ZERO, 0, 0, NULL}
 
 #define st_name         st_dp->d_namep
 #define RT_SOLTAB_NULL  ((struct soltab *)0)

src/gtools/CMakeLists.txt

@@ -1,5 +1,6 @@
 brlcad_adddata(gtransfer.c sample_applications)
 
+brlcad_addexec(brep-audit brep-audit.cpp "libbrep;libbv;librt;libbu")
 brlcad_addexec(ganalyze ganalyze.cpp "libanalyze;librt;libbu;libpkg;${LINENOISE_LIBRARIES}" NO_INSTALL)
 target_include_directories(ganalyze BEFORE PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../libtermio)
 
@@ -20,6 +21,7 @@ add_subdirectory(gsh)
 add_subdirectory(tests)
 
 cmakefiles(
+  brep-audit.cpp
   CMakeLists.txt
 )
 

src/gtools/tests/CMakeLists.txt

@@ -1,5 +1,34 @@
 brlcad_addexec(bigdb bigdb.c "libged;libwdb;librt;libbu" TEST_USESDATA)
 
+brlcad_add_test(
+  NAME brep_audit_surface_trees_all
+  COMMAND
+    brep-audit --surface-trees-all --surface-tree-depth=0
+    ${CMAKE_SOURCE_DIR}/regress/nurbs/sph.g
+)
+
+brlcad_add_test(
+  NAME brep_audit_prep_repeat
+  COMMAND
+    brep-audit --prep-only --prep-repeat=3
+    ${CMAKE_SOURCE_DIR}/regress/nurbs/sph.g sph.brep
+)
+
+brlcad_add_test(
+  NAME brep_audit_ray_depth_compare
+  COMMAND
+    brep-audit --ray-depth-compare --ray-count=4 --feature-ray-count=3
+    ${CMAKE_SOURCE_DIR}/regress/nurbs/sph.g sph.brep
+)
+
+brlcad_add_test(
+  NAME brep_audit_ray_depth_compare_multiface
+  COMMAND
+    brep-audit --ray-depth-compare --ray-count=4 --feature-ray-count=3
+    ${CMAKE_SOURCE_DIR}/src/libbrep/tests/nurbs_csg.g
+    recursive_child_linkage_complex_subtractions.s
+)
+
 brlcad_add_test(NAME slow-bigdb_1gb COMMAND bigdb 1)
 if(TARGET slow-bigdb_1gb)
   set_property(TEST slow-bigdb_1gb PROPERTY SKIP_RETURN_CODE 123)

src/libbrep/BRNode.cpp

@@ -21,6 +21,7 @@
 #include "common.h"
 
 #include <algorithm>
+#include <cfloat>
 
 #include "bu/log.h"
 #include "brep/brnode.h"
@@ -34,32 +35,27 @@ BRNode::BRNode(
     int trim_index,
     int adj_face_index,
     const ON_BoundingBox &node,
-    const ON_BrepFace *face,
     const ON_Interval &t,
-    bool innerTrim,
-    bool checkTrim,
-    bool trimmed) :
-    m_node(node),
-    m_v(),
+    bool innerTrim) :
     m_adj_face_index(adj_face_index),
     m_XIncreasing(false),
     m_Horizontal(false),
     m_Vertical(false),
     m_innerTrim(innerTrim),
-    m_stl(new Stl),
-    m_face(face),
-    m_u(),
     m_trim(curve),
     m_trim_index(trim_index),
     m_t(t),
-    m_checkTrim(checkTrim),
-    m_trimmed(trimmed),
-    m_estimate(),
-    m_slope(0.0),
-    m_bb_diag(0.0),
-    m_start(curve->PointAt(m_t[0])),
-    m_end(curve->PointAt(m_t[1]))
+    m_bbox_min{node.m_min.x, node.m_min.y},
+    m_bbox_max{node.m_max.x, node.m_max.y},
+    m_start{0.0, 0.0},
+    m_end{0.0, 0.0}
 {
+    const ON_3dPoint start = curve->PointAt(m_t[0]);
+    const ON_3dPoint end = curve->PointAt(m_t[1]);
+    m_start[0] = start.x;
+    m_start[1] = start.y;
+    m_end[0] = end.x;
+    m_end[1] = end.y;
     /* check for vertical segments they can be removed from trims
      * above (can't tell direction and don't need
      */
@@ -69,21 +65,6 @@ BRNode::BRNode(
     /*
      * should be okay since we split on Horz/Vert tangents
      */
-    if (m_end[X] < m_start[X]) {
-	m_u[0] = m_end[X];
-	m_u[1] = m_start[X];
-    } else {
-	m_u[0] = m_start[X];
-	m_u[1] = m_end[X];
-    }
-    if (m_end[Y] < m_start[Y]) {
-	m_v[0] = m_end[Y];
-	m_v[1] = m_start[Y];
-    } else {
-	m_v[0] = m_start[Y];
-	m_v[1] = m_end[Y];
-    }
-
     if (NEAR_EQUAL(m_end[X], m_start[X], BREP_UV_DIST_FUZZ)) {
 	m_Vertical = true;
 	if (m_innerTrim) {
@@ -98,299 +79,236 @@ BRNode::BRNode(
 	} else {
 	    m_XIncreasing = false;
 	}
-	m_slope = 0.0;
     } else {
 	if ((m_end[X] - m_start[X]) > 0.0) {
 	    m_XIncreasing = true;
 	} else {
 	    m_XIncreasing = false;
 	}
-	m_slope = (m_end[Y] - m_start[Y]) / (m_end[X] - m_start[X]);
     }
-    m_bb_diag = DIST_PNT_PNT(m_start, m_end);
-}
-
-
-BRNode::BRNode(const ON_BoundingBox &node) :
-    m_node(node),
-    m_v(),
-    m_adj_face_index(-1),
-    m_XIncreasing(false),
-    m_Horizontal(false),
-    m_Vertical(false),
-    m_innerTrim(false),
-    m_stl(new Stl),
-    m_face(NULL),
-    m_u(),
-    m_trim(NULL),
-    m_trim_index(-1),
-    m_t(),
-    m_checkTrim(true),
-    m_trimmed(false),
-    m_estimate(),
-    m_slope(0.0),
-    m_bb_diag(0.0),
-    m_start(ON_3dPoint::UnsetPoint),
-    m_end(ON_3dPoint::UnsetPoint)
-{
-    for (int i = 0; i < 3; i++) {
-	double d = m_node.m_max[i] - m_node.m_min[i];
-	if (NEAR_ZERO(d, ON_ZERO_TOLERANCE)) {
-	    m_node.m_min[i] -= 0.001;
-	    m_node.m_max[i] += 0.001;
-	}
-    }
-    m_start = m_node.m_min;
-    m_end = m_node.m_max;
-}
-
-
-BRNode::~BRNode()
-{
-    /* delete the children */
-    for (size_t i = 0; i < m_stl->m_children.size(); i++) {
-	delete m_stl->m_children[i];
-    }
-
-    delete m_stl;
 }
 
 
 BRNode::BRNode(Deserializer &deserializer, const ON_Brep &brep) :
-    m_node(),
-    m_v(),
     m_adj_face_index(-1),
     m_XIncreasing(false),
     m_Horizontal(false),
     m_Vertical(false),
     m_innerTrim(false),
-    m_stl(new Stl),
-    m_face(NULL),
-    m_u(),
     m_trim(NULL),
     m_trim_index(-1),
     m_t(),
-    m_checkTrim(false),
-    m_trimmed(false),
-    m_estimate(),
-    m_slope(0.0),
-    m_bb_diag(0.0),
-    m_start(ON_3dPoint::UnsetPoint),
-    m_end(ON_3dPoint::UnsetPoint)
+    m_bbox_min{0.0, 0.0},
+    m_bbox_max{0.0, 0.0},
+    m_start{0.0, 0.0},
+    m_end{0.0, 0.0}
 {
-    deserializer.read(m_node);
-    deserializer.read(m_u);
-    deserializer.read(m_v);
     deserializer.read(m_t);
-    deserializer.read(m_start);
-    deserializer.read(m_end);
-    deserializer.read(m_estimate);
-
-    m_slope = deserializer.read_double();
-    m_bb_diag = deserializer.read_double();
+    for (int direction = 0; direction < 2; ++direction) {
+	m_bbox_min[direction] = deserializer.read_double();
+	m_bbox_max[direction] = deserializer.read_double();
+	m_start[direction] = deserializer.read_double();
+	m_end[direction] = deserializer.read_double();
+    }
 
     const uint8_t bool_flags = deserializer.read_uint8();
-    const int face_index = deserializer.read_int32();
     m_trim_index = deserializer.read_int32();
-    m_adj_face_index = deserializer.read_uint32();
-    const std::size_t num_children = deserializer.read_uint32();
+    m_adj_face_index = deserializer.read_int32();
 
     m_XIncreasing = bool_flags & (1 << 0);
     m_Horizontal = bool_flags & (1 << 1);
     m_Vertical = bool_flags & (1 << 2);
     m_innerTrim = bool_flags & (1 << 3);
-    m_checkTrim = bool_flags & (1 << 4);
-    m_trimmed = bool_flags & (1 << 5);
-
-    if (face_index != -1) {
-	if (const ON_BrepFace * const face = brep.m_F.At(face_index))
-	    m_face = face;
-	else
-	    bu_bomb("invalid face index");
-    }
 
     if (m_trim_index != -1) {
 	if (const ON_BrepTrim * const trim = brep.m_T.At(m_trim_index))
 	    m_trim = trim->TrimCurveOf();
 	else
-	    bu_bomb("invalid face index");
+	    bu_bomb("invalid trim index");
     }
-
-    m_stl->m_children.resize(num_children);
-    for (std::vector<const BRNode *>::iterator it = m_stl->m_children.begin(); it != m_stl->m_children.end(); ++it)
-	*it = new BRNode(deserializer, brep);
 }
 
 
 void
 BRNode::serialize(Serializer &serializer) const
 {
-    const uint8_t bool_flags = (m_XIncreasing << 0) | (m_Horizontal << 1) | (m_Vertical << 2) | (m_innerTrim << 3) | (m_checkTrim << 4) | (m_trimmed << 5);
+    const uint8_t bool_flags = (m_XIncreasing << 0) |
+	(m_Horizontal << 1) | (m_Vertical << 2) |
+	(m_innerTrim << 3);
 
-    serializer.write(m_node);
-    serializer.write(m_u);
-    serializer.write(m_v);
     serializer.write(m_t);
-    serializer.write(m_start);
-    serializer.write(m_end);
-    serializer.write(m_estimate);
-
-    serializer.write_double(m_slope);
-    serializer.write_double(m_bb_diag);
+    for (int direction = 0; direction < 2; ++direction) {
+	serializer.write_double(m_bbox_min[direction]);
+	serializer.write_double(m_bbox_max[direction]);
+	serializer.write_double(m_start[direction]);
+	serializer.write_double(m_end[direction]);
+    }
 
     serializer.write_uint8(bool_flags);
-    serializer.write_int32(m_face ? m_face->m_face_index : -1);
     serializer.write_int32(m_trim_index);
-    serializer.write_uint32(m_adj_face_index);
-    serializer.write_uint32(m_stl->m_children.size());
-
-    for (std::vector<const BRNode *>::const_iterator it = m_stl->m_children.begin(); it != m_stl->m_children.end(); ++it)
-	(*it)->serialize(serializer);
+    serializer.write_int32(m_adj_face_index);
 }
 
 
-int
-BRNode::depth() const
+double
+BRNode::curveBBoxDistance(const ON_2dPoint &uv) const
 {
-    int d = 0;
-    for (size_t i = 0; i < m_stl->m_children.size(); i++) {
-	d = 1 + std::max(d, m_stl->m_children[i]->depth());
-    }
-    return d;
+    double du = 0.0;
+    double dv = 0.0;
+    if (uv.x < m_bbox_min[X])
+	du = m_bbox_min[X] - uv.x;
+    else if (uv.x > m_bbox_max[X])
+	du = uv.x - m_bbox_max[X];
+    if (uv.y < m_bbox_min[Y])
+	dv = m_bbox_min[Y] - uv.y;
+    else if (uv.y > m_bbox_max[Y])
+	dv = uv.y - m_bbox_max[Y];
+    return hypot(du, dv);
 }
 
 
-void
-BRNode::getLeaves(std::list<const BRNode *> &out_leaves) const
+double
+BRNode::curveDistance(const ON_2dPoint &uv) const
 {
-    if (!m_stl->m_children.empty()) {
-	for (size_t i = 0; i < m_stl->m_children.size(); i++) {
-	    m_stl->m_children[i]->getLeaves(out_leaves);
+    if (!m_trim)
+	return INFINITY;
+    const size_t sample_count = 8;
+    size_t best_sample = 0;
+    double best_squared = INFINITY;
+    for (size_t sample = 0; sample <= sample_count; ++sample) {
+	const double fraction = (double)sample / (double)sample_count;
+	const ON_3dPoint point = m_trim->PointAt(m_t.ParameterAt(fraction));
+	const double dx = point.x - uv.x;
+	const double dy = point.y - uv.y;
+	const double squared = dx * dx + dy * dy;
+	if (squared < best_squared) {
+	    best_squared = squared;
+	    best_sample = sample;
 	}
-    } else {
-	out_leaves.push_back(this);
     }
-}
 
-
-const BRNode *
-BRNode::closer(const ON_3dPoint &pt, const BRNode *left, const BRNode *right) const
-{
-    double ldist = pt.DistanceTo(left->m_estimate);
-    double rdist = pt.DistanceTo(right->m_estimate);
-    TRACE("\t" << ldist << " < " << rdist);
-    if (ldist < rdist) {
-	return left;
-    } else {
-	return right;
+    double low = best_sample ?
+	(double)(best_sample - 1) / (double)sample_count : 0.0;
+    double high = best_sample < sample_count ?
+	(double)(best_sample + 1) / (double)sample_count : 1.0;
+    const double golden = 0.6180339887498948482;
+    double first = high - golden * (high - low);
+    double second = low + golden * (high - low);
+    auto squared_distance = [this, &uv](double fraction) {
+	const ON_3dPoint point = m_trim->PointAt(m_t.ParameterAt(fraction));
+	const double dx = point.x - uv.x;
+	const double dy = point.y - uv.y;
+	return dx * dx + dy * dy;
+    };
+    double first_squared = squared_distance(first);
+    double second_squared = squared_distance(second);
+    for (int iteration = 0; iteration < 24; ++iteration) {
+	best_squared = std::min(best_squared,
+	    std::min(first_squared, second_squared));
+	if (first_squared < second_squared) {
+	    high = second;
+	    second = first;
+	    second_squared = first_squared;
+	    first = high - golden * (high - low);
+	    first_squared = squared_distance(first);
+	} else {
+	    low = first;
+	    first = second;
+	    first_squared = second_squared;
+	    second = low + golden * (high - low);
+	    second_squared = squared_distance(second);
+	}
     }
+    best_squared = std::min(best_squared,
+	std::min(first_squared, second_squared));
+    return sqrt(std::max(0.0, best_squared));
 }
 
 
 bool
-BRNode::isTrimmed(const ON_2dPoint &uv, double &trimdist) const
+BRNode::crossesAbove(const ON_2dPoint &uv, double &distance,
+	double endpoint_tolerance) const
 {
-    point_t bmin, bmax;
-    BRNode::GetBBox(bmin, bmax);
-    if ((bmin[X] <= uv[X]) && (uv[X] <= bmax[X]))   /* if check trim and in BBox */
-    {
-	fastf_t v = getCurveEstimateOfV(uv[X], 0.0000001);
-	trimdist = v - uv[Y];
-	if (uv[Y] <= v) {
-	    if (m_XIncreasing) {
-		return true;
-	    } else {
-		return false;
-	    }
-	} else if (uv[Y] > v) {
-	    if (!m_XIncreasing) {
-		return true;
-	    } else {
-		return false;
-	    }
-	} else {
-	    return true;
-	}
-    } else {
-	trimdist = -1.0;
-	if (m_trimmed) {
-	    return true;
-	} else {
+    distance = -1.0;
+    if (!m_trim || m_Vertical)
+	return false;
+    const double first_u = m_start[X];
+    const double second_u = m_end[X];
+    const bool exact =
+	((first_u <= uv.x && uv.x < second_u) ||
+	 (second_u <= uv.x && uv.x < first_u));
+    double query_u = uv.x;
+    if (!exact) {
+	if (!(endpoint_tolerance > 0.0))
+	    return false;
+	/* A closed trim owns both coincident curve-domain endpoints.  Its exact
+	 * half-open test already gives that vertex to one incident branch; the
+	 * imported-loop gap allowance must not admit the other branch too. */
+	if (m_trim->IsClosed())
+	    return false;
+	const double minimum_u = std::min(first_u, second_u);
+	const double minimum_t = first_u < second_u ? m_t[0] : m_t[1];
+	const ON_Interval curve_domain = m_trim->Domain();
+	const double parameter_scale = std::max(1.0,
+	    std::max(fabs(curve_domain.Min()), fabs(curve_domain.Max())));
+	const double parameter_fuzz = 64.0 * DBL_EPSILON * parameter_scale;
+	const bool minimum_is_trim_endpoint =
+	    fabs(minimum_t - curve_domain.Min()) <= parameter_fuzz ||
+	    fabs(minimum_t - curve_domain.Max()) <= parameter_fuzz;
+	const double endpoint_fuzz = std::max(BREP_UV_DIST_FUZZ,
+	    std::max(0.0, endpoint_tolerance));
+	if (uv.x < minimum_u &&
+		minimum_u - uv.x <= endpoint_fuzz &&
+		minimum_is_trim_endpoint)
+	    query_u = minimum_u;
+	else
 	    return false;
-	}
     }
+    const double curve_v = getCurveEstimateOfV(query_u, 1.0e-7);
+    distance = curve_v - uv.y;
+    return distance >= 0.0;
 }
 
 
-ON_2dPoint
-BRNode::getClosestPointEstimate(const ON_3dPoint &pt) const
+bool
+BRNode::crossesRight(const ON_2dPoint &uv, double &distance,
+	double endpoint_tolerance) const
 {
-    ON_Interval u, v;
-    return getClosestPointEstimate(pt, u, v);
-}
-
-
-ON_2dPoint
-BRNode::getClosestPointEstimate(const ON_3dPoint &pt, ON_Interval &u, ON_Interval &v) const
-{
-    if (isLeaf()) {
-	double uvs[5][2] = {{m_u.Min(), m_v.Min()},  /* include the corners for an easy refinement */
-			    {m_u.Max(), m_v.Min()},
-			    {m_u.Max(), m_v.Max()},
-			    {m_u.Min(), m_v.Max()},
-			    {m_u.Mid(), m_v.Mid()}
-	}; /* include the estimate */
-	ON_3dPoint corners[5];
-	const ON_Surface *surf = m_face->SurfaceOf();
-
-	u = m_u;
-	v = m_v;
-
-	/* ??? should we pass these in from SurfaceTree::curveBBox()
-	 * to avoid this recalculation?
-	 */
-	if (!surf->EvPoint(uvs[0][0], uvs[0][1], corners[0]) ||
-	    !surf->EvPoint(uvs[1][0], uvs[1][1], corners[1]) ||
-	    !surf->EvPoint(uvs[2][0], uvs[2][1], corners[2]) ||
-	    !surf->EvPoint(uvs[3][0], uvs[3][1], corners[3]))
-	{
-	    throw std::exception(); /* FIXME */
-	}
-	corners[4] = BRNode::m_estimate;
-
-	/* find the point on the curve closest to pt */
-	size_t mini = 0;
-	double mindist = pt.DistanceTo(corners[mini]);
-	double tmpdist;
-	for (size_t i = 1; i < 5; i++) {
-	    tmpdist = pt.DistanceTo(corners[i]);
-	    TRACE("\t" << mindist << " < " << tmpdist);
-	    if (tmpdist < mindist) {
-		mini = i;
-		mindist = tmpdist;
-	    }
-	}
-	TRACE("Closest: " << mindist << "; " << PT2(uvs[mini]));
-	return ON_2dPoint(uvs[mini][0], uvs[mini][1]);
-    } else {
-	if (!m_stl->m_children.empty()) {
-	    const BRNode *closestNode = m_stl->m_children[0];
-	    for (size_t i = 1; i < m_stl->m_children.size(); i++) {
-		closestNode = closer(pt, closestNode, m_stl->m_children[i]);
-	    }
-	    return closestNode->getClosestPointEstimate(pt, u, v);
-	} else {
-	    throw std::exception();
-	}
+    distance = -1.0;
+    if (!m_trim || m_Horizontal)
+	return false;
+    const double first_v = m_start[Y];
+    const double second_v = m_end[Y];
+    const bool exact =
+	((first_v <= uv.y && uv.y < second_v) ||
+	 (second_v <= uv.y && uv.y < first_v));
+    double query_v = uv.y;
+    if (!exact) {
+	if (!(endpoint_tolerance > 0.0))
+	    return false;
+	if (m_trim->IsClosed())
+	    return false;
+	const double minimum_v = std::min(first_v, second_v);
+	const double minimum_t = first_v < second_v ? m_t[0] : m_t[1];
+	const ON_Interval curve_domain = m_trim->Domain();
+	const double parameter_scale = std::max(1.0,
+	    std::max(fabs(curve_domain.Min()), fabs(curve_domain.Max())));
+	const double parameter_fuzz = 64.0 * DBL_EPSILON * parameter_scale;
+	const bool minimum_is_trim_endpoint =
+	    fabs(minimum_t - curve_domain.Min()) <= parameter_fuzz ||
+	    fabs(minimum_t - curve_domain.Max()) <= parameter_fuzz;
+	const double endpoint_fuzz = std::max(BREP_UV_DIST_FUZZ,
+	    std::max(0.0, endpoint_tolerance));
+	if (uv.y < minimum_v &&
+		minimum_v - uv.y <= endpoint_fuzz &&
+		minimum_is_trim_endpoint)
+	    query_v = minimum_v;
+	else
+	    return false;
     }
-}
-
-
-fastf_t
-BRNode::getLinearEstimateOfV(fastf_t u) const
-{
-    fastf_t v = m_start[Y] + m_slope * (u - m_start[X]);
-    return v;
+    const double curve_u = getCurveEstimateOfU(query_v, 1.0e-7);
+    distance = curve_u - uv.x;
+    return distance >= 0.0;
 }
 
 
@@ -401,13 +319,13 @@ BRNode::getCurveEstimateOfV(fastf_t u, fastf_t tol) const
     double Ta, Tb;
 
     if (m_start[X] < m_end[X]) {
-	VMOVE(A, m_start);
-	VMOVE(B, m_end);
+	VSET(A, m_start[X], m_start[Y], 0.0);
+	VSET(B, m_end[X], m_end[Y], 0.0);
 	Ta = m_t[0];
 	Tb = m_t[1];
     } else {
-	VMOVE(A, m_end);
-	VMOVE(B, m_start);
+	VSET(A, m_end[X], m_end[Y], 0.0);
+	VSET(B, m_start[X], m_start[Y], 0.0);
 	Ta = m_t[1];
 	Tb = m_t[0];
     }
@@ -493,13 +411,13 @@ BRNode::getCurveEstimateOfU(fastf_t v, fastf_t tol) const
     double Ta, Tb;
 
     if (m_start[Y] < m_end[Y]) {
-	VMOVE(A, m_start);
-	VMOVE(B, m_end);
+	VSET(A, m_start[X], m_start[Y], 0.0);
+	VSET(B, m_end[X], m_end[Y], 0.0);
 	Ta = m_t[0];
 	Tb = m_t[1];
     } else {
-	VMOVE(A, m_end);
-	VMOVE(B, m_start);
+	VSET(A, m_end[X], m_end[Y], 0.0);
+	VSET(B, m_start[X], m_start[Y], 0.0);
 	Ta = m_t[1];
 	Tb = m_t[0];
     }

src/libbrep/CMakeLists.txt

@@ -9,6 +9,7 @@ set(
   Subcurve.cpp
   Subsurface.cpp
   boolean.cpp
+  cobb.cpp
   cdt/cdt.cpp
   cdt/fast.cpp
   cdt/surf.cpp

src/libbrep/cdt/cdt.cpp

@@ -26,17 +26,426 @@
  */
 
 #include "common.h"
+#include <algorithm>
 #include <iostream>
+#include <map>
 #include <numeric>
 #include <queue>
 #include <string>
+#include <utility>
+#include <vector>
 #include "bg/chull.h"
+#include "bg/tri_pt.h"
 #include "bg/tri_tri.h"
 #include "./cdt.h"
 
 #define BREP_PLANAR_TOL 0.05
 #define MAX_TRIANGULATION_ATTEMPTS 5
 
+/* cpolyedge_t containers historically used the default pointer comparator.
+ * That is adequate for membership, but it must not determine the order of
+ * geometry-changing operations: allocator layout can then change the order in
+ * which singular trims are subdivided and, consequently, the final surface
+ * sampling.  Every singular segment has an authoritative STEP trim index,
+ * parameter interval, and pair of face-local topology vertices, which together
+ * provide a stable processing order. */
+static bool
+singular_edge_process_less(const cpolyedge_t *a, const cpolyedge_t *b)
+{
+    if (a == b)
+	return false;
+    if (!a)
+	return true;
+    if (!b)
+	return false;
+    if (a->trim_ind != b->trim_ind)
+	return a->trim_ind < b->trim_ind;
+    if (a->trim_start < b->trim_start)
+	return true;
+    if (b->trim_start < a->trim_start)
+	return false;
+    if (a->trim_end < b->trim_end)
+	return true;
+    if (b->trim_end < a->trim_end)
+	return false;
+    if (a->v2d[0] != b->v2d[0])
+	return a->v2d[0] < b->v2d[0];
+    return a->v2d[1] < b->v2d[1];
+}
+
+/* Correct triangle winding only after the complete mesh proves closed and
+ * manifold and misorientation is its sole defect.  bg_trimesh_sync establishes
+ * a consistent orientation for each connected component, but its seed may
+ * choose either of the two possible component orientations.  Preserve the
+ * OpenNURBS face orientation by selecting the result that changes the fewest
+ * source triangles in each component.  A tied component has no authoritative
+ * majority, so leave the mesh unchanged rather than guessing.
+ *
+ * This function is transactional: faces and face normals are not modified
+ * unless the candidate subsequently passes the complete solid test. */
+static int
+closed_mesh_orientation_sync(int *faces, int face_count,
+	fastf_t *vertices, int vertex_count, int *face_normals)
+{
+    if (!faces || !vertices || face_count < 4 || vertex_count < 4)
+	return -1;
+
+    struct bg_trimesh_solid_errors errors = BG_TRIMESH_SOLID_ERRORS_INIT_NULL;
+    const int initial_invalid = bg_trimesh_solid2(vertex_count, face_count,
+	vertices, faces, &errors);
+    if (!initial_invalid) {
+	bg_free_trimesh_solid_errors(&errors);
+	return 0;
+    }
+
+    const bool only_misoriented = errors.misoriented.count > 0 &&
+	errors.degenerate.count == 0 && errors.unmatched.count == 0 &&
+	errors.excess.count == 0;
+    bg_free_trimesh_solid_errors(&errors);
+    if (!only_misoriented)
+	return -1;
+
+    std::vector<int> original(faces, faces + 3 * face_count);
+    std::vector<int> candidate(original);
+    if (bg_trimesh_sync(candidate.data(), candidate.data(), face_count) <= 0)
+	return -1;
+
+    /* Determine the connected closed-shell components from shared edges. */
+    std::vector<std::vector<int>> adjacent((size_t)face_count);
+    std::map<std::pair<int, int>, int> first_edge_face;
+    for (int face = 0; face < face_count; ++face) {
+	for (int edge = 0; edge < 3; ++edge) {
+	    int vertex_a = original[(size_t)face * 3 + edge];
+	    int vertex_b = original[(size_t)face * 3 + (edge + 1) % 3];
+	    if (vertex_a > vertex_b)
+		std::swap(vertex_a, vertex_b);
+	    const std::pair<int, int> key(vertex_a, vertex_b);
+	    auto first = first_edge_face.find(key);
+	    if (first == first_edge_face.end()) {
+		first_edge_face[key] = face;
+		continue;
+	    }
+	    adjacent[(size_t)face].push_back(first->second);
+	    adjacent[(size_t)first->second].push_back(face);
+	}
+    }
+
+    std::vector<int> component_id((size_t)face_count, -1);
+    std::vector<std::vector<int>> components;
+    for (int seed = 0; seed < face_count; ++seed) {
+	if (component_id[(size_t)seed] >= 0)
+	    continue;
+	const int current_id = (int)components.size();
+	components.push_back(std::vector<int>());
+	std::queue<int> work;
+	work.push(seed);
+	component_id[(size_t)seed] = current_id;
+	while (!work.empty()) {
+	    const int face = work.front();
+	    work.pop();
+	    components[(size_t)current_id].push_back(face);
+	    for (int neighbor : adjacent[(size_t)face]) {
+		if (component_id[(size_t)neighbor] >= 0)
+		    continue;
+		component_id[(size_t)neighbor] = current_id;
+		work.push(neighbor);
+	    }
+	}
+    }
+
+    for (const std::vector<int> &component : components) {
+	size_t changed = 0;
+	for (int face : component) {
+	    const size_t offset = (size_t)face * 3;
+	    if (candidate[offset] != original[offset] ||
+		    candidate[offset + 1] != original[offset + 1] ||
+		    candidate[offset + 2] != original[offset + 2])
+		changed++;
+	}
+	if (changed * 2 == component.size())
+	    return -1;
+	if (changed * 2 > component.size()) {
+	    for (int face : component)
+		std::swap(candidate[(size_t)face * 3],
+		    candidate[(size_t)face * 3 + 1]);
+	}
+    }
+
+    if (bg_trimesh_solid2(vertex_count, face_count, vertices,
+	    candidate.data(), NULL) != 0)
+	return -1;
+
+    int changed_count = 0;
+    for (int face = 0; face < face_count; ++face) {
+	const size_t offset = (size_t)face * 3;
+	if (candidate[offset] == original[offset] &&
+		candidate[offset + 1] == original[offset + 1] &&
+		candidate[offset + 2] == original[offset + 2])
+	    continue;
+	changed_count++;
+	if (face_normals)
+	    std::swap(face_normals[offset], face_normals[offset + 1]);
+    }
+
+    std::copy(candidate.begin(), candidate.end(), faces);
+    return changed_count;
+}
+
+/* A failed local face remesh can leave a few open triangle islands in addition
+ * to the useful tessellation of the source solid.  Those islands are removed
+ * only under one of two proofs: every retained component is independently
+ * closed and solid, or (for a partial display mesh) every point of an island
+ * is already covered within tolerance by retained triangles from the same
+ * source BREP face.  Valid disjoint shells are retained, and no change is made
+ * if neither proof succeeds.
+ *
+ * The operation is transactional with respect to the caller's face arrays.
+ * Vertex arrays may retain now-unused entries; that is legal for a BoT and
+ * keeps the face-normal and point ownership maps stable. */
+static int
+closed_mesh_component_filter(int *faces, int *face_count,
+	fastf_t *vertices, int vertex_count, int *face_normals,
+	const int *source_faces = NULL, double overlap_tolerance = 0.0,
+	bool source_is_solid = true, bool emit_diagnostics = true)
+{
+    if (!faces || !face_count || !vertices || *face_count < 2 ||
+	    vertex_count < 4)
+	return 0;
+
+    /* Validate all indices before making the first transactional edit. */
+    for (int face = 0; face < *face_count; ++face) {
+	for (int corner = 0; corner < 3; ++corner) {
+	    const int vertex = faces[(size_t)face * 3 + corner];
+	    if (vertex < 0 || vertex >= vertex_count)
+		return 0;
+	}
+    }
+
+    /* Distinct face-mesh vertices may collapse to the same shared output
+     * point.  Cull the resulting zero-area index triangles before component
+     * discovery; otherwise their repeated edge can falsely bridge two
+     * components in this audit even though rt_bot_split correctly ignores
+     * that bridge. */
+    const int original_face_count = *face_count;
+    std::vector<int> compact_faces;
+    std::vector<int> compact_normals;
+    std::vector<int> compact_sources;
+    compact_faces.reserve((size_t)original_face_count * 3);
+    if (face_normals)
+	compact_normals.reserve((size_t)original_face_count * 3);
+    if (source_faces)
+	compact_sources.reserve((size_t)original_face_count);
+    for (int face = 0; face < original_face_count; ++face) {
+	const size_t offset = (size_t)face * 3;
+	if (faces[offset] == faces[offset + 1] ||
+		faces[offset + 1] == faces[offset + 2] ||
+		faces[offset + 2] == faces[offset])
+	    continue;
+	for (int corner = 0; corner < 3; ++corner) {
+	    compact_faces.push_back(faces[offset + corner]);
+	    if (face_normals)
+		compact_normals.push_back(face_normals[offset + corner]);
+	}
+	if (source_faces)
+	    compact_sources.push_back(source_faces[face]);
+    }
+    int removed_count = original_face_count -
+	(int)(compact_faces.size() / 3);
+    if (removed_count) {
+	std::copy(compact_faces.begin(), compact_faces.end(), faces);
+	if (face_normals)
+	    std::copy(compact_normals.begin(), compact_normals.end(),
+		face_normals);
+	*face_count = (int)(compact_faces.size() / 3);
+	if (source_faces)
+	    source_faces = compact_sources.data();
+    }
+    if (*face_count < 2)
+	return removed_count;
+
+    std::map<std::pair<int, int>, std::vector<int>> edge_faces;
+    for (int face = 0; face < *face_count; ++face) {
+	for (int corner = 0; corner < 3; ++corner) {
+	    int first = faces[(size_t)face * 3 + corner];
+	    int second = faces[(size_t)face * 3 + (corner + 1) % 3];
+	    if (first > second)
+		std::swap(first, second);
+	    edge_faces[std::make_pair(first, second)].push_back(face);
+	}
+    }
+
+    std::vector<int> component_id((size_t)*face_count, -1);
+    std::vector<std::vector<int>> components;
+    for (int seed = 0; seed < *face_count; ++seed) {
+	if (component_id[(size_t)seed] >= 0)
+	    continue;
+	const int current_id = (int)components.size();
+	components.push_back(std::vector<int>());
+	std::queue<int> work;
+	work.push(seed);
+	component_id[(size_t)seed] = current_id;
+	while (!work.empty()) {
+	    const int face = work.front();
+	    work.pop();
+	    components[(size_t)current_id].push_back(face);
+	    for (int corner = 0; corner < 3; ++corner) {
+		int first = faces[(size_t)face * 3 + corner];
+		int second =
+		    faces[(size_t)face * 3 + (corner + 1) % 3];
+		if (first > second)
+		    std::swap(first, second);
+		const std::vector<int> &edge_neighbors =
+		    edge_faces[std::make_pair(first, second)];
+		/* An edge used by more than two triangles is non-manifold.  It
+		 * must not join an otherwise closed shell to a local remesh
+		 * artifact during component discovery.  Restricting adjacency
+		 * to exactly two uses also matches the topological connectivity
+		 * required by a valid solid mesh. */
+		if (edge_neighbors.size() != 2)
+		    continue;
+		for (int neighbor : edge_neighbors) {
+		    if (component_id[(size_t)neighbor] >= 0)
+			continue;
+		    component_id[(size_t)neighbor] = current_id;
+		    work.push(neighbor);
+		}
+	    }
+	}
+    }
+    if (components.size() < 2)
+	return removed_count;
+
+    std::vector<bool> retain(components.size(), false);
+    size_t retained_face_count = 0;
+    for (size_t component = 0;
+	    source_is_solid && component < components.size(); ++component) {
+	std::vector<int> component_faces;
+	component_faces.reserve(components[component].size() * 3);
+	for (int face : components[component]) {
+	    const size_t offset = (size_t)face * 3;
+	    component_faces.push_back(faces[offset]);
+	    component_faces.push_back(faces[offset + 1]);
+	    component_faces.push_back(faces[offset + 2]);
+	}
+	if (component_faces.size() >= 12 &&
+		bg_trimesh_solid2(vertex_count,
+		    (int)components[component].size(), vertices,
+		    component_faces.data(), NULL) == 0) {
+	    retain[component] = true;
+	    retained_face_count += components[component].size();
+	}
+    }
+    if (retained_face_count == (size_t)*face_count)
+	return removed_count;
+
+    /* A partial display tessellation may not itself be closed (for example,
+     * another BREP face may have independently failed to tessellate).  We can
+     * still remove a disconnected local-remesh island when every one of its
+     * vertices and its centroid lies on retained triangles from the same
+     * source face.  This is an exact redundancy proof within the tessellation
+     * tolerance, not a largest-component guess. */
+    if (!retained_face_count && source_faces &&
+	    overlap_tolerance > 0.0) {
+	size_t largest = 0;
+	for (size_t component = 1; component < components.size(); ++component) {
+	    if (components[component].size() > components[largest].size())
+		largest = component;
+	}
+	double minimum_failed_distance = DBL_MAX;
+	int failed_source_face = -1;
+	const auto point_covered =
+	    [&](const point_t point, int source_face) {
+		double minimum_distance = DBL_MAX;
+		for (int candidate : components[largest]) {
+		    if (source_faces[candidate] != source_face)
+			continue;
+		    point_t triangle[3];
+		    for (int corner = 0; corner < 3; ++corner) {
+			const int vertex =
+			    faces[(size_t)candidate * 3 + corner];
+			VMOVE(triangle[corner],
+			    &vertices[(size_t)vertex * 3]);
+		    }
+		    const double distance = bg_tri_closest_pt(NULL, point,
+			triangle[0], triangle[1], triangle[2]);
+		    minimum_distance = std::min(minimum_distance, distance);
+		    if (distance <= overlap_tolerance)
+			return true;
+		}
+		minimum_failed_distance = minimum_distance;
+		failed_source_face = source_face;
+		return false;
+	    };
+	bool redundant = true;
+	for (size_t component = 0;
+		component < components.size() && redundant; ++component) {
+	    if (component == largest)
+		continue;
+	    for (int face : components[component]) {
+		point_t centroid = VINIT_ZERO;
+		for (int corner = 0; corner < 3; ++corner) {
+		    const int vertex = faces[(size_t)face * 3 + corner];
+		    point_t point;
+		    VMOVE(point, &vertices[(size_t)vertex * 3]);
+		    if (!point_covered(point, source_faces[face])) {
+			redundant = false;
+			break;
+		    }
+		    VADD2(centroid, centroid, point);
+		}
+		VSCALE(centroid, centroid, 1.0 / 3.0);
+		if (redundant &&
+			!point_covered(centroid, source_faces[face]))
+		    redundant = false;
+		if (!redundant)
+		    break;
+	    }
+	}
+	if (!redundant) {
+	    if (!emit_diagnostics)
+		return removed_count;
+	    bu_log("CDT retained an open disconnected component because source "
+		"face %d was not redundantly covered (distance %.17g, "
+		"tolerance %.17g)\n", failed_source_face,
+		minimum_failed_distance, overlap_tolerance);
+	    return removed_count;
+	}
+	retain[largest] = true;
+	retained_face_count = components[largest].size();
+    }
+    if (!retained_face_count)
+	return removed_count;
+
+    std::vector<int> candidate_faces;
+    std::vector<int> candidate_normals;
+    candidate_faces.reserve(retained_face_count * 3);
+    if (face_normals)
+	candidate_normals.reserve(retained_face_count * 3);
+    for (int face = 0; face < *face_count; ++face) {
+	if (!retain[(size_t)component_id[(size_t)face]])
+	    continue;
+	const size_t offset = (size_t)face * 3;
+	for (int corner = 0; corner < 3; ++corner) {
+	    candidate_faces.push_back(faces[offset + corner]);
+	    if (face_normals)
+		candidate_normals.push_back(face_normals[offset + corner]);
+	}
+    }
+    if (!source_faces &&
+	    bg_trimesh_solid2(vertex_count, (int)retained_face_count,
+		vertices, candidate_faces.data(), NULL) != 0)
+	return removed_count;
+
+    const int removed = removed_count +
+	*face_count - (int)retained_face_count;
+    std::copy(candidate_faces.begin(), candidate_faces.end(), faces);
+    if (face_normals)
+	std::copy(candidate_normals.begin(), candidate_normals.end(),
+	    face_normals);
+    *face_count = (int)retained_face_count;
+    return removed;
+}
+
 // TODO - get rid of all BN_TOL_DIST-only tolerances - if the object is
 // very small, that distance is too big (e.g. for linearity testing).
 
@@ -493,8 +902,12 @@ ON_Brep_CDT_Tessellate(struct ON_Brep_CDT_State *s_cdt, int face_cnt, int *faces
 	    int cnt = 0;
 	    wq = &w1;
 	    nq = &w2;
-	    nq->push(*(s_cdt->unsplit_singular_edges.begin()));
-	    s_cdt->unsplit_singular_edges.erase(s_cdt->unsplit_singular_edges.begin());
+	    std::set<cpolyedge_t *>::iterator first =
+		std::min_element(s_cdt->unsplit_singular_edges.begin(),
+			s_cdt->unsplit_singular_edges.end(),
+			singular_edge_process_less);
+	    nq->push(*first);
+	    s_cdt->unsplit_singular_edges.erase(first);
 	    while (cnt < 6) {
 		cnt = 0;
 		tmpq = wq;
@@ -504,8 +917,12 @@ ON_Brep_CDT_Tessellate(struct ON_Brep_CDT_State *s_cdt, int face_cnt, int *faces
 		    cpolyedge_t *ce = wq->front();
 		    wq->pop();
 		    std::set<cpolyedge_t *> nedges = split_singular_seg(s_cdt, ce, 0);
-		    std::set<cpolyedge_t *>::iterator n_it;
-		    for (n_it = nedges.begin(); n_it != nedges.end(); n_it++) {
+		    std::vector<cpolyedge_t *> ordered_edges(nedges.begin(),
+			    nedges.end());
+		    std::sort(ordered_edges.begin(), ordered_edges.end(),
+			    singular_edge_process_less);
+		    for (std::vector<cpolyedge_t *>::const_iterator n_it =
+			    ordered_edges.begin(); n_it != ordered_edges.end(); ++n_it) {
 			nq->push(*n_it);
 			cnt++;
 		    }
@@ -558,6 +975,8 @@ ON_Brep_CDT_Tessellate(struct ON_Brep_CDT_State *s_cdt, int face_cnt, int *faces
 	trimesh_error_report(s_cdt, valid_fcnt, valid_vcnt, valid_faces, valid_vertices, &se);
     }
 
+    bg_free_trimesh_solid_errors(&se);
+
     bu_free(valid_faces, "faces");
     bu_free(valid_vertices, "vertices");
 
@@ -610,7 +1029,7 @@ ON_Brep_CDT_Mesh(
     std::set<ON_3dPoint *> vfnormals;
     std::set<ON_3dPoint *> flip_normals;
     for (size_t fi = 0; fi < active_faces.size(); fi++) {
-	cdt_mesh_t *fmesh = &s_cdt->fmeshes[(int)fi];
+	cdt_mesh_t *fmesh = &s_cdt->fmeshes[active_faces[fi]];
 	RTree<size_t, double, 3>::Iterator tree_it;
 	fmesh->tris_tree.GetFirst(tree_it);
 	size_t t_ind;
@@ -699,8 +1118,10 @@ ON_Brep_CDT_Mesh(
     // Iterate over faces, adding points and faces to BoT container.  Note: all
     // 3D points should be geometrically unique in this final container.
     int face_cnt = 0;
+    std::vector<int> output_face_ids;
+    output_face_ids.reserve(triangle_cnt);
     for (size_t fi = 0; fi < active_faces.size(); fi++) {
-	cdt_mesh_t *fmesh = &s_cdt->fmeshes[(int)fi];
+	cdt_mesh_t *fmesh = &s_cdt->fmeshes[active_faces[fi]];
 	RTree<size_t, double, 3>::Iterator tree_it;
 	fmesh->tris_tree.GetFirst(tree_it);
 	size_t t_ind;
@@ -725,11 +1146,149 @@ ON_Brep_CDT_Mesh(
 		}
 	    }
 
+	    output_face_ids.push_back(fmesh->f_id);
 	    face_cnt++;
 	    ++tree_it;
 	}
     }
 
+    /* A complete, topologically solid BREP supplies authoritative shell
+     * orientation.  If tessellation produced an otherwise closed/manifold
+     * mesh with inconsistent triangle winding, repair that winding only when
+     * the fully synchronized candidate validates as a solid.  Partial/open
+     * meshes retain their original per-face ordering. */
+    if (!exp_face_cnt || !exp_faces) {
+	const bool source_is_solid = s_cdt->orig_brep &&
+	    s_cdt->orig_brep->IsSolid();
+	const int removed = closed_mesh_component_filter(*faces, fcnt,
+	    *vertices, *vcnt, face_normals ? *face_normals : NULL,
+	    output_face_ids.data(), std::max(s_cdt->absmin,
+		ON_ZERO_TOLERANCE), source_is_solid);
+	if (removed > 0) {
+	    bu_log("%s: removed %d redundant open triangle artifact%s after "
+		"transactional component validation\n",
+		s_cdt->name ? s_cdt->name : "BREP", removed,
+		removed == 1 ? "" : "s");
+	}
+	if (source_is_solid) {
+	    const int synchronized = closed_mesh_orientation_sync(*faces, *fcnt,
+		*vertices, *vcnt, face_normals ? *face_normals : NULL);
+	    if (synchronized > 0) {
+		bu_log("%s: synchronized %d triangle orientations after complete closed-mesh validation\n",
+		    s_cdt->name ? s_cdt->name : "BREP", synchronized);
+	    }
+	}
+    }
+
+    return 0;
+}
+
+int
+cdt_test_spurious_components(void)
+{
+    fastf_t vertices[] = {
+	0.0, 0.0, 0.0,
+	1.0, 0.0, 0.0,
+	0.0, 1.0, 0.0,
+	0.0, 0.0, 1.0,
+	2.0, 0.0, 0.0,
+	3.0, 0.0, 0.0,
+	2.0, 1.0, 0.0,
+	4.0, 0.0, 0.0
+    };
+    int tetra_and_open[] = {
+	0, 2, 1, 0, 1, 3, 0, 3, 2, 1, 2, 3,
+	4, 5, 6
+    };
+    int face_count = 5;
+    if (closed_mesh_component_filter(tetra_and_open, &face_count, vertices,
+	    8, NULL) != 1 || face_count != 4)
+	return 1;
+
+    int two_open_components[] = {0, 1, 2, 4, 5, 6};
+    face_count = 2;
+    if (closed_mesh_component_filter(two_open_components, &face_count,
+	    vertices, 8, NULL) != 0 || face_count != 2)
+	return 2;
+
+    fastf_t two_tetra_vertices[] = {
+	0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+	0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
+	3.0, 0.0, 0.0, 4.0, 0.0, 0.0,
+	3.0, 1.0, 0.0, 3.0, 0.0, 1.0
+    };
+    int two_tetra_faces[] = {
+	0, 2, 1, 0, 1, 3, 0, 3, 2, 1, 2, 3,
+	4, 6, 5, 4, 5, 7, 4, 7, 6, 5, 6, 7
+    };
+    face_count = 8;
+    if (closed_mesh_component_filter(two_tetra_faces, &face_count,
+	    two_tetra_vertices, 8, NULL) != 0 || face_count != 8)
+	return 3;
+
+    fastf_t open_duplicate_vertices[] = {
+	0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+	1.0, 1.0, 0.0, 0.0, 1.0, 0.0,
+	0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+	1.0, 1.0, 0.0
+    };
+    int open_duplicate_faces[] = {
+	0, 1, 2, 0, 2, 3, 4, 5, 6
+    };
+    const int open_duplicate_sources[] = {7, 7, 7};
+    face_count = 3;
+    if (closed_mesh_component_filter(open_duplicate_faces, &face_count,
+	    open_duplicate_vertices, 7, NULL, open_duplicate_sources,
+	    1.0e-12) != 1 || face_count != 2)
+	return 4;
+
+    open_duplicate_vertices[20] = 0.1;
+    int distinct_faces[] = {0, 1, 2, 0, 2, 3, 4, 5, 6};
+    face_count = 3;
+    if (closed_mesh_component_filter(distinct_faces, &face_count,
+	    open_duplicate_vertices, 7, NULL, open_duplicate_sources,
+	    1.0e-6, true, false) != 0 || face_count != 3)
+	return 5;
+
+    /* A redundant triangle attached to a closed shell through an edge with
+     * three uses is not part of that manifold shell.  This is the compact
+     * form of the NIST spherical-cap shading regression: treating the
+     * non-manifold edge as ordinary adjacency hid the extra component from
+     * the cleanup audit even though downstream BoT topology detected it. */
+    fastf_t nonmanifold_vertices[] = {
+	0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
+	0.0, 1.0, 0.0, 0.0, 0.0, 1.0,
+	1.0, 0.0, 0.0
+    };
+    int nonmanifold_faces[] = {
+	0, 2, 1, 0, 1, 3, 0, 3, 2, 1, 2, 3,
+	0, 2, 4
+    };
+    face_count = 5;
+    if (closed_mesh_component_filter(nonmanifold_faces, &face_count,
+	    nonmanifold_vertices, 5, NULL) != 1 || face_count != 4)
+	return 6;
+
+    /* Face-local vertices can collapse to one shared output vertex.  The
+     * resulting index-degenerate triangle must not survive into a BoT or
+     * distort component discovery. */
+    int degenerate_faces[] = {
+	0, 2, 1, 0, 1, 3, 0, 3, 2, 1, 2, 3,
+	4, 4, 1
+    };
+    face_count = 5;
+    if (closed_mesh_component_filter(degenerate_faces, &face_count,
+	    nonmanifold_vertices, 5, NULL) != 1 || face_count != 4)
+	return 7;
+
+    int invalid_faces[] = {0, 1, 2, 3, 3, 99};
+    const int invalid_original[] = {0, 1, 2, 3, 3, 99};
+    face_count = 2;
+    if (closed_mesh_component_filter(invalid_faces, &face_count,
+	    nonmanifold_vertices, 5, NULL) != 0 || face_count != 2 ||
+	    !std::equal(invalid_faces, invalid_faces + 6, invalid_original))
+	return 8;
+
     return 0;
 }
 

src/libbrep/cdt/cdt_util.cpp

@@ -861,4 +861,3 @@ int ON_Brep_CDT_VList(
 // c-file-style: "stroustrup"
 // End:
 // ex: shiftwidth=4 tabstop=8
-

src/libbrep/cdt/edge.cpp

@@ -1675,4 +1675,3 @@ finalize_rtrees(struct ON_Brep_CDT_State *s_cdt)
 // c-file-style: "stroustrup"
 // End:
 // ex: shiftwidth=4 tabstop=8
-

src/libbrep/cdt/fast.cpp

@@ -29,6 +29,7 @@
 #include "common.h"
 
 #include <vector>
+#include <atomic>
 #include <list>
 #include <map>
 #include <stack>
@@ -57,6 +58,8 @@
 #include "assert.h"
 
 #include "vmath.h"
+#include "bu/parallel.h"
+#include "bu/env.h"
 
 #include "bu/cv.h"
 #include "bu/opt.h"
@@ -290,6 +293,55 @@ getEdgePoints(ON_BrepTrim &trim,
 
     }
 
+    /* Boundary points shared by two faces must be taken from their common
+     * 3-D BREP edge, not independently from the two surface lifts.  Valid
+     * trims may differ from that edge by their asserted edge tolerance; using
+     * the lifts directly opens a visible crack in the fast shaded-display
+     * mesh.  A trim and its edge are not required to use identical parameter
+     * domains, so verify the normalized correspondence at every sample before
+     * changing any point.  If it is not valid for the entire trim, leave the
+     * trim on the existing general path rather than producing a mixture of
+     * surface and edge points. */
+    const ON_BrepEdge *edge = trim.Edge();
+    const ON_Curve *edge_curve = edge ? edge->EdgeCurveOf() : NULL;
+    const ON_Interval edge_domain = edge_curve ? edge_curve->Domain() :
+	ON_Interval::EmptyInterval;
+    const double edge_tolerance = edge ? std::max(tol->dist,
+	edge->m_tolerance) : tol->dist;
+
+    bool edge_parameterization_matches = edge_curve && range.IsIncreasing() &&
+	edge_domain.IsIncreasing();
+    if (edge_parameterization_matches) {
+	for (std::map<double, BrepTrimPoint *>::iterator sample =
+		param_points->begin(); sample != param_points->end(); ++sample) {
+	    BrepTrimPoint *point = sample->second;
+	    if (!point || !point->p3d) {
+		edge_parameterization_matches = false;
+		break;
+	    }
+	    double fraction = range.NormalizedParameterAt(point->t);
+	    if (trim.m_bRev3d)
+		fraction = 1.0 - fraction;
+	    const ON_3dPoint edge_point = edge_curve->PointAt(
+		edge_domain.ParameterAt(fraction));
+	    if (!edge_point.IsValid() ||
+		    edge_point.DistanceTo(*point->p3d) > edge_tolerance) {
+		edge_parameterization_matches = false;
+		break;
+	    }
+	}
+    }
+    if (edge_parameterization_matches) {
+	for (std::map<double, BrepTrimPoint *>::iterator sample =
+		param_points->begin(); sample != param_points->end(); ++sample) {
+	    BrepTrimPoint *point = sample->second;
+	    double fraction = range.NormalizedParameterAt(point->t);
+	    if (trim.m_bRev3d)
+		fraction = 1.0 - fraction;
+	    *point->p3d = edge_curve->PointAt(edge_domain.ParameterAt(fraction));
+	}
+    }
+
     return param_points;
 }
 
@@ -2298,6 +2350,55 @@ bg_CDT(std::vector<int> &faces, std::vector<fastf_t> &pnt_norms, std::vector<fas
 }
 
 
+struct fast_cdt_face_result {
+    std::vector<int> faces;
+    std::vector<fastf_t> norms;
+    std::vector<fastf_t> pnts;
+};
+
+struct fast_cdt_parallel_state {
+    const ON_Brep *brep;
+    const struct bg_tess_tol *ttol;
+    const struct bn_tol *tol;
+    std::vector<fast_cdt_face_result> *results;
+    std::atomic<int> next_face;
+    std::atomic<size_t> result_bytes;
+    std::atomic<bool> over_budget;
+    size_t result_budget;
+};
+
+static void
+fast_cdt_face_worker(int UNUSED(cpu), void *data)
+{
+    fast_cdt_parallel_state *state = (fast_cdt_parallel_state *)data;
+    const int face_cnt = state->brep->m_F.Count();
+
+    for (;;) {
+	if (state->over_budget.load())
+	    return;
+	const int face_index = state->next_face.fetch_add(1);
+	if (face_index >= face_cnt)
+	    return;
+
+	const ON_BrepFace& face = state->brep->m_F[face_index];
+	bg_CDT((*state->results)[face_index].faces,
+	       (*state->results)[face_index].norms,
+	       (*state->results)[face_index].pnts,
+	       face, state->ttol, state->tol);
+
+	const fast_cdt_face_result &result = (*state->results)[face_index];
+	const size_t result_bytes = result.faces.size() * sizeof(int) +
+	    (result.norms.size() + result.pnts.size()) * sizeof(fastf_t);
+	if (state->result_bytes.fetch_add(result_bytes) + result_bytes > state->result_budget) {
+	    state->over_budget = true;
+	    (*state->results)[face_index].faces.clear();
+	    (*state->results)[face_index].norms.clear();
+	    (*state->results)[face_index].pnts.clear();
+	    return;
+	}
+    }
+}
+
 int
 brep_cdt_fast(int **faces, int *face_cnt, vect_t **pnt_norms, point_t **pnts, int *pntcnt,
 	const ON_Brep *brep, int index, const struct bg_tess_tol *ttol, const struct bn_tol *tol)
@@ -2326,10 +2427,30 @@ brep_cdt_fast(int **faces, int *face_cnt, vect_t **pnt_norms, point_t **pnts, in
     std::vector<fastf_t> all_pnts;
 
     if (index == -1) {
-        for (index = 0; index < brep->m_F.Count(); index++) {
-            const ON_BrepFace& face = brep->m_F[index];
-            bg_CDT(all_faces, all_norms, all_pnts, face, ttol, tol);
-        }
+	/* Each face produces an independent set of indices and points.  Keep
+	 * those results separate while workers run, then concatenate them in
+	 * face order.  This avoids locking the hot CDT path and keeps output
+	 * stable across runs.  The BRep trim-user cleanup below remains serial. */
+	std::vector<fast_cdt_face_result> face_results(brep->m_F.Count());
+	ssize_t available = bu_mem(BU_MEM_AVAIL, NULL);
+	size_t result_budget = (available > 0) ? (size_t)available / 4 :
+	    (size_t)512 * 1024 * 1024;
+	fast_cdt_parallel_state state = {
+	    brep, ttol, tol, &face_results, 0, 0, false, result_budget
+	};
+	bu_parallel(fast_cdt_face_worker, 0, &state);
+	if (state.over_budget.load()) {
+	    bu_log("brep_cdt_fast: tessellation result exceeded memory budget (%zu bytes)\n", result_budget);
+	} else {
+	    for (size_t fi = 0; fi < face_results.size(); fi++) {
+		const size_t point_offset = all_pnts.size() / 3;
+		for (size_t i = 0; i < face_results[fi].faces.size(); i++) {
+		    all_faces.push_back((int)point_offset + face_results[fi].faces[i]);
+		}
+		all_norms.insert(all_norms.end(), face_results[fi].norms.begin(), face_results[fi].norms.end());
+		all_pnts.insert(all_pnts.end(), face_results[fi].pnts.begin(), face_results[fi].pnts.end());
+	    }
+	}
     } else if (index < brep->m_F.Count()) {
         const ON_BrepFaceArray& brep_faces = brep->m_F;
         if (index < brep_faces.Count()) {
@@ -2387,4 +2508,3 @@ brep_cdt_fast(int **faces, int *face_cnt, vect_t **pnt_norms, point_t **pnts, in
 // c-file-style: "stroustrup"
 // End:
 // ex: shiftwidth=4 tabstop=8
-

src/libbrep/cdt/mesh.cpp

@@ -490,6 +490,27 @@ cpolygon_t::add_ordered_edge(const struct edge2d_t &e)
     return nedge;
 }
 
+/* Return a stable boundary-loop starting edge.  poly is a pointer set, so
+ * dereferencing begin() makes traversal order depend on allocator layout.
+ * Geometry algorithms which walk next/prev from that arbitrary entry can then
+ * assign samples or exercise bounded refinements in a different order. */
+cpolyedge_t *
+cpolygon_t::first_edge() const
+{
+    cpolyedge_t *first = NULL;
+    for (std::set<cpolyedge_t *>::const_iterator edge = poly.begin();
+	    edge != poly.end(); ++edge) {
+	cpolyedge_t *candidate = *edge;
+	if (!candidate)
+	    continue;
+	if (!first || candidate->v2d[0] < first->v2d[0] ||
+		(candidate->v2d[0] == first->v2d[0] &&
+		 candidate->v2d[1] < first->v2d[1]))
+	    first = candidate;
+    }
+    return first;
+}
+
 void
 cpolygon_t::remove_ordered_edge(const struct edge2d_t &e)
 {
@@ -1150,7 +1171,17 @@ cpolygon_t::cdt(triangulation_t ttype)
     int *opoly = (int *)bu_calloc(poly.size()+1, sizeof(int), "polygon points");
 
     size_t vcnt = 1;
-    cpolyedge_t *pe = (*poly.begin());
+    /* bg_nested_poly_triangulate may select a different valid diagonal when
+     * the same closed boundary is cyclically rotated.  Downstream face-mesh
+     * stitching requires that choice to be independent of allocator layout,
+     * so use the same stable topological start as the LSCM projection. */
+    cpolyedge_t *pe = first_edge();
+    last_cdt_start_vertex = pe ? pe->v2d[0] : -1;
+    if (!pe) {
+	bu_free(opoly, "polygon points");
+	bu_free(bgp_2d, "free libbg 2d points array)");
+	return false;
+    }
     cpolyedge_t *first = pe;
     cpolyedge_t *next = pe->next;
 
@@ -3354,6 +3385,33 @@ loop_to_bgpoly(cpolygon_t *loop)
     return opoly;
 }
 
+static bool
+point_on_polygon_boundary(const point2d_t *points, int point_index,
+	const int *polygon, size_t polygon_point_count, double tolerance_sq)
+{
+    if (!points || point_index < 0 || !polygon || polygon_point_count < 2)
+	return false;
+    const double px = points[point_index][X];
+    const double py = points[point_index][Y];
+    for (size_t edge = 0; edge + 1 < polygon_point_count; ++edge) {
+	const int first = polygon[edge];
+	const int second = polygon[edge + 1];
+	const double ax = points[first][X];
+	const double ay = points[first][Y];
+	const double dx = points[second][X] - ax;
+	const double dy = points[second][Y] - ay;
+	const double length_sq = dx * dx + dy * dy;
+	double parameter = length_sq > DBL_EPSILON ?
+	    ((px - ax) * dx + (py - ay) * dy) / length_sq : 0.0;
+	parameter = std::max(0.0, std::min(1.0, parameter));
+	const double ex = px - (ax + parameter * dx);
+	const double ey = py - (ay + parameter * dy);
+	if (ex * ex + ey * ey <= tolerance_sq)
+	    return true;
+    }
+    return false;
+}
+
 bool
 cdt_mesh_t::cdt()
 {
@@ -3440,8 +3498,10 @@ cdt_mesh_t::cdt()
 	}
     }
 
-    // Build Steiner array, filtering out any points that fall inside a hole.
-    // Sampled interior points inside trimmed-away hole regions confuse detria.
+    // Build the Steiner array, filtering points in trimmed-away holes and
+    // points numerically on a constrained boundary.  Detria rejects a
+    // Steiner point on a constrained segment; the boundary already supplies
+    // the exact sample, so retaining that duplicate adds no geometry.
     // Pre-build per-hole 2D polygon arrays once (reused for each Steiner point test).
     std::vector<std::vector<double>> hole_polys_flat; // pairs of (x,y) stored flat
     std::vector<size_t> hole_polys_npts;
@@ -3460,8 +3520,24 @@ cdt_mesh_t::cdt()
 
     std::vector<int> steiner_vec;
     steiner_vec.reserve(m_interior_pnts.size());
+    const double scaled_u_range = (umax - umin) * uscale;
+    const double scaled_v_range = (vmax - vmin) * vscale;
+    const double boundary_scale = std::max(1.0,
+	hypot(scaled_u_range, scaled_v_range));
+    const double boundary_tolerance =
+	sqrt(DBL_EPSILON) * boundary_scale;
+    const double boundary_tolerance_sq =
+	boundary_tolerance * boundary_tolerance;
     for (auto p_it = m_interior_pnts.begin(); p_it != m_interior_pnts.end(); p_it++) {
 	int idx = (int)*p_it;
+	bool on_boundary = point_on_polygon_boundary(bgp_2d, idx, opoly,
+	    outer_loop.poly.size() + 1, boundary_tolerance_sq);
+	for (int hi = 0; hi < holes_cnt && !on_boundary; ++hi) {
+	    on_boundary = point_on_polygon_boundary(bgp_2d, idx,
+		holes_array[hi], holes_npts[hi], boundary_tolerance_sq);
+	}
+	if (on_boundary)
+	    continue;
 	bool in_hole = false;
 	for (int hi = 0; hi < holes_cnt && !in_hole; hi++) {
 	    point2d_t test_pnt;
@@ -3824,7 +3900,77 @@ cdt_mesh_t::repair()
 
     // Now that the out-and-out problem triangles have been handled,
     // remesh near singularities to try and produce more reasonable
-    // triangles.
+    // triangles.  This is a quality refinement, not a topology repair.  Keep
+    // it transactional: process_seed_tri may have replaced several patches
+    // before a later singular seed proves unmeshable.  Returning with those
+    // partial mutations used to leak isolated triangles into an otherwise
+    // valid face mesh (notably the NIST MBE PMI 6 spherical cap).
+    boundary_edges_stale = true;
+    boundary_edges_update();
+    const bool pre_singularity_mesh_valid = problem_edges.empty();
+    std::vector<triangle_t> pre_singularity_triangle_store;
+    std::vector<size_t> pre_singularity_active_triangles;
+    decltype(v2edges) pre_singularity_v2edges;
+    decltype(v2tris) pre_singularity_v2tris;
+    decltype(edges2tris) pre_singularity_edges2tris;
+    decltype(uedges2tris) pre_singularity_uedges2tris;
+    decltype(boundary_edges) pre_singularity_boundary_edges;
+    decltype(problem_edges) pre_singularity_problem_edges;
+    if (pre_singularity_mesh_valid && has_singularities) {
+	/* Preserve the complete indexed mesh state rather than rebuilding it
+	 * through tri_add.  Coincident singularity triangles may intentionally
+	 * have the same three 3-D vertex indices but different orientations;
+	 * tri_add's duplicate filter would discard one and make rollback itself
+	 * non-transactional. */
+	pre_singularity_triangle_store = tris_vect;
+	pre_singularity_v2edges = v2edges;
+	pre_singularity_v2tris = v2tris;
+	pre_singularity_edges2tris = edges2tris;
+	pre_singularity_uedges2tris = uedges2tris;
+	pre_singularity_boundary_edges = boundary_edges;
+	pre_singularity_problem_edges = problem_edges;
+	RTree<size_t, double, 3>::Iterator snapshot_it;
+	tris_tree.GetFirst(snapshot_it);
+	while (!snapshot_it.IsNull()) {
+	    pre_singularity_active_triangles.push_back(*snapshot_it);
+	    ++snapshot_it;
+	}
+    }
+    const auto restore_pre_singularity_mesh = [&]() {
+	if (!pre_singularity_mesh_valid ||
+		pre_singularity_triangle_store.empty() ||
+		pre_singularity_active_triangles.empty())
+	    return false;
+	tris_vect = pre_singularity_triangle_store;
+	tris_tree.RemoveAll();
+	for (size_t triangle_index : pre_singularity_active_triangles) {
+	    if (triangle_index >= tris_vect.size()) return false;
+	    triangle_t &triangle = tris_vect[triangle_index];
+	    triangle.m = this;
+	    ON_3dPoint *point = pnts[triangle.v[0]];
+	    ON_BoundingBox bounds(*point, *point);
+	    for (int vertex = 1; vertex < 3; ++vertex) {
+		point = pnts[triangle.v[vertex]];
+		bounds.Set(*point, true);
+	    }
+	    const double minimum[3] = {bounds.Min().x, bounds.Min().y,
+		bounds.Min().z};
+	    const double maximum[3] = {bounds.Max().x, bounds.Max().y,
+		bounds.Max().z};
+	    tris_tree.Insert(minimum, maximum, triangle_index);
+	}
+	v2edges = pre_singularity_v2edges;
+	v2tris = pre_singularity_v2tris;
+	edges2tris = pre_singularity_edges2tris;
+	uedges2tris = pre_singularity_uedges2tris;
+	boundary_edges = pre_singularity_boundary_edges;
+	problem_edges = pre_singularity_problem_edges;
+	seed_tris.clear();
+	new_tris.clear();
+	boundary_edges_stale = false;
+	bounding_box_stale = true;
+	return true;
+    };
 
     if (has_singularities) {
 	std::vector<triangle_t> s_tris = this->singularity_triangles();
@@ -3839,6 +3985,11 @@ cdt_mesh_t::repair()
 		bool pseed = process_seed_tri(seed, false, deg, NULL);
 
 		if (!pseed || seed_tris.size() >= st_size) {
+		    if (restore_pre_singularity_mesh()) {
+			bu_log("Face %d: retained the valid pre-refinement mesh after singularity quality remeshing made no progress\n",
+			    f_id);
+			return true;
+		    }
 		    std::cerr << f_id << ":  Error - failed to process refinement seed triangle!\n";
 		    struct bu_vls fname = BU_VLS_INIT_ZERO;
 		    bu_vls_sprintf(&fname, "%d-failed_seed.plot3", f_id);
@@ -3849,7 +4000,6 @@ cdt_mesh_t::repair()
 		    serialize(bu_vls_cstr(&fname));
 		    bu_vls_free(&fname);
 		    return false;
-		    break;
 		}
 
 		st_size = seed_tris.size();
@@ -4597,13 +4747,13 @@ cdt_mesh_t::serialize(const char *fname)
 	sfile << m_it->first << "," << m_it->second << "\n";
     }
 
-    sfile << "TRIANGLES_VECT" << tris_vect.size() << "\n";
+    sfile << "TRIANGLES_VECT " << tris_vect.size() << "\n";
     std::vector<triangle_t>::iterator t_it;
     for (t_it = tris_vect.begin(); t_it != tris_vect.end(); t_it++) {
 	sfile << (*t_it).v[0] << "," << (*t_it).v[1] << "," << (*t_it).v[2] << "," << (*t_it).ind << "\n";
     }
 
-    sfile << "TRIANGLES_TREE" << tris_tree.Count() << "\n";
+    sfile << "TRIANGLES_TREE " << tris_tree.Count() << "\n";
     RTree<size_t, double, 3>::Iterator tree_it;
     size_t t_ind;
     triangle_t tri;
@@ -4679,9 +4829,11 @@ cdt_mesh_t::deserialize(const char *fname)
     if (std::getline(sfile,switch_line)) {
 	if (switch_line == std::string("V1")) {
 	    version = 1;
+	} else if (switch_line == std::string("V2")) {
+	    version = 2;
 	}
     }
-    if (version < 1 || version > 1) {
+    if (version < 1 || version > 2) {
 	std::cerr << "Invalid deserialization file - format version " << switch_line << "\n";
 	return false;
     }
@@ -4716,11 +4868,27 @@ cdt_mesh_t::deserialize(const char *fname)
     problem_edges.clear();
 
     while (std::getline(sfile,switch_line)) {
-	std::cout << switch_line << "\n";
 	size_t spos = switch_line.find_first_of(' ');
-	std::string dtype = switch_line.substr(0, spos);
-	switch_line.erase(0, spos+1);
-	long lcnt = std::stol(switch_line);
+	std::string dtype;
+	std::string count_text;
+	if (spos != std::string::npos) {
+	    dtype = switch_line.substr(0, spos);
+	    count_text = switch_line.substr(spos + 1);
+	} else if (switch_line.compare(0, 14, "TRIANGLES_VECT") == 0) {
+	    /* V2 snapshots written before the delimiter fix concatenated these
+	     * two record names and their counts.  Accept those diagnostics so
+	     * existing failure captures remain replayable. */
+	    dtype = "TRIANGLES_VECT";
+	    count_text = switch_line.substr(14);
+	} else if (switch_line.compare(0, 14, "TRIANGLES_TREE") == 0) {
+	    dtype = "TRIANGLES_TREE";
+	    count_text = switch_line.substr(14);
+	} else {
+	    std::cerr << "Malformed serialization record: " << switch_line
+		<< "\nSerialization import failed.\n";
+	    return false;
+	}
+	long lcnt = std::stol(count_text);
 
 	if (dtype == std::string("POINTS")) {
 	    for (long i = 0; i < lcnt; i++) {
@@ -4781,7 +4949,8 @@ cdt_mesh_t::deserialize(const char *fname)
 	    continue;
 	}
 
-	if (dtype == std::string("TRIANGLES_VECT")) {
+	if (dtype == std::string("TRIANGLES") ||
+		dtype == std::string("TRIANGLES_VECT")) {
 	    for (long i = 0; i < lcnt; i++) {
 		std::string tline;
 		std::getline(sfile,tline);
@@ -4796,9 +4965,17 @@ cdt_mesh_t::deserialize(const char *fname)
 		long v2 = std::stol(v2str);
 		long v3 = std::stol(v3str);
 		triangle_t tri(v1, v2, v3);
-		// The tree is loaded separately - just do the basic population
-		tri.ind = tris_vect.size();
-		tris_vect.push_back(tri);
+		if (dtype == std::string("TRIANGLES")) {
+		    /* V1 stored only active triangles.  Rebuild the spatial index and
+		     * adjacency maps through the normal insertion path. */
+		    tri.m = this;
+		    tri_add(tri);
+		} else {
+		    /* V2 stores inactive vector entries as well; its following tree
+		     * record identifies and indexes the active subset. */
+		    tri.ind = tris_vect.size();
+		    tris_vect.push_back(tri);
+		}
 	    }
 	    continue;
 	}
@@ -5285,10 +5462,17 @@ cdt_mesh_t::lscm_reproject(cpolygon_t *polygon)
     // ── Step 1: Walk the boundary loop ───────────────────────────────────────
     // Collect exactly poly.size() unique boundary vertices (one per edge,
     // using each edge's start vertex) so LSCMParameterization receives N
-    // distinct vertices without a closing repeat.
+    // distinct vertices without a closing repeat.  polygon->poly is a set of
+    // pointers, so its begin() depends on heap layout.  The choice of starting
+    // edge rotates the convex LSCM boundary and can change the triangulation
+    // of difficult patches.  Choose the lexicographically first topological
+    // edge instead, making the result independent of allocation order and
+    // unrelated inputs such as the requested output object name.
     std::vector<int32_t> bnd_loop;
     {
-	cpolyedge_t *pe  = *polygon->poly.begin();
+	cpolyedge_t *pe = polygon->first_edge();
+	if (!pe)
+	    return false;
 	cpolyedge_t *cur = pe;
 	do {
 	    bnd_loop.push_back((int32_t)cur->v2d[0]);
@@ -6172,6 +6356,100 @@ void cdt_mesh_t::polygon_print_3d(cpolygon_t *polygon)
     std::cout << "\n";
 }
 
+/* Exercise the allocator-independent polygon start used by the actual CDT
+ * call.  A square is deliberately co-circular, so rotating its boundary input
+ * can select the opposite valid diagonal unless the triangulator receives a
+ * stable starting vertex.  Keep this compact guard here because cpolygon_t is
+ * an internal type whose symbols are not exported from libbrep. */
+int
+cdt_test_boundary_start(void)
+{
+    const auto initialize_regular_polygon = [](cpolygon_t &polygon,
+	    int point_count, int first_allocated_edge) {
+	for (int point = 0; point < point_count; ++point) {
+	    const double angle = 2.0 * ON_PI * point / point_count;
+	    ON_2dPoint p(cos(angle), sin(angle));
+	    polygon.add_point(p, point);
+	}
+	for (int offset = 0; offset < point_count; ++offset) {
+	    const int start = (first_allocated_edge + offset) % point_count;
+	    edge2d_t edge(start, (start + 1) % point_count);
+	    polygon.add_ordered_edge(edge);
+	}
+    };
+    const auto canonical_triangles = [](const cpolygon_t &polygon) {
+	std::set<std::vector<long> > result;
+	for (std::set<triangle_t>::const_iterator triangle =
+		polygon.tris.begin(); triangle != polygon.tris.end(); ++triangle) {
+	    std::vector<long> vertices;
+	    vertices.push_back(triangle->v[0]);
+	    vertices.push_back(triangle->v[1]);
+	    vertices.push_back(triangle->v[2]);
+	    std::sort(vertices.begin(), vertices.end());
+	    result.insert(vertices);
+	}
+	return result;
+    };
+    const auto release_edges = [](cpolygon_t &polygon) {
+	for (std::set<cpolyedge_t *>::iterator edge = polygon.poly.begin();
+		edge != polygon.poly.end(); ++edge)
+	    delete *edge;
+	polygon.poly.clear();
+    };
+
+    for (int point_count = 4; point_count <= 12; ++point_count) {
+	for (int first_allocated_edge = 1;
+		first_allocated_edge < point_count; ++first_allocated_edge) {
+	    cpolygon_t edge_zero_first;
+	    cpolygon_t rotated_first;
+	    initialize_regular_polygon(edge_zero_first, point_count, 0);
+	    initialize_regular_polygon(rotated_first, point_count,
+		first_allocated_edge);
+	    for (cpolyedge_t *edge : edge_zero_first.poly) {
+		if (!edge || edge->trim_ind != -1 || edge->loop_type != 0 ||
+			edge->defines_spnt || edge->split_status != 0 ||
+			edge->eseg != NULL) {
+		    release_edges(edge_zero_first);
+		    release_edges(rotated_first);
+		    return 1;
+		}
+	    }
+	    const bool valid = edge_zero_first.cdt() && rotated_first.cdt();
+	    const std::set<std::vector<long> > first =
+		canonical_triangles(edge_zero_first);
+	    const std::set<std::vector<long> > rotated =
+		canonical_triangles(rotated_first);
+	    release_edges(edge_zero_first);
+	    release_edges(rotated_first);
+	    if (!valid || edge_zero_first.last_cdt_start_vertex != 0 ||
+		    rotated_first.last_cdt_start_vertex != 0 ||
+		    first.size() != static_cast<size_t>(point_count - 2) ||
+		    first != rotated)
+		return 1;
+	}
+    }
+    return 0;
+}
+
+int
+cdt_test_boundary_steiner_filter(void)
+{
+    point2d_t points[6] = {
+	{0.0, 0.0}, {2.0, 0.0}, {2.0, 2.0}, {0.0, 2.0},
+	{1.0, 0.0}, {1.0, 1.0}
+    };
+    const int polygon[] = {0, 1, 2, 3, 0};
+    const double tolerance = sqrt(DBL_EPSILON) * sqrt(8.0);
+    if (!point_on_polygon_boundary(points, 4, polygon, 5,
+	    tolerance * tolerance))
+	return 1;
+    if (point_on_polygon_boundary(points, 5, polygon, 5,
+	    tolerance * tolerance))
+	return 1;
+    return 0;
+}
+
+
 // PImpl exposure of some mesh operations for use in tests
 struct cdt_bmesh_impl {
     cdt_mesh_t fmesh;
@@ -6199,8 +6477,7 @@ cdt_bmesh_deserialize(const char *fname, struct cdt_bmesh *m)
 {
     if (!fname || !m) return -1;
     if (!bu_file_exists(fname, NULL)) return -1;
-    m->i->fmesh.deserialize(fname);
-    return 0;
+    return m->i->fmesh.deserialize(fname) ? 0 : -1;
 }
 
 int

src/libbrep/cdt/mesh.h

@@ -441,7 +441,15 @@ class cpolyedge_t
 	    polygon = NULL;
 	    prev = NULL;
 	    next = NULL;
+	    trim_ind = -1;
+	    loop_type = 0;
 	    defines_spnt = false;
+	    spnt = ON_2dPoint::UnsetPoint;
+	    trim_start = 0.0;
+	    trim_end = 0.0;
+	    split_status = 0;
+	    v1_dist = 0.0;
+	    v2_dist = 0.0;
 	    eseg = NULL;
 	};
 
@@ -453,7 +461,15 @@ class cpolyedge_t
 	    polygon = NULL;
 	    prev = NULL;
 	    next = NULL;
+	    trim_ind = -1;
+	    loop_type = 0;
 	    defines_spnt = false;
+	    spnt = ON_2dPoint::UnsetPoint;
+	    trim_start = 0.0;
+	    trim_end = 0.0;
+	    split_status = 0;
+	    v1_dist = 0.0;
+	    v2_dist = 0.0;
 	    eseg = NULL;
 	};
 
@@ -506,6 +522,7 @@ class cpolygon_t
 	cpolyedge_t *add_edge(const struct uedge2d_t &e);
 	void remove_edge(const struct uedge2d_t &e);
 	std::set<cpolyedge_t *> replace_edges(std::set<uedge_t> &new_edges, std::set<uedge_t> &old_edges);
+	cpolyedge_t *first_edge() const;
 
 	/* Means to update the point array if we're incrementally building. orig_index should
 	 * identify the same point in the parent's index, so cdt() knows what triangles to
@@ -561,6 +578,7 @@ class cpolygon_t
 	std::set<triangle_t> visited_triangles;
 	std::set<uedge2d_t> active_edges;
 	std::set<uedge2d_t> self_isect_edges;
+	long last_cdt_start_vertex = -1;
 	ON_Plane tplane;
 	ON_Plane fit_plane;
 	ON_3dVector pdir;

src/libbrep/cdt/surf.cpp

@@ -1050,4 +1050,3 @@ GetInteriorPoints(struct ON_Brep_CDT_State *s_cdt, int face_index)
 // c-file-style: "stroustrup"
 // End:
 // ex: shiftwidth=4 tabstop=8
-

src/libbrep/cobb.cpp

@@ -0,0 +1,319 @@
+/*                       C O B B . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2013-2026 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * Based off of code from Ayam:
+ *
+ * This software is copyrighted (c) 1998-2012 by
+ * Randolf Schultz (randolf.schultz@gmail.com).
+ * All rights reserved.
+ *
+ * The author hereby grants permission to use, copy, modify, distribute,
+ * and license this software and its documentation for any purpose, provided
+ * that existing copyright notices are retained in all copies and that this
+ * notice is included verbatim in any distributions. No written agreement,
+ * license, or royalty fee is required for any of the authorized uses.
+ * Modifications to this software may be copyrighted by their authors
+ * and need not follow the licensing terms described here, provided that
+ * the new terms are clearly indicated on the first page of each file where
+ * they apply.
+ *
+ * IN NO EVENT SHALL THE AUTHOR OR DISTRIBUTORS BE LIABLE TO ANY PARTY
+ * FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES
+ * ARISING OUT OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY
+ * DERIVATIVES THEREOF, EVEN IF THE AUTHOR HAS BEEN ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ *
+ * THE AUTHOR AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES,
+ * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE
+ * IS PROVIDED ON AN "AS IS" BASIS, AND THE AUTHOR AND DISTRIBUTORS HAVE
+ * NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
+ * MODIFICATIONS.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ */
+/** @file cobb.cpp
+ *
+ * J. E. Cobb's six-patch rational Bezier sphere, with both the historical
+ * independent-face topology and a shared solid topology for seam tests.
+ */
+
+#include "common.h"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
+#include "brep/cobb.h"
+#include "opennurbs.h"
+
+
+ON_BezierSurface *
+ON_Brep_CobbSphereFace(double rotation_x, double rotation_z)
+{
+    const double t = sqrt(3.0);
+    const double d = sqrt(2.0);
+    const double s = sqrt(6.0);
+    const double cvs[5][5][4] = {
+	{
+	    {4.0*(1.0-t), 4.0*(1.0-t), 4.0*(1.0-t), 4.0*(3.0-t)},
+	    {-d, d*(t-4.0), d*(t-4.0), d*(3.0*t-2.0)},
+	    {0.0, 4.0*(1.0-2.0*t)/3.0, 4.0*(1.0-2.0*t)/3.0,
+		4.0*(5.0-t)/3.0},
+	    {d, d*(t-4.0), d*(t-4.0), d*(3.0*t-2.0)},
+	    {4.0*(t-1.0), 4.0*(1.0-t), 4.0*(1.0-t), 4.0*(3.0-t)}
+	},
+	{
+	    {d*(t-4.0), -d, d*(t-4.0), d*(3.0*t-2.0)},
+	    {(2.0-3.0*t)/2.0, (2.0-3.0*t)/2.0, -(t+6.0)/2.0,
+		(t+6.0)/2.0},
+	    {0.0, d*(2.0*t-7.0)/3.0, -5.0*s/3.0, d*(t+6.0)/3.0},
+	    {(3.0*t-2.0)/2.0, (2.0-3.0*t)/2.0, -(t+6.0)/2.0,
+		(t+6.0)/2.0},
+	    {d*(4.0-t), -d, d*(t-4.0), d*(3.0*t-2.0)}
+	},
+	{
+	    {4.0*(1.0-2.0*t)/3.0, 0.0, 4.0*(1.0-2.0*t)/3.0,
+		4.0*(5.0-t)/3.0},
+	    {d*(2.0*t-7.0)/3.0, 0.0, -5.0*s/3.0, d*(t+6.0)/3.0},
+	    {0.0, 0.0, 4.0*(t-5.0)/3.0, 4.0*(5.0*t-1.0)/9.0},
+	    {-d*(2.0*t-7.0)/3.0, 0.0, -5.0*s/3.0, d*(t+6.0)/3.0},
+	    {-4.0*(1.0-2.0*t)/3.0, 0.0, 4.0*(1.0-2.0*t)/3.0,
+		4.0*(5.0-t)/3.0}
+	},
+	{
+	    {d*(t-4.0), d, d*(t-4.0), d*(3.0*t-2.0)},
+	    {(2.0-3.0*t)/2.0, -(2.0-3.0*t)/2.0, -(t+6.0)/2.0,
+		(t+6.0)/2.0},
+	    {0.0, -d*(2.0*t-7.0)/3.0, -5.0*s/3.0, d*(t+6.0)/3.0},
+	    {(3.0*t-2.0)/2.0, -(2.0-3.0*t)/2.0, -(t+6.0)/2.0,
+		(t+6.0)/2.0},
+	    {d*(4.0-t), d, d*(t-4.0), d*(3.0*t-2.0)}
+	},
+	{
+	    {4.0*(1.0-t), -4.0*(1.0-t), 4.0*(1.0-t), 4.0*(3.0-t)},
+	    {-d, -d*(t-4.0), d*(t-4.0), d*(3.0*t-2.0)},
+	    {0.0, -4.0*(1.0-2.0*t)/3.0, 4.0*(1.0-2.0*t)/3.0,
+		4.0*(5.0-t)/3.0},
+	    {d, -d*(t-4.0), d*(t-4.0), d*(3.0*t-2.0)},
+	    {4.0*(t-1.0), -4.0*(1.0-t), 4.0*(1.0-t), 4.0*(3.0-t)}
+	}
+    };
+
+    ON_BezierSurface *surface = new ON_BezierSurface(3, true, 5, 5);
+    surface->ReserveCVCapacity(100);
+    for (int i = 0; i < 5; ++i) {
+	for (int j = 0; j < 5; ++j)
+	    surface->SetCV(i, j, ON_4dPoint(cvs[i][j][0], cvs[i][j][1],
+		cvs[i][j][2], cvs[i][j][3]));
+    }
+
+    const ON_3dPoint center(0.0, 0.0, 0.0);
+    surface->Rotate(rotation_x * ON_PI / 180.0,
+	ON_3dVector(1.0, 0.0, 0.0), center);
+    surface->Rotate(rotation_z * ON_PI / 180.0,
+	ON_3dVector(0.0, 0.0, 1.0), center);
+    return surface;
+}
+
+
+static ON_NurbsSurface *
+cobb_nurbs_face(double rotation_x, double rotation_z, double radius,
+	const ON_3dPoint &origin)
+{
+    ON_BezierSurface *bezier = ON_Brep_CobbSphereFace(rotation_x,
+	rotation_z);
+    if (!bezier)
+	return NULL;
+    ON_NurbsSurface *surface = ON_NurbsSurface::New();
+    if (!bezier->GetNurbForm(*surface)) {
+	delete surface;
+	delete bezier;
+	return NULL;
+    }
+    delete bezier;
+    surface->Scale(radius);
+    surface->Translate(ON_3dVector(origin.x, origin.y, origin.z));
+    return surface;
+}
+
+
+static bool
+cobb_face_outward(const ON_Surface *surface, const ON_3dPoint &origin)
+{
+    const ON_Interval u = surface->Domain(0);
+    const ON_Interval v = surface->Domain(1);
+    ON_3dPoint point;
+    ON_3dVector du;
+    ON_3dVector dv;
+    if (!surface->Ev1Der(u.Mid(), v.Mid(), point, du, dv))
+	return true;
+    return ON_CrossProduct(du, dv) * (point - origin) > 0.0;
+}
+
+
+ON_Brep *
+ON_Brep_CobbSphereUnsewn(double radius, const ON_3dPoint &origin)
+{
+    if (!(radius > 0.0) || !std::isfinite(radius) || !origin.IsValid())
+	return NULL;
+
+    const double rotations[6][2] = {
+	{0.0, 0.0}, {90.0, 0.0}, {180.0, 0.0},
+	{270.0, 0.0}, {90.0, 90.0}, {90.0, 270.0}
+    };
+    ON_Brep *brep = ON_Brep::New();
+    for (int i = 0; i < 6; ++i) {
+	ON_NurbsSurface *surface = cobb_nurbs_face(rotations[i][0],
+	    rotations[i][1], radius, origin);
+	if (!surface) {
+	    delete brep;
+	    return NULL;
+	}
+	ON_BrepFace *face = brep->NewFace(*surface);
+	if (!face) {
+	    delete surface;
+	    delete brep;
+	    return NULL;
+	}
+	face->m_bRev = !cobb_face_outward(surface, origin);
+	delete surface;
+    }
+    brep->Standardize();
+    brep->Compact();
+    return brep;
+}
+
+
+ON_Brep *
+ON_Brep_CobbSphereSewn(double radius, const ON_3dPoint &origin)
+{
+    if (!(radius > 0.0) || !std::isfinite(radius) || !origin.IsValid())
+	return NULL;
+
+    const double rotations[6][2] = {
+	{0.0, 0.0}, {90.0, 0.0}, {180.0, 0.0},
+	{270.0, 0.0}, {90.0, 90.0}, {90.0, 270.0}
+    };
+    const double coordinate_scale = std::max(std::max(radius, 1.0),
+	std::max(std::max(fabs(origin.x), fabs(origin.y)), fabs(origin.z)));
+    const double corner_tolerance = std::max(ON_ZERO_TOLERANCE,
+	128.0 * std::numeric_limits<double>::epsilon() * coordinate_scale);
+    ON_Brep *brep = ON_Brep::New();
+    for (int face_index = 0; face_index < 6; ++face_index) {
+	ON_NurbsSurface *surface = cobb_nurbs_face(rotations[face_index][0],
+	    rotations[face_index][1], radius, origin);
+	if (!surface) {
+	    delete brep;
+	    return NULL;
+	}
+
+	const ON_Interval u = surface->Domain(0);
+	const ON_Interval v = surface->Domain(1);
+	const ON_3dPoint corners[4] = {
+	    surface->PointAt(u.Min(), v.Min()),
+	    surface->PointAt(u.Max(), v.Min()),
+	    surface->PointAt(u.Max(), v.Max()),
+	    surface->PointAt(u.Min(), v.Max())
+	};
+	int vid[4] = {-1, -1, -1, -1};
+	for (int corner = 0; corner < 4; ++corner) {
+	    for (int vertex = 0; vertex < brep->m_V.Count(); ++vertex) {
+		if (corners[corner].DistanceTo(brep->m_V[vertex].Point()) <=
+			corner_tolerance) {
+		    vid[corner] = vertex;
+		    break;
+		}
+	    }
+	}
+
+	/* Oriented outer-loop corner pairs for south, east, north, west. */
+	const int side_corners[4][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}};
+	int eid[4] = {-1, -1, -1, -1};
+	bool reversed[4] = {false, false, false, false};
+	for (int side = 0; side < 4; ++side) {
+	    const int from = vid[side_corners[side][0]];
+	    const int to = vid[side_corners[side][1]];
+	    if (from < 0 || to < 0)
+		continue;
+	    for (int edge = 0; edge < brep->m_E.Count(); ++edge) {
+		const int edge_from = brep->m_E[edge].m_vi[0];
+		const int edge_to = brep->m_E[edge].m_vi[1];
+		if ((edge_from == from && edge_to == to) ||
+			(edge_from == to && edge_to == from)) {
+		    eid[side] = edge;
+		    reversed[side] = edge_from != from;
+		    break;
+		}
+	    }
+	}
+
+	const bool outward = cobb_face_outward(surface, origin);
+	ON_BrepFace *face = brep->NewFace(surface, vid, eid, reversed);
+	if (!face) {
+	    delete surface;
+	    delete brep;
+	    return NULL;
+	}
+	face->m_bRev = !outward;
+	for (int corner = 0; corner < 4; ++corner)
+	    brep->m_V[vid[corner]].m_tolerance = 0.0;
+	for (int side = 0; side < 4; ++side)
+	    brep->m_E[eid[side]].m_tolerance = 0.0;
+    }
+
+    for (int trim = 0; trim < brep->m_T.Count(); ++trim) {
+	brep->m_T[trim].m_tolerance[0] = 0.0;
+	brep->m_T[trim].m_tolerance[1] = 0.0;
+    }
+    brep->Standardize();
+    brep->Compact();
+    return brep;
+}
+
+
+double
+ON_Brep_CobbSphereMaxRadialError(const ON_Brep *brep, double radius,
+	const ON_3dPoint &origin, int samples_per_direction)
+{
+    if (!brep || !(radius > 0.0) || !std::isfinite(radius) ||
+	    !origin.IsValid() || samples_per_direction < 2)
+	return std::numeric_limits<double>::infinity();
+
+    double maximum_error = 0.0;
+    for (int si = 0; si < brep->m_S.Count(); ++si) {
+	const ON_Surface *surface = brep->m_S[si];
+	if (!surface)
+	    return std::numeric_limits<double>::infinity();
+	const ON_Interval u = surface->Domain(0);
+	const ON_Interval v = surface->Domain(1);
+	for (int i = 0; i < samples_per_direction; ++i) {
+	    const double up = u.ParameterAt((double)i /
+		(double)(samples_per_direction - 1));
+	    for (int j = 0; j < samples_per_direction; ++j) {
+		const double vp = v.ParameterAt((double)j /
+		    (double)(samples_per_direction - 1));
+		const double sample_radius = surface->PointAt(up, vp).
+		    DistanceTo(origin);
+		maximum_error = std::max(maximum_error,
+		    fabs(sample_radius - radius));
+	    }
+	}
+    }
+    return maximum_error;
+}
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C++
+ * c-basic-offset: 4
+ * indent-tabs-mode: t
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */

src/libbrep/tests/CMakeLists.txt

@@ -5,8 +5,29 @@ find_package(Threads REQUIRED)
 brlcad_addexec(test_brep_ppx ppx.cpp "libbrep" TEST)
 brlcad_add_test(NAME brep_ppx COMMAND test_brep_ppx)
 
-brlcad_addexec(test_brep_assembly assembly.cpp "libbrep" TEST)
-brlcad_add_test(NAME brep_assembly COMMAND test_brep_assembly)
+# Replay serialized face-mesh failures without rerunning a complete BRep
+# tessellation.  mesh.cpp emits this format at the exact repair failure site,
+# making it useful both for diagnosis and for compact regression fixtures.
+brlcad_addexec(test_brep_cdt_mesh brep_cdt_mesh.cpp "libbrep" TEST)
+brlcad_add_test(
+  NAME brep_cdt_mesh_single_bad_edge
+  COMMAND test_brep_cdt_mesh "${CMAKE_CURRENT_SOURCE_DIR}/single_bad_edge.cdtmesh"
+)
+
+# A cyclic rotation of the same co-circular boundary must not change the CDT
+# diagonal merely because its cpolyedge allocations have different addresses.
+brlcad_addexec(test_brep_cdt_boundary_start brep_cdt_boundary_start.cpp "libbrep" TEST)
+brlcad_add_test(NAME brep_cdt_boundary_start COMMAND test_brep_cdt_boundary_start)
+
+# An interior sample exactly on a constrained face boundary must not be passed
+# to the triangulator a second time as a Steiner point.
+brlcad_addexec(test_brep_cdt_boundary_steiner brep_cdt_boundary_steiner.cpp "libbrep" TEST)
+brlcad_add_test(NAME brep_cdt_boundary_steiner COMMAND test_brep_cdt_boundary_steiner)
+
+# A closed source solid may not retain redundant open local-remesh islands
+# after either retained-solid or same-source-face coverage validation.
+brlcad_addexec(test_brep_cdt_spurious_components brep_cdt_spurious_components.cpp "libbrep" TEST)
+brlcad_add_test(NAME brep_cdt_spurious_components COMMAND test_brep_cdt_spurious_components)
 
 # Standalone CDT hole-triangulation tests: verify that bg_nested_poly_triangulate
 # produces no triangles whose centroid falls inside a hole polygon, for two
@@ -26,13 +47,20 @@ endif()
 brlcad_addexec(test_pullback_context pullback_context.cpp "libbrep;Threads::Threads" TEST)
 brlcad_add_test(NAME brep_pullback_context COMMAND test_pullback_context)
 
+# The adaptive SurfaceTree depth budget starts only after closed seams and all
+# internal NURBS knot boundaries have established well-behaved subdomains.
+brlcad_addexec(test_surfacetree_knots surfacetree_knots.cpp "libbrep" TEST)
+brlcad_add_test(NAME brep_surfacetree_knots COMMAND test_surfacetree_knots)
+
 cmakefiles(
   CMakeLists.txt
   ayam_hyperbolid.3dm
-  assembly.cpp
   brep_cdt_invalid.cpp
   brep_cdt_hole_face17.c
   brep_cdt_hole_face27.c
+  brep_cdt_boundary_start.cpp
+  brep_cdt_boundary_steiner.cpp
+  brep_cdt_spurious_components.cpp
   brep_cdt_mesh.cpp
   cyl_hole.g
   cyls_holes.g
@@ -43,6 +71,7 @@ cmakefiles(
   plate_mode_samples_rhino.stp
   pullback_context.cpp
   single_bad_edge.cdtmesh
+  surfacetree_knots.cpp
   test_curve_intersect.cpp
   test_point_intersect.cpp
   vsp.stp

src/libbrep/tests/brep_cdt_boundary_start.cpp

@@ -0,0 +1,25 @@
+/*                    B R E P _ C D T _ B O U N D A R Y _ S T A R T . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2026 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ */
+
+#include "common.h"
+
+#include "bu/app.h"
+#include "brep/cdt.h"
+
+int
+main(int argc, const char **argv)
+{
+    bu_setprogname(argv[0]);
+    if (argc != 1)
+	return 1;
+
+    return cdt_test_boundary_start();
+}

src/libbrep/tests/brep_cdt_boundary_steiner.cpp

@@ -0,0 +1,19 @@
+/*                   B R E P _ C D T _ B O U N D A R Y _ S T E I N E R . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2026 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ */
+
+#include "common.h"
+#include "brep/cdt.h"
+
+int
+main()
+{
+    return cdt_test_boundary_steiner_filter();
+}

src/libbrep/tests/brep_cdt_mesh.cpp

@@ -31,6 +31,7 @@ main(int argc, const char **argv)
 
     if (argc != 2) {
 	std::cerr << "brep_cdt_mesh <serialization_file>\n";
+	return 1;
     }
     struct cdt_bmesh *fmesh;
     if (cdt_bmesh_create(&fmesh)) return -1;
@@ -57,4 +58,3 @@ main(int argc, const char **argv)
 // c-file-style: "stroustrup"
 // End:
 // ex: shiftwidth=4 tabstop=8
-

src/libbrep/tests/brep_cdt_spurious_components.cpp

@@ -0,0 +1,19 @@
+/*                     B R E P _ C D T _ S P U R I O U S _ C O M P O N E N T S . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2026 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ */
+
+#include "common.h"
+#include "brep/cdt.h"
+
+int
+main()
+{
+    return cdt_test_spurious_components();
+}

src/libbrep/tests/surfacetree_knots.cpp

@@ -0,0 +1,121 @@
+/*                 S U R F A C E T R E E _ K N O T S . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2026 United States Government as represented by
+ * the U.S. Army Research Laboratory.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
+ */
+
+#include "common.h"
+
+#include <cmath>
+#include <cstdio>
+#include <list>
+#include <vector>
+
+#include "brep.h"
+#include "brep/surfacetree.h"
+
+
+int
+main()
+{
+    ON_NurbsSurface surface(3, false, 3, 3, 8, 7);
+    if (!surface.MakeClampedUniformKnotVector(0) ||
+	    !surface.MakeClampedUniformKnotVector(1)) {
+	std::printf("FAIL: knot vector construction\n");
+	return 1;
+    }
+    const ON_Interval domain[2] = {surface.Domain(0), surface.Domain(1)};
+    for (int u = 0; u < surface.CVCount(0); ++u) {
+	for (int v = 0; v < surface.CVCount(1); ++v) {
+	    const double uf = (double)u / (surface.CVCount(0) - 1);
+	    const double vf = (double)v / (surface.CVCount(1) - 1);
+	    const double x = domain[0].ParameterAt(uf);
+	    const double y = domain[1].ParameterAt(vf);
+	    const double z = 0.05 * std::sin(0.7 * u) * std::cos(0.6 * v);
+	    if (!surface.SetCV(u, v, ON_3dPoint(x, y, z))) {
+		std::printf("FAIL: control point construction\n");
+		return 1;
+	    }
+	}
+    }
+    if (!surface.IsValid()) {
+	std::printf("FAIL: invalid test surface\n");
+	return 1;
+    }
+
+    ON_Brep brep;
+    ON_BrepFace *face = brep.NewFace(surface);
+    if (!face) {
+	std::printf("FAIL: face construction\n");
+	return 1;
+    }
+    brlcad::SurfaceTree tree(face, false, 0);
+    if (!tree.Valid()) {
+	std::printf("FAIL: depth-zero SurfaceTree construction reason=%d\n",
+	    (int)tree.Failure());
+	return 1;
+    }
+
+    const ON_Surface *tree_surface = face->SurfaceOf();
+    std::vector<double> spans[2];
+    size_t expected_leaves = 1;
+    for (int direction = 0; direction < 2; ++direction) {
+	const int span_count = tree_surface->SpanCount(direction);
+	if (span_count < 2) {
+	    std::printf("FAIL: direction %d has only %d spans\n", direction,
+		span_count);
+	    return 1;
+	}
+	spans[direction].resize((size_t)span_count + 1);
+	if (!tree_surface->GetSpanVector(direction,
+		spans[direction].data())) {
+	    std::printf("FAIL: direction %d span extraction\n", direction);
+	    return 1;
+	}
+	expected_leaves *= (size_t)span_count;
+    }
+
+    std::list<const brlcad::BBNode *> leaves;
+    tree.getLeaves(leaves);
+    if (leaves.size() != expected_leaves) {
+	std::printf("FAIL: depth-zero leaves=%zu expected=%zu\n", leaves.size(),
+	    expected_leaves);
+	return 1;
+    }
+    for (std::list<const brlcad::BBNode *>::const_iterator leaf = leaves.begin();
+	    leaf != leaves.end(); ++leaf) {
+	const ON_Interval interval[2] = {(*leaf)->m_u, (*leaf)->m_v};
+	for (int direction = 0; direction < 2; ++direction) {
+	    for (size_t split = 1;
+		    split + 1 < spans[direction].size(); ++split) {
+		const double knot = spans[direction][split];
+		if (knot > interval[direction].Min() &&
+			knot < interval[direction].Max()) {
+		    std::printf("FAIL: leaf crosses direction %d knot %.17g "
+			"in [%.17g %.17g]\n", direction, knot,
+			interval[direction].Min(), interval[direction].Max());
+		    return 1;
+		}
+	    }
+	}
+    }
+    std::printf("SurfaceTree structural knots: spans=%d/%d leaves=%zu "
+	"adaptive-depth=0\n", tree_surface->SpanCount(0),
+	tree_surface->SpanCount(1), leaves.size());
+    return 0;
+}
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C++
+ * c-basic-offset: 4
+ * indent-tabs-mode: t
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */

src/libbu/dir.c

@@ -578,7 +578,7 @@ static const char *
 vdir(char *result, size_t len, va_list args)
 {
     struct bu_vls vls = BU_VLS_INIT_ZERO;
-    static char buf[MAXPATHLEN] = {0};
+    static THREADLOCAL char buf[MAXPATHLEN] = {0};
     const char *cpath;
     uintptr_t arg;
 

src/libbu/heap.c

@@ -214,13 +214,19 @@ bu_heap_get(size_t sz)
 
     /* init */
     if (heap->count == 0) {
-
-	if (registered++ == 0) {
+	/*
+	 * Different per-CPU bins can be initialized concurrently.  Serialize
+	 * the process-global environment check and atexit registration rather
+	 * than racing on the lazy initialization flag.
+	 */
+	bu_semaphore_acquire(BU_SEM_SYSCALL);
+	if (!registered) {
+	    registered = 1;
 	    ret = getenv("BU_HEAP_PRINT");
-	    if ((++registered == 2) && (ret && atoi(ret) > 0)) {
+	    if (ret && atoi(ret) > 0)
 		atexit(heap_print);
-	    }
 	}
+	bu_semaphore_release(BU_SEM_SYSCALL);
 
 	heap->count++;
 	heap->pages = (char **)bu_malloc(1 * sizeof(char *), "heap malloc pages[]");
@@ -252,7 +258,14 @@ bu_heap_put(void *ptr, size_t sz)
 	return;
     }
 
-    /* TODO: actually do something useful :) */
+	/* TODO: This bump allocator cannot safely reuse individual objects yet.
+	 * Allocations are selected by bu_parallel_id(), but put has no owner or
+	 * page metadata and objects may be destroyed by a different thread.
+	 * Any implementation must define cross-thread return, page liveness and
+	 * reclamation/compaction semantics before adding per-bin free lists.
+	 * Bulk-lifetime clients should use an owning arena or compact container
+	 * instead of relying on this process-lifetime storage.
+	 */
 
     return;
 }

src/libged/brep/plot.cpp

@@ -678,18 +678,17 @@ drawBBNode(const SurfaceTree* st, struct bu_list *vlfree, struct bv_vlblock *vbp
 	    return;
 	}
     } else {
-	if (!node->get_children().empty()) {
-	    for (std::vector<BBNode*>::const_iterator childnode = node->get_children().begin(); childnode
-		 != node->get_children().end(); childnode++) {
-		drawBBNode(st, vlfree, vbp, *childnode);
-	    }
-	}
+	for (std::size_t child_index = 0; child_index < node->childCount();
+		++child_index)
+	    drawBBNode(st, vlfree, vbp, node->child(child_index));
     }
 }
 
 static void
 plotFaceFromSurfaceTree(const SurfaceTree* st, struct bu_list *vlfree, struct bv_vlblock *vbp, int UNUSED(isocurveres), int UNUSED(gridres))
 {
+    if (!st || !st->Valid())
+	return;
     const BBNode *root = st->getRootNode();
     drawBBNode(st, vlfree, vbp, root);
 }

src/librt/prep.cpp

@@ -75,6 +75,7 @@ rt_i_internal_create(void)
     struct rt_i_internal *ip;
     BU_GET(ip, struct rt_i_internal);
     memset(ip, 0, sizeof(struct rt_i_internal));
+    ip->rti_brep_surface_tree_depth = RT_BREP_DEFAULT_SURFACE_TREE_DEPTH;
 
     for (int i=0; i < RT_DBNHASH; i++) {
 	BU_LIST_INIT(&(ip->rti_solidheads[i]));

src/librt/primitives/brep/brep_local.h

@@ -25,18 +25,204 @@
 #ifndef LIBRT_PRIMITIVES_BREP_BREP_LOCAL_H
 #define LIBRT_PRIMITIVES_BREP_BREP_LOCAL_H
 
+#include <vector>
+
+
+struct brep_edge_span {
+    ON_BezierCurve curve;
+    ON_BoundingBox bbox;
+    ON_Interval edge_domain;
+    int edge_index = -1;
+};
+
+
+struct brep_trim_span {
+    ON_BezierCurve curve;
+    ON_Interval trim_domain;
+};
+
+
+struct brep_edge_trim_cell {
+    ON_Interval edge_domain;
+    ON_Interval trim_domain;
+    size_t edge_span = 0;
+    size_t trim_span_begin = 0;
+    size_t trim_span_count = 0;
+    int trim_index = -1;
+};
+
+
+struct brep_edge_record {
+    ON_BoundingBox bbox;
+    double bbox_coordinate_scale = 0.0;
+    size_t span_begin = 0;
+    size_t span_count = 0;
+    double tolerance = ON_UNSET_VALUE;
+    double model_tolerance = 0.0;
+    double declared_tolerance = ON_UNSET_VALUE;
+    double measured_discrepancy = ON_UNSET_VALUE;
+    double discrepancy_lower_bound = ON_UNSET_VALUE;
+    double discrepancy_upper_bound = ON_UNSET_VALUE;
+    double discrepancy_bound_tolerance = ON_UNSET_VALUE;
+    size_t correspondence_cells = 0;
+    size_t correspondence_depth = 0;
+    size_t discrepancy_bound_cells = 0;
+    size_t discrepancy_bound_depth = 0;
+    size_t trim_cell_begin[2] = {0, 0};
+    size_t trim_cell_count[2] = {0, 0};
+    int edge_index = -1;
+    int face_index[2] = {-1, -1};
+    bool supported = false;
+    bool bbox_coordinate_scale_valid = false;
+    bool discrepancy_measured = false;
+    bool correspondence_screened = false;
+    bool correspondence_supported = false;
+    bool correspondence_exhausted = false;
+    bool discrepancy_bounded = false;
+    bool discrepancy_bound_exhausted = false;
+    bool discrepancy_endpoints_certified = false;
+    bool frame_interval_supported = false;
+    bool discrepancy_sample_authorized = false;
+    int discrepancy_proof_class = 0;
+    bool discrepancy_authorized = false;
+    bool tolerance_inferred = false;
+};
+
+
+struct brep_surface_control_layout {
+    int stride[2] = {0, 0};
+    int order[2] = {0, 0};
+    bool direct = false;
+    bool rational = false;
+};
+
+
+struct brep_surface_span {
+    ON_BezierSurface surface;
+    brep_surface_control_layout control_layout;
+    ON_BoundingBox bbox;
+    double bbox_coordinate_scale = 0.0;
+    ON_Interval surface_domain[2];
+    int face_index = -1;
+    int span_index = -1;
+    unsigned int singular_side_mask = 0;
+    bool control_points_valid = false;
+    bool bbox_coordinate_scale_valid = false;
+};
+
+
+struct brep_face_record {
+    ON_BoundingBox bbox;
+    double bbox_coordinate_scale = 0.0;
+    size_t span_begin = 0;
+    size_t span_count = 0;
+    size_t span_tree_root = 0;
+    int face_index = -1;
+    int nurb_form_status = 0;
+    bool status2_revolution_singular_map = false;
+    bool supported = false;
+    bool bbox_coordinate_scale_valid = false;
+    bool span_tree_available = false;
+};
+
+
+/* A balanced hierarchy over one face's contiguous Bezier spans.  Traversal
+ * retains span order so bounded root collection stays deterministic. */
+struct brep_prepared_span_node {
+    ON_BoundingBox bbox;
+    double bbox_coordinate_scale = 0.0;
+    size_t span_count = 0;
+    bool cullable = false;
+    bool bbox_coordinate_scale_valid = false;
+};
+
+
+/* A balanced hierarchy over contiguous face records.  Keeping the original
+ * order avoids perturbing the bounded prepared-root collection. */
+struct brep_prepared_face_node {
+    ON_BoundingBox bbox;
+    double bbox_coordinate_scale = 0.0;
+    size_t face_begin = 0;
+    size_t span_count = 0;
+    size_t supported_face_count = 0;
+    size_t reparameterized_face_count = 0;
+    size_t unsupported_face_count = 0;
+    size_t left_child = 0;
+    size_t right_child = 0;
+    bool cullable = false;
+    bool leaf = false;
+    bool bbox_coordinate_scale_valid = false;
+};
+
+
+/* A compact immutable hierarchy over BREP face bounds.  Leaves retain their
+ * source index, allowing callers to restore the original face order after
+ * spatial selection. */
+struct brep_face_bvh_node {
+    ON_BoundingBox bbox;
+    double bbox_coordinate_scale = 0.0;
+    size_t left_child = 0;
+    size_t right_child = 0;
+    size_t index = 0;
+    bool leaf = false;
+    bool bbox_coordinate_scale_valid = false;
+};
+
+
+/* One oriented arc of the link obtained by intersecting a sufficiently small
+ * sphere about a manifold vertex with the incident faces.  outgoing is the
+ * tangent of edge_index directed away from the vertex.  face_index carries
+ * this arc from outgoing to the next arc's outgoing tangent, with solid
+ * directions on its left. */
+struct brep_vertex_arc {
+    ON_3dVector outgoing;
+    ON_3dVector outward_normal;
+    double clockwise_sweep = 0.0;
+    int edge_index = -1;
+    int face_index = -1;
+};
+
+
+struct brep_vertex_record {
+    ON_3dPoint point;
+    std::vector<brep_vertex_arc> arcs;
+    int vertex_index = -1;
+    bool planar = false;
+    bool supported = false;
+};
+
 
 /**
  * The b-rep specific data structure for caching the prepared
  * acceleration data structure.
  */
 struct brep_specific {
-    ON_Brep* brep;
-    BrepBoundingVolume* bvh;
-    int is_solid;
-    int plate_mode;
-    int plate_mode_nocos;
-    double plate_mode_thickness;
+    ON_Brep *brep = NULL;
+    BrepBoundingVolume *bvh = NULL;
+    std::vector<brep_face_bvh_node> legacy_face_nodes;
+    size_t legacy_face_root = (size_t)-1;
+    std::vector<const brlcad::CurveTree *> ctrees;
+    std::vector<brep_edge_record> edge_records;
+    std::vector<brep_edge_span> edge_spans;
+    std::vector<brep_edge_trim_cell> edge_trim_cells;
+    std::vector<brep_trim_span> edge_trim_spans;
+    std::vector<brep_face_record> face_records;
+    std::vector<brep_surface_span> surface_spans;
+    std::vector<brep_prepared_span_node> prepared_span_nodes;
+    std::vector<brep_prepared_face_node> prepared_face_nodes;
+    std::vector<brep_face_bvh_node> prepared_face_bvh_nodes;
+    size_t prepared_face_bvh_root = (size_t)-1;
+    std::vector<brep_vertex_record> vertex_records;
+    int surface_tree_depth_limit = 0;
+    int surface_tree_maximum_depth = 0;
+    size_t surface_tree_nodes = 0;
+    size_t surface_tree_leaves = 0;
+    size_t curve_tree_leaves = 0;
+    size_t surface_tree_build_microseconds = 0;
+    int is_solid = 0;
+    int plate_mode = 0;
+    int plate_mode_nocos = 0;
+    double plate_mode_thickness = 0.0;
 };
 
 #endif /* LIBRT_PRIMITIVES_BREP_BREP_LOCAL_H */

src/librt/primitives/bspline/bspline.cpp

@@ -26,9 +26,6 @@
  */
 /** @} */
 
-/* define to use new brep ray-tracing */
-#define CONVERT_TO_BREP 1
-
 /* define to display old nurbs wireframe plot */
 //#define OLD_WIREFRAME 1
 
@@ -50,9 +47,7 @@
 #include "rt/geom.h"
 #include "raytrace.h"
 
-#ifdef CONVERT_TO_BREP
-#  include "brep/defines.h"
-#endif /* CONVERT_TO_BREP */
+#include "brep/defines.h"
 
 #include "../../librt_private.h"
 
@@ -61,39 +56,16 @@
 extern "C" {
 #endif
 
-struct nurb_specific {
-    struct nurb_specific *next;	/* next surface in the solid */
-    struct face_g_snurb *srf;	/* Original surface description */
-    struct bu_list bez_hd;	/* List of Bezier face_g_snurbs */
-};
+extern void rt_nurb_brep(ON_Brep **b, const struct rt_db_internal *ip, const struct bn_tol *tol);
 
-struct nurb_hit {
-    struct nurb_hit *next;
-    struct nurb_hit *prev;
-    fastf_t hit_dist;	/* Distance from the r_pt to surface */
-    point_t hit_point;	/* intersection point */
-    vect_t hit_normal;	/* Surface normal */
-    fastf_t hit_uv[2];	/* Surface parametric u, v */
-    int hit_surfno;	/* index for the surface hit */
-    char *hit_private;	/* Store current nurb root */
-};
-
-#define NULL_HIT (struct nurb_hit *)0
-
-
-#ifdef CONVERT_TO_BREP
-    extern void rt_nurb_brep(ON_Brep **b, const struct rt_db_internal *ip, const struct bn_tol *tol);
-
-    extern int rt_brep_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip);
-    extern void rt_brep_print(const struct soltab *stp);
-    extern int rt_brep_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead);
-    extern void rt_brep_norm(struct hit *hitp, struct soltab *stp, struct xray *rp);
-    extern void rt_brep_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp);
-    extern void rt_brep_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp);
-    extern void rt_brep_free(struct soltab *stp);
-    extern int rt_brep_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct bg_tess_tol *ttol, const struct bn_tol *tol, const struct bview *UNUSED(info));
-
-#endif /* CONVERT_TO_BREP */
+extern int rt_brep_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip);
+extern void rt_brep_print(const struct soltab *stp);
+extern int rt_brep_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead);
+extern void rt_brep_norm(struct hit *hitp, struct soltab *stp, struct xray *rp);
+extern void rt_brep_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp);
+extern void rt_brep_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp);
+extern void rt_brep_free(struct soltab *stp);
+extern int rt_brep_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct bg_tess_tol *ttol, const struct bn_tol *tol, const struct bview *UNUSED(info));
 
 
 static int
@@ -121,8 +93,6 @@ rt_nurb_grans(struct face_g_snurb *srf)
 int
 rt_nurb_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct bn_tol *UNUSED(tol)) {
     int i;
-    struct nurb_specific *nurbs = NULL;
-    struct nurb_specific *next;
     struct rt_nurb_internal *sip;
 
     sip = (struct rt_nurb_internal *) ip->idb_ptr;
@@ -131,49 +101,28 @@ rt_nurb_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct
 
     for (i = 0; i < sip->nsrf; i++) {
 	struct face_g_snurb * s;
-	struct nurb_specific * n;
-
-	BU_GET(n, struct nurb_specific);
-
-	/* Store off the original face_g_snurb */
-	s = nmg_nurb_scopy(sip->srfs[i]);
-	NMG_CK_SNURB(s);
-	nmg_nurb_s_bound(s, s->min_pt, s->max_pt);
-
-	n->srf = s;
-	BU_LIST_INIT(&n->bez_hd);
+	struct bu_list bez_hd;
+	BU_LIST_INIT(&bez_hd);
 
 	/* Grind up the original surf into a list of Bezier face_g_snurbs */
-	(void)nmg_nurb_bezier(&n->bez_hd, sip->srfs[i]);
+	(void)nmg_nurb_bezier(&bez_hd, sip->srfs[i]);
 
 	/* Compute bounds of each Bezier face_g_snurb */
-	for (BU_LIST_FOR(s, face_g_snurb, &n->bez_hd)) {
+	for (BU_LIST_FOR(s, face_g_snurb, &bez_hd)) {
 	    NMG_CK_SNURB(s);
 	    nmg_nurb_s_bound(s, s->min_pt, s->max_pt);
 	    VMINMAX((*min), (*max), s->min_pt);
 	    VMINMAX((*min), (*max), s->max_pt);
 	}
 
-	n->next = nurbs;
-	nurbs = n;
-    }
-    /* zero thickness will get missed by the raytracer */
-    BBOX_NONDEGEN((*min), (*max), SMALL_FASTF);
-
-    for (; nurbs != (struct nurb_specific *)0; nurbs = next) {
-	struct face_g_snurb *s;
-
-	next = nurbs->next;
-
-	/* There is a linked list of surfaces to free for each nurb */
-	while (BU_LIST_WHILE (s, face_g_snurb, &nurbs->bez_hd)) {
+	while (BU_LIST_WHILE(s, face_g_snurb, &bez_hd)) {
 	    NMG_CK_SNURB(s);
 	    BU_LIST_DEQUEUE(&(s->l));
 	    nmg_nurb_free_snurb(s);
 	}
-	nmg_nurb_free_snurb(nurbs->srf);	/* original surf */
-	BU_PUT(nurbs, struct nurb_specific);
     }
+    /* zero thickness will get missed by the raytracer */
+    BBOX_NONDEGEN((*min), (*max), SMALL_FASTF);
 
     return 0;
 }
@@ -187,20 +136,14 @@ rt_nurb_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const struct
 int
 rt_nurb_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
 {
-#ifndef CONVERT_TO_BREP
-    struct nurb_specific *nurbs = NULL;
-#endif
-
     struct rt_nurb_internal *sip;
     const struct bn_tol *tol = &rtip->rti_tol;
 
     sip = (struct rt_nurb_internal *) ip->idb_ptr;
     RT_NURB_CK_MAGIC(sip);
 
-#ifdef CONVERT_TO_BREP
     ON_Brep *brep = ON_Brep::New();
     rt_nurb_brep(&brep, ip, tol);
-    sip->brep = brep;
 
     struct rt_brep_internal bi;
     bi.magic = RT_BREP_INTERNAL_MAGIC;
@@ -211,83 +154,13 @@ rt_nurb_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip)
     di.idb_ptr = (void *)&bi;
 
     return rt_brep_prep(stp, &di, rtip);
-#else
-    fastf_t los = tol->dist;
-
-    if (los < SMALL_FASTF)
-	los = SMALL_FASTF;
-
-    for (i = 0; i < sip->nsrf; i++) {
-	struct face_g_snurb * s;
-	struct nurb_specific * n;
-
-	BU_GET(n, struct nurb_specific);
-
-	/* Store off the original face_g_snurb */
-	s = nmg_nurb_scopy(sip->srfs[i]);
-	NMG_CK_SNURB(s);
-	nmg_nurb_s_bound(s, s->min_pt, s->max_pt);
-
-	n->srf = s;
-	BU_LIST_INIT(&n->bez_hd);
-
-	/* Grind up the original surf into a list of Bezier face_g_snurbs */
-	(void)nmg_nurb_bezier(&n->bez_hd, sip->srfs[i]);
-
-	/* Compute bounds of each Bezier face_g_snurb */
-	for (BU_LIST_FOR(s, face_g_snurb, &n->bez_hd)) {
-	    NMG_CK_SNURB(s);
-	    nmg_nurb_s_bound(s, s->min_pt, s->max_pt);
-	    VMINMAX(stp->st_min, stp->st_max, s->min_pt);
-	    VMINMAX(stp->st_min, stp->st_max, s->max_pt);
-	}
-
-	n->next = nurbs;
-	nurbs = n;
-    }
-
-    stp->st_specific = (void *)nurbs;
-
-    /* zero thickness will get missed by the raytracer */
-    BBOX_NONDEGEN(stp->st_min, stp->st_max, los);
-
-    VADD2SCALE(stp->st_center, stp->st_max, stp->st_min, 0.5);
-    {
-	fastf_t f, dx, dy, dz;
-	dx = (stp->st_max[0] - stp->st_min[0])/2;
-	f = dx;
-	dy = (stp->st_max[1] - stp->st_min[1])/2;
-	if (dy > f)  f = dy;
-	dz = (stp->st_max[2] - stp->st_min[2])/2;
-	if (dz > f)  f = dz;
-	stp->st_aradius = f;
-	stp->st_bradius = sqrt(dx*dx + dy*dy + dz*dz);
-    }
-
-    return 0;
-#endif /* CONVERT_TO_BREP */
 }
 
 
 void
 rt_nurb_print(const struct soltab *stp)
 {
-#ifdef CONVERT_TO_BREP
     return rt_brep_print(stp);
-#else
-    struct nurb_specific *nurb =
-	(struct nurb_specific *)stp->st_specific;
-
-    if (nurb == (struct nurb_specific *)0) {
-	bu_log("rt_nurb_print: no surfaces\n");
-	return;
-    }
-
-    for (; nurb != (struct nurb_specific *)0; nurb = nurb->next) {
-	/* XXX There is a linked list of Bezier surfaces to print here too */
-	nmg_nurb_s_print("NURB", nurb->srf);
-    }
-#endif /* CONVERT_TO_BREP */
 }
 
 
@@ -302,139 +175,7 @@ rt_nurb_print(const struct soltab *stp)
 int
 rt_nurb_shot(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
 {
-#ifdef CONVERT_TO_BREP
     return rt_brep_shot(stp, rp, ap, seghead);
-#else
-    struct nurb_specific *nurb =
-	(struct nurb_specific *)stp->st_specific;
-    struct seg *segp;
-    const struct bn_tol *tol = &ap->a_rt_i->rti_tol;
-    point_t p1, p2, p3, p4;
-    vect_t dir1, dir2;
-    plane_t plane1, plane2;
-    struct nurb_hit * hit;
-    struct nurb_hit hit_list;
-    vect_t invdir;
-    int hit_num;
-    int surface = 0;
-
-    VINVDIR(invdir, rp->r_dir);
-
-    /* Create two orthogonal Planes their intersection contains the
-     * ray so we can project the surface into a 2 dimensional problem
-     */
-
-    bn_vec_ortho(dir1, rp->r_dir);
-    VCROSS(dir2, rp->r_dir, dir1);
-
-    VMOVE(p1, rp->r_pt);
-    VADD2(p2, rp->r_pt, rp->r_dir);
-    VADD2(p3, rp->r_pt, dir1);
-    VADD2(p4, rp->r_pt, dir2);
-
-    /* Note: the equation of the plane in BRL-CAD is
-     * Ax + By + Cz = D represented by [A B C D]
-     */
-    bg_make_plane_3pnts(plane1, p1, p3, p2, tol);
-    bg_make_plane_3pnts(plane2, p1, p2, p4, tol);
-
-    /* make sure that the hit_list is zero */
-
-    hit_list.next = (struct nurb_hit *)0;
-    hit_list.prev = (struct nurb_hit *)0;
-    hit_list.hit_dist = 0;
-    VSET(hit_list.hit_point, 0.0, 0.0, 0.0);
-    VSET(hit_list.hit_normal, 0.0, 0.0, 0.0);
-    hit_list.hit_uv[0] = 	hit_list.hit_uv[1] = 0.0;
-    hit_list.hit_surfno = -1;
-    hit_list.hit_private = (char *)0;
-
-    while (nurb != (struct nurb_specific *) 0) {
-	struct face_g_snurb * s;
-	struct nmg_nurb_uv_hit *hp;
-
-	for (BU_LIST_FOR (s, face_g_snurb, &nurb->bez_hd)) {
-	    if (!rt_in_rpp(rp, invdir, s->min_pt, s->max_pt))
-		continue;
-
-#define UV_TOL 1.0e-6	/* Paper says 1.0e-4 is reasonable for 1k images, not close up */
-	    hp = nmg_nurb_intersect(
-		s, plane1, plane2, UV_TOL, NULL);
-	    while (hp != (struct nmg_nurb_uv_hit *)0) {
-		struct nmg_nurb_uv_hit * o;
-
-		if (RT_G_DEBUG & RT_DEBUG_SPLINE)
-		    bu_log("hit at %d %d sub = %d u = %f v = %f\n",
-			   ap->a_x, ap->a_y, hp->sub, hp->u, hp->v);
-
-		hit = (struct nurb_hit *)
-		    rt_conv_uv(nurb, rp, hp);
-		hit->hit_surfno = surface;
-
-		o = hp;
-		hp = hp->next;
-		bu_free((char *)o,
-			"rt_nurb_shot:nmg_nurb_uv_hit structure");
-
-		rt_nurb_add_hit(&hit_list, hit, tol);
-	    }
-
-	    /* surfaces are incrementally numbered based on the
-	     * face_g_snurb list for all nurbs in this nurb_specific.
-	     */
-	    surface++;
-	}
-	nurb = nurb->next;
-	/* Insert Trimming routines here */
-    }
-
-    /* Convert hits to segments for rt */
-
-    hit_num = 0;
-
-    while (hit_list.next != NULL_HIT) {
-	struct nurb_hit * h1, * h2;
-
-	RT_GET_SEG(segp, ap->a_resource);
-
-	h1 = (struct nurb_hit *) rt_return_nurb_hit(&hit_list);
-	h2 = (struct nurb_hit *) rt_return_nurb_hit(&hit_list);
-
-	segp->seg_stp = stp;
-	segp->seg_in.hit_dist = h1->hit_dist;
-	VMOVE(segp->seg_in.hit_point, h1->hit_point);
-	segp->seg_in.hit_vpriv[0] = h1->hit_uv[0];
-	segp->seg_in.hit_vpriv[1] = h1->hit_uv[1];
-	segp->seg_in.hit_vpriv[2] = 0;
-	segp->seg_in.hit_private = h1->hit_private;
-	segp->seg_in.hit_surfno = h1->hit_surfno;
-	hit_num++;
-
-	if (h2 != NULL_HIT) {
-	    segp->seg_out.hit_dist = h2->hit_dist;
-	    VMOVE(segp->seg_out.hit_point, h2->hit_point);
-	    segp->seg_out.hit_vpriv[0] = h2->hit_uv[0];
-	    segp->seg_out.hit_vpriv[1] = h2->hit_uv[1];
-	    segp->seg_out.hit_private = h2->hit_private;
-	    bu_free((char *)h2, "rt_nurb_shot: nurb hit");
-	    hit_num++;
-	} else {
-	    segp->seg_out.hit_dist = h1->hit_dist + .01;
-	    VJOIN1(segp->seg_out.hit_point,
-		   rp->r_pt, segp->seg_out.hit_dist, rp->r_dir);
-	    segp->seg_out.hit_vpriv[0] = h1->hit_uv[0];
-	    segp->seg_out.hit_vpriv[1] = h1->hit_uv[1];
-	    segp->seg_out.hit_vpriv[2] = 1;
-	    segp->seg_out.hit_private = h1->hit_private;
-	}
-
-	bu_free((char *)h1, "rt_nurb_shot:nurb hit");
-
-	BU_LIST_INSERT(&(seghead->l), &(segp->l));
-    }
-
-    return hit_num;	/* not hit */
-#endif /* CONVERT_TO_BREP */
 }
 
 
@@ -461,27 +202,7 @@ rt_nurb_vshot(struct soltab *stp[], struct xray *rp[], struct seg *segp, int n,
 void
 rt_nurb_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
 {
-#ifdef CONVERT_TO_BREP
     return rt_brep_norm(hitp, stp, rp);
-#else
-/*	struct nurb_specific *nurb =
-	(struct nurb_specific *)stp->st_specific; */
-
-    struct face_g_snurb *n = (struct face_g_snurb *) hitp->hit_private;
-    fastf_t u = hitp->hit_vpriv[0];
-    fastf_t v = hitp->hit_vpriv[1];
-    fastf_t norm[4];
-
-    nmg_nurb_s_norm(n, u, v, norm);
-
-    VMOVE(hitp->hit_normal, norm);
-
-    if (hitp->hit_vpriv[2] == 1) {
-	VREVERSE(hitp->hit_normal, norm);
-    }
-
-    return;
-#endif /* CONVERT_TO_BREP */
 }
 
 
@@ -491,27 +212,7 @@ rt_nurb_norm(struct hit *hitp, struct soltab *stp, struct xray *rp)
 void
 rt_nurb_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
 {
-#ifdef CONVERT_TO_BREP
     return rt_brep_curve(cvp, hitp, stp);
-#else
-/*	struct nurb_specific *nurb =
-	(struct nurb_specific *)stp->st_specific; */
-    struct face_g_snurb *srf = (struct face_g_snurb *) hitp->hit_private;
-    fastf_t u, v;
-
-    if (srf->order[0] <= 2 && srf->order[1] <= 2) {
-	cvp->crv_c1 = cvp->crv_c2 = 0;
-
-	/* any tangent direction */
-	bn_vec_ortho(cvp->crv_pdir, hitp->hit_normal);
-	return;
-    }
-
-    u = hitp->hit_vpriv[0];
-    v = hitp->hit_vpriv[1];
-
-    nmg_nurb_curvature((struct nmg_curvature *)cvp, srf, u, v);
-#endif /* CONVERT_TO_BREP */
 }
 
 
@@ -524,48 +225,14 @@ rt_nurb_curve(struct curvature *cvp, struct hit *hitp, struct soltab *stp)
 void
 rt_nurb_uv(struct application *ap, struct soltab *stp, struct hit *hitp, struct uvcoord *uvp)
 {
-#ifdef CONVERT_TO_BREP
     return rt_brep_uv(ap, stp, hitp, uvp);
-#else
-/*	struct nurb_specific *nurb =
-	(struct nurb_specific *)stp->st_specific; */
-    uvp->uv_u = hitp->hit_vpriv[0];
-    uvp->uv_v = hitp->hit_vpriv[1];
-    return;
-#endif /* CONVERT_TO_BREP */
 }
 
 
 void
 rt_nurb_free(struct soltab *stp)
 {
-#ifdef CONVERT_TO_BREP
     return rt_brep_free(stp);
-#else
-
-    struct nurb_specific *nurb =
-	(struct nurb_specific *)stp->st_specific;
-    struct nurb_specific *next;
-
-    if (nurb == (struct nurb_specific *)0)
-	bu_bomb("rt_nurb_free: no surfaces\n");
-
-    for (; nurb != (struct nurb_specific *)0; nurb = next) {
-	struct face_g_snurb *s;
-
-	next = nurb->next;
-
-	/* There is a linked list of surfaces to free for each nurb */
-	while (BU_LIST_WHILE (s, face_g_snurb, &nurb->bez_hd)) {
-	    NMG_CK_SNURB(s);
-	    BU_LIST_DEQUEUE(&(s->l));
-	    nmg_nurb_free_snurb(s);
-	}
-	nmg_nurb_free_snurb(nurb->srf);	/* original surf */
-	BU_PUT(nurb, struct nurb_specific);
-    }
-    return;
-#endif /* CONVERT_TO_BREP */
 }
 
 
@@ -695,7 +362,7 @@ rt_nurb_plot(struct bu_list *vhead, struct rt_db_internal *ip, const struct bg_t
     return rt_brep_plot(vhead, &di, ttol, tol, NULL);
 #else
     return 0;
-#endif /* CONVERT_TO_BREP */
+#endif /* NEW_WIREFRAME */
 }
 
 

src/librt/primitives/extrude/extrude.c

@@ -174,10 +174,11 @@ rt_extrude_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const str
     vert_count = skt->vert_count;
     /* count how many additional vertices we will need for arc centers */
     for (i = 0; i < skt->curve.count; i++) {
-	struct carc_seg *csg=(struct carc_seg *)skt->curve.segment[i];
+	uint32_t *magic = (uint32_t *)skt->curve.segment[i];
 
-	if (csg->magic != CURVE_CARC_MAGIC)
+	if (*magic != CURVE_CARC_MAGIC)
 	    continue;
+	struct carc_seg *csg=(struct carc_seg *)skt->curve.segment[i];
 
 	if (csg->radius <= 0.0)
 	    continue;
@@ -249,12 +250,13 @@ rt_extrude_bbox(struct rt_db_internal *ip, point_t *min, point_t *max, const str
      * beyond the listed vertices
      */
     for (i = 0; i < skt->curve.count; i++) {
-	struct carc_seg *csg=(struct carc_seg *)skt->curve.segment[i];
-	struct carc_seg *csg_extr=(struct carc_seg *)extr.crv.segment[i];
+	uint32_t *magic = (uint32_t *)skt->curve.segment[i];
 	point_t center;
 
-	if (csg->magic != CURVE_CARC_MAGIC)
+	if (*magic != CURVE_CARC_MAGIC)
 	    continue;
+	struct carc_seg *csg=(struct carc_seg *)skt->curve.segment[i];
+	struct carc_seg *csg_extr=(struct carc_seg *)extr.crv.segment[i];
 
 	if (csg->radius <= 0.0) {
 	    /* full circle */
@@ -450,10 +452,11 @@ rt_extrude_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip
     vert_count = skt->vert_count;
     /* count how many additional vertices we will need for arc centers */
     for (i = 0; i < skt->curve.count; i++) {
-	struct carc_seg *csg=(struct carc_seg *)skt->curve.segment[i];
+	uint32_t *magic = (uint32_t *)skt->curve.segment[i];
 
-	if (csg->magic != CURVE_CARC_MAGIC)
+	if (*magic != CURVE_CARC_MAGIC)
 	    continue;
+	struct carc_seg *csg=(struct carc_seg *)skt->curve.segment[i];
 
 	if (csg->radius <= 0.0)
 	    continue;
@@ -525,12 +528,13 @@ rt_extrude_prep(struct soltab *stp, struct rt_db_internal *ip, struct rt_i *rtip
      * beyond the listed vertices
      */
     for (i = 0; i < skt->curve.count; i++) {
-	struct carc_seg *csg=(struct carc_seg *)skt->curve.segment[i];
-	struct carc_seg *csg_extr=(struct carc_seg *)extr->crv.segment[i];
+	uint32_t *magic = (uint32_t *)skt->curve.segment[i];
 	point_t center;
 
-	if (csg->magic != CURVE_CARC_MAGIC)
+	if (*magic != CURVE_CARC_MAGIC)
 	    continue;
+	struct carc_seg *csg=(struct carc_seg *)skt->curve.segment[i];
+	struct carc_seg *csg_extr=(struct carc_seg *)extr->crv.segment[i];
 
 	if (csg->radius <= 0.0) {
 	    /* full circle */

src/librt/tests/CMakeLists.txt

@@ -188,6 +188,91 @@ brlcad_add_test(NAME rt_primitive_brep COMMAND rt_primitive_brep)
 brlcad_addexec(rt_brep_io brep_io.cpp "librt;libbrep;libbu;Threads::Threads;${M_LIBRARY}" TEST)
 brlcad_add_test(NAME rt_brep_io COMMAND rt_brep_io)
 
+# Directed analytic primitive vs converted-BRep ray intersections
+brlcad_addexec(rt_brep_ray brep_ray.cpp "librt;libbrep;libwdb;libbu;${M_LIBRARY}" TEST)
+brlcad_add_test(NAME rt_brep_ray COMMAND rt_brep_ray --directed-only)
+
+# Mixed collapsed-pole and regular BREP ray events
+brlcad_addexec(rt_brep_mixed_pole brep_mixed_pole.cpp librt;libbrep;libbu;${M_LIBRARY} TEST)
+brlcad_add_test(NAME rt_brep_mixed_pole COMMAND rt_brep_mixed_pole)
+
+# Internal Bezier-span duplicate BREP ray events
+brlcad_addexec(rt_brep_adjacent_span brep_adjacent_span.cpp "librt;libbrep;libbu;${M_LIBRARY}" TEST)
+brlcad_add_test(NAME rt_brep_adjacent_span COMMAND rt_brep_adjacent_span)
+
+# Internal Bezier-span junction BREP ray events
+brlcad_addexec(rt_brep_span_junction brep_span_junction.cpp "librt;libbrep;libbu;${M_LIBRARY}" TEST)
+brlcad_add_test(NAME rt_brep_span_junction COMMAND rt_brep_span_junction)
+foreach(span_junction_face_index 0 1 2 3 4)
+  brlcad_add_test(
+    NAME rt_brep_span_junction_face_${span_junction_face_index}
+    COMMAND rt_brep_span_junction ${span_junction_face_index})
+endforeach()
+
+brlcad_add_test(
+  NAME rt_brep_ray_pole
+  COMMAND rt_brep_ray --pole-only --surface-tree-depth=2
+)
+brlcad_add_test(NAME rt_brep_ray_crofton COMMAND rt_brep_ray --crofton-only)
+brlcad_add_test(
+  NAME rt_brep_ray_regular_stream
+  COMMAND rt_brep_ray --regular-stream-only --surface-tree-depth=0
+)
+brlcad_add_test(
+  NAME rt_brep_ray_crofton_hard
+  COMMAND rt_brep_ray --crofton-hard-only
+)
+brlcad_add_test(NAME rt_brep_ray_seam COMMAND rt_brep_ray --seam-only)
+brlcad_add_test(
+  NAME rt_brep_ray_seam_authority
+  COMMAND rt_brep_ray --seam-authority-only
+)
+brlcad_add_test(
+  NAME rt_brep_ray_throughput
+  COMMAND rt_brep_ray --throughput-only
+)
+brlcad_add_test(
+  NAME rt_brep_ray_restriction
+  COMMAND rt_brep_ray --restriction-only
+)
+brlcad_add_test(NAME rt_brep_ray_endpoint COMMAND rt_brep_ray --endpoint-only)
+brlcad_add_test(
+  NAME rt_brep_ray_nonisoparametric
+  COMMAND rt_brep_ray --nonisoparametric-only
+)
+brlcad_add_test(
+  NAME rt_brep_ray_contact_trim
+  COMMAND rt_brep_ray --contact-trim-only
+)
+brlcad_add_test(
+  NAME rt_brep_ray_defect
+  COMMAND rt_brep_ray --defect-only
+)
+brlcad_add_test(NAME rt_brep_ray_interval COMMAND rt_brep_ray --interval-only)
+brlcad_add_test(
+  NAME rt_brep_ray_local_root
+  COMMAND rt_brep_ray --local-root-only
+)
+brlcad_add_test(NAME rt_brep_ray_cleanup COMMAND rt_brep_ray --cleanup-only)
+brlcad_add_test(
+  NAME rt_brep_ray_grazing_root
+  COMMAND rt_brep_ray --grazing-root-only
+)
+brlcad_add_test(
+  NAME rt_brep_ray_trim_interval
+  COMMAND rt_brep_ray --trim-interval-only
+)
+brlcad_add_test(
+  NAME rt_brep_ray_trim_classifier
+  COMMAND rt_brep_ray --trim-classifier-only
+)
+brlcad_add_test(
+  NAME rt_brep_ray_source_union
+  COMMAND rt_brep_ray --source-union-only
+)
+brlcad_add_test(NAME rt_brep_ray_fold COMMAND rt_brep_ray --fold-only)
+brlcad_add_test(NAME rt_brep_ray_affine COMMAND rt_brep_ray --affine-only)
+
 # stand-alone Crofton estimator tests
 brlcad_addexec(rt_crofton crofton.c "librt;libwdb;libbu;${M_LIBRARY}" TEST)
 brlcad_add_test(NAME rt_crofton COMMAND rt_crofton)
@@ -216,6 +301,7 @@ set(
   arb_intersect.g
   binary_attribute.c
   binunif.c
+  brep_mixed_pole.cpp
   brep_io.cpp
   bound_internal.c
   crofton.c
@@ -224,6 +310,7 @@ set(
   datum_annot_semantics.c
   dsp_raytrace_check.c
   dsp_sa_vol_check.c
+  brep_ray.cpp
   primitive_brep.cpp
   reprep_prune.c
   brep_boolean_tests.g

src/librt/tests/brep_adjacent_span.cpp

@@ -0,0 +1,755 @@
+/*             B R E P _ A D J A C E N T _ S P A N . C P P
+ * BRL-CAD
+ *
+ * Internal Bezier-span BREP ray regression.
+ */
+
+#include "common.h"
+
+#include <cmath>
+#include <cstdio>
+#include <vector>
+
+#include "bu/app.h"
+#include "bu/malloc.h"
+#include "raytrace.h"
+#include "rt/geom.h"
+#include "brep.h"
+#include "brep/cobb.h"
+#include "../librt_private.h"
+
+#include "../primitives/brep/brep_local.h"
+
+static const int ADJACENT_SPAN_SURFACE_TREE_DEPTH = 0;
+static const double ADJACENT_SPAN_RADIUS = 10.0;
+static const double ADJACENT_SPAN_RADIAL_OFFSET = 8.0;
+static const double ADJACENT_SPAN_ORIGIN_DISTANCE = 20.0;
+static const int ADJACENT_SPAN_HIERARCHY_SPAN_COUNT = 16;
+static const double ADJACENT_SPAN_ENTRY_DISTANCE =
+    ADJACENT_SPAN_ORIGIN_DISTANCE - std::sqrt(
+	ADJACENT_SPAN_RADIUS * ADJACENT_SPAN_RADIUS -
+	ADJACENT_SPAN_RADIAL_OFFSET * ADJACENT_SPAN_RADIAL_OFFSET);
+static const double ADJACENT_SPAN_EXIT_DISTANCE =
+    ADJACENT_SPAN_ORIGIN_DISTANCE + std::sqrt(
+	ADJACENT_SPAN_RADIUS * ADJACENT_SPAN_RADIUS -
+	ADJACENT_SPAN_RADIAL_OFFSET * ADJACENT_SPAN_RADIAL_OFFSET);
+
+
+static struct soltab *
+prepare_solid(struct rt_i *rtip, struct rt_db_internal *intern, int type)
+{
+    struct soltab *stp = (struct soltab *)bu_calloc(1, sizeof(struct soltab),
+	"adjacent span ray soltab");
+    stp->l.magic = RT_SOLTAB_MAGIC;
+    stp->l2.magic = RT_SOLTAB2_MAGIC;
+    stp->st_rtip = rtip;
+    struct directory *directory = (struct directory *)bu_calloc(1,
+	sizeof(struct directory), "adjacent span ray directory");
+    directory->d_magic = RT_DIR_MAGIC;
+    directory->d_namep = (char *)"adjacent_span_ray.s";
+    stp->st_dp = directory;
+    stp->st_id = type;
+    stp->st_meth = &OBJ[type];
+
+    if (OBJ[type].ft_prep(stp, intern, rtip)) {
+	if (stp->st_specific && stp->st_meth && stp->st_meth->ft_free)
+	    stp->st_meth->ft_free(stp);
+	bu_free((void *)stp->st_dp, "adjacent span ray directory");
+	bu_free(stp, "adjacent span ray soltab");
+	return NULL;
+    }
+    return stp;
+}
+
+
+static void
+free_solid(struct soltab *stp)
+{
+    if (!stp)
+	return;
+    if (stp->st_meth && stp->st_meth->ft_free)
+	stp->st_meth->ft_free(stp);
+    bu_free((void *)stp->st_dp, "adjacent span ray directory");
+    bu_free(stp, "adjacent span ray soltab");
+}
+
+
+static void
+free_rtip(struct rt_i *rtip, struct resource *resource)
+{
+    if (!rtip || !resource)
+	return;
+    rt_clean_resource_basic(rtip, resource);
+    BU_PTBL_SET(&rtip->rti_resources, 0, NULL);
+    rt_i_destroy(rtip);
+}
+
+
+static ON_Brep *
+make_adjacent_span_sphere()
+{
+    const ON_3dPoint center(0.0, 0.0, 0.0);
+    ON_Brep *brep = ON_Brep_CobbSphereSewn(ADJACENT_SPAN_RADIUS, center);
+    if (!brep)
+	return NULL;
+    for (int face_index = 0; face_index < brep->m_F.Count(); ++face_index) {
+	const int surface_index = brep->m_F[face_index].m_si;
+	ON_NurbsSurface *surface = surface_index >= 0 &&
+	    surface_index < brep->m_S.Count() ? ON_NurbsSurface::Cast(
+	brep->m_S[surface_index]) : NULL;
+	if (!surface) {
+	    delete brep;
+	    return NULL;
+	}
+	const ON_Interval domain = surface->Domain(1);
+	const double span_width = 0.25 * (domain.Max() - domain.Min());
+	if (!domain.IsIncreasing() || !(span_width > 0.0) ||
+		!surface->InsertKnot(1, domain.Min() + span_width) ||
+		!surface->InsertKnot(1, domain.Mid())) {
+	    delete brep;
+	    return NULL;
+	}
+    }
+    return brep;
+}
+
+
+static ON_Brep *
+make_hierarchy_span_sphere()
+{
+    const ON_3dPoint center(0.0, 0.0, 0.0);
+    ON_Brep *brep = ON_Brep_CobbSphereSewn(ADJACENT_SPAN_RADIUS, center);
+    if (!brep)
+	return NULL;
+    for (int face_index = 0; face_index < brep->m_F.Count(); ++face_index) {
+	const int surface_index = brep->m_F[face_index].m_si;
+	ON_NurbsSurface *surface = surface_index >= 0 &&
+	    surface_index < brep->m_S.Count() ? ON_NurbsSurface::Cast(
+	brep->m_S[surface_index]) : NULL;
+	if (!surface) {
+	    delete brep;
+	    return NULL;
+	}
+	const ON_Interval domain = surface->Domain(1);
+	const double span_width = (domain.Max() - domain.Min()) /
+	    ADJACENT_SPAN_HIERARCHY_SPAN_COUNT;
+	if (!domain.IsIncreasing() || !(span_width > 0.0)) {
+	    delete brep;
+	    return NULL;
+	}
+	for (int span = 1; span < ADJACENT_SPAN_HIERARCHY_SPAN_COUNT;
+		++span) {
+	    if (!surface->InsertKnot(1, domain.Min() + span * span_width)) {
+		delete brep;
+		return NULL;
+	    }
+	}
+    }
+    return brep;
+}
+
+
+static struct soltab *
+prepare_brep_solid(struct rt_i *rtip, ON_Brep *brep)
+{
+    if (!rtip || !brep)
+	return NULL;
+    struct rt_brep_internal brep_internal = {};
+    brep_internal.magic = RT_BREP_INTERNAL_MAGIC;
+    brep_internal.brep = brep;
+    struct rt_db_internal database_internal;
+    RT_DB_INTERNAL_INIT(&database_internal);
+    database_internal.idb_major_type = DB5_MAJORTYPE_BRLCAD;
+    database_internal.idb_type = ID_BREP;
+    database_internal.idb_meth = &OBJ[ID_BREP];
+    database_internal.idb_ptr = &brep_internal;
+    return prepare_solid(rtip, &database_internal, ID_BREP);
+}
+
+
+static bool
+bbox_contains(const ON_BoundingBox &outer, const ON_BoundingBox &inner)
+{
+    if (!outer.IsValid() || !inner.IsValid())
+	return false;
+    for (int axis = 0; axis < 3; ++axis) {
+	if (outer.m_min[axis] > inner.m_min[axis] ||
+		outer.m_max[axis] < inner.m_max[axis])
+	    return false;
+    }
+    return true;
+}
+
+
+static bool
+prepared_span_tree_node_valid(const brep_specific &specific,
+    size_t node_index, size_t span_begin, size_t span_count)
+{
+    if (!span_count || node_index >= specific.prepared_span_nodes.size() ||
+	span_begin > specific.surface_spans.size() ||
+	span_count > specific.surface_spans.size() - span_begin)
+	return false;
+    const brep_prepared_span_node &node =
+	specific.prepared_span_nodes[node_index];
+    if (node.span_count != span_count || !node.cullable || !node.bbox.IsValid())
+	return false;
+    if (span_count == 1)
+	return bbox_contains(node.bbox, specific.surface_spans[span_begin].bbox);
+    const size_t left_count = span_count / 2;
+    const size_t left_child = node_index + 1;
+    const size_t right_child = left_child + 2 * left_count - 1;
+    if (right_child >= specific.prepared_span_nodes.size() ||
+	!prepared_span_tree_node_valid(specific, left_child, span_begin,
+	    left_count) ||
+	!prepared_span_tree_node_valid(specific, right_child,
+	    span_begin + left_count, span_count - left_count))
+	return false;
+    return bbox_contains(node.bbox,
+	specific.prepared_span_nodes[left_child].bbox) &&
+	bbox_contains(node.bbox,
+	specific.prepared_span_nodes[right_child].bbox);
+}
+
+
+static bool
+prepared_span_tree_available(const struct soltab *solid)
+{
+    const brep_specific *specific = solid ?
+	static_cast<const brep_specific *>(solid->st_specific) : NULL;
+    if (!specific)
+	return false;
+    size_t multi_span_faces = 0;
+    bool multi_level_tree = false;
+    for (std::vector<brep_face_record>::const_iterator record_it =
+	    specific->face_records.begin();
+	    record_it != specific->face_records.end(); ++record_it) {
+	const brep_face_record &record = *record_it;
+	if (!record.supported || record.span_count < 2)
+	    continue;
+	multi_span_faces++;
+	if (record.span_count >= 3)
+	    multi_level_tree = true;
+	if (!record.span_tree_available ||
+		record.span_tree_root >= specific->prepared_span_nodes.size())
+	    return false;
+	if (!prepared_span_tree_node_valid(*specific, record.span_tree_root,
+		record.span_begin, record.span_count))
+	    return false;
+    }
+    return multi_span_faces > 0 && multi_level_tree;
+}
+
+
+static bool
+disable_prepared_span_tree(struct soltab *solid)
+{
+    brep_specific *specific = solid ?
+	static_cast<brep_specific *>(solid->st_specific) : NULL;
+    if (!specific)
+	return false;
+    for (std::vector<brep_face_record>::iterator record_it =
+	    specific->face_records.begin();
+	    record_it != specific->face_records.end(); ++record_it) {
+	record_it->span_tree_root = 0;
+	record_it->span_tree_available = false;
+    }
+    specific->prepared_span_nodes.clear();
+    return true;
+}
+
+
+static bool
+disable_prepared_face_bvh(struct soltab *solid)
+{
+    brep_specific *specific = solid ?
+	static_cast<brep_specific *>(solid->st_specific) : NULL;
+    if (!specific)
+	return false;
+    specific->prepared_face_bvh_nodes.clear();
+    specific->prepared_face_bvh_root = (size_t)-1;
+    return true;
+}
+
+
+static bool
+disable_prepared_face_tree(struct soltab *solid)
+{
+    brep_specific *specific = solid ?
+	static_cast<brep_specific *>(solid->st_specific) : NULL;
+    if (!specific || !disable_prepared_face_bvh(solid))
+	return false;
+    specific->prepared_face_nodes.clear();
+    return true;
+}
+
+
+static bool
+prepared_face_tree_node_valid(const brep_specific &specific,
+    size_t node_index, size_t face_begin, size_t face_count)
+{
+    if (!face_count || node_index >= specific.prepared_face_nodes.size() ||
+	face_begin > specific.face_records.size() ||
+	face_count > specific.face_records.size() - face_begin)
+	return false;
+    const brep_prepared_face_node &node =
+	specific.prepared_face_nodes[node_index];
+    if (node.face_begin != face_begin)
+	return false;
+    if (face_count == 1) {
+	const brep_face_record &record = specific.face_records[face_begin];
+	const size_t supported = record.supported ? 1 : 0;
+	const size_t reparameterized = record.supported &&
+	    record.nurb_form_status == 2 ? 1 : 0;
+	const size_t unsupported = record.supported ? 0 : 1;
+	const bool cullable = record.supported && record.span_count &&
+	    record.bbox.IsValid();
+	return node.leaf && node.span_count ==
+	    (record.supported ? record.span_count : 0) &&
+	    node.supported_face_count == supported &&
+	    node.reparameterized_face_count == reparameterized &&
+	    node.unsupported_face_count == unsupported &&
+	    node.cullable == cullable && (!record.supported ||
+	    !record.bbox.IsValid() ||
+	    bbox_contains(node.bbox, record.bbox));
+    }
+    const size_t left_count = face_count / 2;
+    const size_t right_count = face_count - left_count;
+    if (node.leaf || node.left_child >= specific.prepared_face_nodes.size() ||
+	node.right_child >= specific.prepared_face_nodes.size() ||
+	!prepared_face_tree_node_valid(specific, node.left_child, face_begin,
+	    left_count) || !prepared_face_tree_node_valid(specific,
+	    node.right_child, face_begin + left_count, right_count))
+	return false;
+    const brep_prepared_face_node &left =
+	specific.prepared_face_nodes[node.left_child];
+    const brep_prepared_face_node &right =
+	specific.prepared_face_nodes[node.right_child];
+    const bool left_complete = !left.supported_face_count || left.cullable;
+    const bool right_complete = !right.supported_face_count || right.cullable;
+    const bool cullable = node.span_count && left_complete && right_complete &&
+	node.bbox.IsValid();
+    return node.span_count == left.span_count + right.span_count &&
+	node.supported_face_count == left.supported_face_count +
+	right.supported_face_count && node.reparameterized_face_count ==
+	left.reparameterized_face_count + right.reparameterized_face_count &&
+	node.unsupported_face_count == left.unsupported_face_count +
+	right.unsupported_face_count && node.cullable == cullable &&
+	(!left.bbox.IsValid() || bbox_contains(node.bbox, left.bbox)) &&
+	(!right.bbox.IsValid() || bbox_contains(node.bbox, right.bbox));
+}
+
+
+static bool
+prepared_face_tree_available(const struct soltab *solid)
+{
+    const brep_specific *specific = solid ?
+	static_cast<const brep_specific *>(solid->st_specific) : NULL;
+    return specific && !specific->face_records.empty() &&
+	!specific->prepared_face_nodes.empty() &&
+	prepared_face_tree_node_valid(*specific, 0, 0,
+	    specific->face_records.size());
+}
+
+
+static bool
+prepared_face_bvh_node_valid(const brep_specific &specific,
+    size_t node_index, std::vector<bool> &seen_faces)
+{
+    if (node_index >= specific.prepared_face_bvh_nodes.size())
+	return false;
+    const brep_face_bvh_node &node =
+	specific.prepared_face_bvh_nodes[node_index];
+    if (!node.bbox.IsValid())
+	return false;
+    if (node.leaf) {
+	if (node.index >= specific.face_records.size() ||
+		seen_faces[node.index])
+	    return false;
+	const brep_face_record &record = specific.face_records[node.index];
+	if (!record.supported || !record.span_count || !record.bbox.IsValid() ||
+		!bbox_contains(node.bbox, record.bbox))
+	    return false;
+	seen_faces[node.index] = true;
+	return true;
+    }
+    if (node.left_child <= node_index || node.right_child <= node_index ||
+	node.left_child >= specific.prepared_face_bvh_nodes.size() ||
+	node.right_child >= specific.prepared_face_bvh_nodes.size())
+	return false;
+    const brep_face_bvh_node &left =
+	specific.prepared_face_bvh_nodes[node.left_child];
+    const brep_face_bvh_node &right =
+	specific.prepared_face_bvh_nodes[node.right_child];
+    return bbox_contains(node.bbox, left.bbox) &&
+	bbox_contains(node.bbox, right.bbox) &&
+	prepared_face_bvh_node_valid(specific, node.left_child, seen_faces) &&
+	prepared_face_bvh_node_valid(specific, node.right_child, seen_faces);
+}
+
+
+static bool
+prepared_face_bvh_available(const struct soltab *solid)
+{
+    const brep_specific *specific = solid ?
+	static_cast<const brep_specific *>(solid->st_specific) : NULL;
+    if (!specific || specific->prepared_face_bvh_nodes.empty() ||
+	specific->prepared_face_bvh_root >=
+	specific->prepared_face_bvh_nodes.size())
+	return false;
+    std::vector<bool> seen_faces(specific->face_records.size(), false);
+    if (!prepared_face_bvh_node_valid(*specific,
+	specific->prepared_face_bvh_root, seen_faces))
+	return false;
+    for (size_t record_index = 0;
+	record_index < specific->face_records.size(); ++record_index) {
+	if (specific->face_records[record_index].supported !=
+		seen_faces[record_index])
+	    return false;
+    }
+    return true;
+}
+
+
+static int
+shoot_trace(struct soltab *stp, struct rt_i *rtip, struct resource *resource,
+    const point_t origin, const vect_t direction,
+    struct rt_brep_shot_trace &trace)
+{
+    struct application application;
+    struct seg seghead;
+    struct xray ray;
+    RT_APPLICATION_INIT(&application);
+    application.a_rt_i = rtip;
+    application.a_resource = resource;
+    VMOVE(ray.r_pt, origin);
+    VMOVE(ray.r_dir, direction);
+    ray.magic = RT_RAY_MAGIC;
+    BU_LIST_INIT(&seghead.l);
+    const int hits = _rt_brep_shot_trace(stp, &ray, &application, &seghead,
+	&trace);
+
+    struct seg *segment;
+    while (BU_LIST_WHILE(segment, seg, &seghead.l)) {
+	BU_LIST_DEQUEUE(&segment->l);
+	RT_FREE_SEG(segment, resource);
+    }
+    return hits;
+}
+
+
+static void
+report_trace(const struct rt_brep_shot_trace &trace)
+{
+    std::printf("  selected/fallback/hits=%zu/%d/%zu boxes/roots/events="
+	"%zu/%zu/%zu stream=%zu/%zu/%zu/%zu/%zu stage=%d\n",
+	trace.prepared_production_selected, trace.prepared_production_fallback,
+	trace.prepared_production_hits, trace.stored_surface_boxes,
+	trace.stored_local_roots, trace.stored_physical_events,
+	trace.physical_event_regular_stream_attempts,
+	trace.physical_event_regular_stream_certified,
+	trace.physical_event_regular_stream_components,
+	trace.physical_event_regular_stream_boxes,
+	trace.physical_event_regular_stream_roots,
+	trace.physical_event_regular_stream_failure_stage);
+    for (size_t root_index = 0; root_index < trace.stored_local_roots;
+	 ++root_index) {
+	const struct rt_brep_trace_local_root &root =
+	    trace.local_roots[root_index];
+	std::printf("  root[%zu] face/span=%d/%d uv=(%.17g,%.17g) t=%.17g "
+	    "normal/direction=%.17g/%d\n", root_index, root.face_index,
+	    root.span_index, root.uv[0], root.uv[1], root.dist,
+	    root.normal_dot, root.direction);
+    }
+    for (size_t box_index = 0; box_index < trace.stored_surface_boxes;
+	 ++box_index) {
+	const struct rt_brep_trace_surface_box &box =
+	    trace.surface_boxes[box_index];
+	std::printf("  box[%zu] face/span=%d/%d uv=[%.17g,%.17g]x"
+	    "[%.17g,%.17g] t=[%.17g,%.17g] state/sign=%d/%d\n",
+	    box_index, box.face_index, box.span_index, box.uv_min[0],
+	    box.uv_max[0], box.uv_min[1], box.uv_max[1], box.t_min,
+	    box.t_max, box.disposition, box.determinant_sign);
+    }
+    for (size_t event_index = 0;
+	 event_index < trace.stored_physical_events; ++event_index) {
+	const struct rt_brep_trace_physical_event &event =
+	    trace.physical_events[event_index];
+	std::printf("  event[%zu] face/span=%d/%d t=[%.17g,%.17g] "
+	    "distance=%.17g root/boxes=%zu/%zu certificate/direction=%d/%d\n",
+	    event_index, event.face_index, event.span_index, event.t_min,
+	    event.t_max, event.dist, event.source_root, event.source_box_count,
+	    event.certificate, event.direction);
+    }
+}
+
+
+static bool
+trace_has_expected_events(const struct rt_brep_shot_trace &trace, int hits,
+    double distance_tolerance)
+{
+    if (hits != 2 || trace.prepared_production_selected != 1 ||
+	trace.prepared_production_fallback != RT_BREP_PREPARED_FALLBACK_NONE ||
+	trace.prepared_production_hits != 2 || trace.physical_event_complete != 1 ||
+	trace.physical_event_regular_stream_attempts != 1 ||
+	trace.physical_event_regular_stream_certified != 1 ||
+	trace.physical_event_regular_stream_components != 2 ||
+	trace.physical_event_regular_stream_boxes !=
+	trace.stored_surface_boxes ||
+	trace.physical_event_regular_stream_roots !=
+	trace.stored_local_roots ||
+	trace.physical_event_regular_stream_failure_stage ||
+	trace.stored_local_roots != 4 || trace.stored_physical_events != 2 ||
+	trace.physical_event_regular != 2 || trace.physical_event_unresolved ||
+	trace.physical_event_direction_mismatches ||
+	trace.physical_event_state_failures ||
+	trace.physical_event_material_segments != 1 ||
+	trace.physical_event_subminimum_contacts ||
+	trace.physical_event_tolerance_ambiguous)
+	return false;
+
+    const double expected_distance[2] = {
+	ADJACENT_SPAN_ENTRY_DISTANCE, ADJACENT_SPAN_EXIT_DISTANCE
+    };
+    const int expected_direction[2] = {
+	RT_BREP_TRACE_ENTERING, RT_BREP_TRACE_LEAVING
+    };
+    size_t source_boxes = 0;
+    for (size_t event_index = 0; event_index < 2; ++event_index) {
+	const struct rt_brep_trace_physical_event &event =
+	    trace.physical_events[event_index];
+	if (event.certificate != RT_BREP_TRACE_EVENT_REGULAR_INTERIOR ||
+		event.source_kind != RT_BREP_TRACE_EVENT_SOURCE_LOCAL_ROOT ||
+		event.source_root >= trace.stored_local_roots ||
+		event.source_box >= trace.stored_surface_boxes ||
+		event.source_box_count < 2 || !event.determinant_sign ||
+		event.direction != expected_direction[event_index] ||
+		fabs(event.dist - expected_distance[event_index]) >
+		distance_tolerance || event.dist < event.t_min ||
+		event.dist > event.t_max)
+	    return false;
+	source_boxes += event.source_box_count;
+    }
+    if (source_boxes != trace.stored_surface_boxes)
+	return false;
+    for (size_t box_index = 0; box_index < trace.stored_surface_boxes;
+	 ++box_index) {
+	const struct rt_brep_trace_surface_box &box =
+	    trace.surface_boxes[box_index];
+	if (box.disposition != RT_BREP_TRACE_BOX_RESOLVED_REGULAR_STREAM ||
+		!box.determinant_sign)
+	    return false;
+    }
+    return true;
+}
+
+
+static bool
+trace_has_hierarchy_pruning(const struct soltab *solid,
+    const struct rt_brep_shot_trace &trace, int hits,
+    double distance_tolerance)
+{
+    const brep_specific *specific = solid ?
+	static_cast<const brep_specific *>(solid->st_specific) : NULL;
+    if (!specific || !trace_has_expected_events(trace, hits,
+	distance_tolerance))
+	return false;
+    size_t largest_face_span_count = 0;
+    for (std::vector<brep_face_record>::const_iterator record_it =
+	    specific->face_records.begin();
+	    record_it != specific->face_records.end(); ++record_it) {
+	const brep_face_record &record = *record_it;
+	if (record.supported && record.span_count > largest_face_span_count)
+	    largest_face_span_count = record.span_count;
+    }
+    const size_t prepared_span_count = specific->surface_spans.size();
+    return largest_face_span_count >= ADJACENT_SPAN_HIERARCHY_SPAN_COUNT &&
+	trace.prepared_surface_spans == prepared_span_count &&
+	trace.candidate_surface_spans + trace.excluded_surface_spans ==
+	prepared_span_count && trace.excluded_surface_spans >
+	3 * prepared_span_count / 4 && trace.candidate_surface_spans <
+	prepared_span_count / 4 && trace.surface_coefficient_expansion_avoided >=
+	trace.excluded_surface_spans;
+}
+
+
+int
+main(int argc, char **argv)
+{
+    bu_setprogname(argv[0]);
+    struct rt_i *rtip = rt_dirbuild_inmem(NULL, 0, NULL, 0);
+    if (!rtip) {
+	std::printf("FAIL: in-memory raytrace setup\n");
+	return 1;
+    }
+    rtip->rti_tol.magic = BN_TOL_MAGIC;
+    rtip->rti_tol.dist = 0.0005;
+    rtip->rti_tol.dist_sq = rtip->rti_tol.dist * rtip->rti_tol.dist;
+    rtip->rti_tol.perp = 1.0e-6;
+    rtip->rti_tol.para = 1.0 - rtip->rti_tol.perp;
+    if (!_rt_brep_set_surface_tree_depth(rtip,
+	ADJACENT_SPAN_SURFACE_TREE_DEPTH)) {
+	std::printf("FAIL: BREP surface tree setup\n");
+	rt_i_destroy(rtip);
+	return 1;
+    }
+
+    struct resource resource = {};
+    rt_init_resource(&resource, 0, rtip);
+    ON_Brep *brep = make_adjacent_span_sphere();
+    if (!brep || !brep->IsSolid()) {
+	std::printf("FAIL: adjacent-span sphere construction\n");
+	delete brep;
+	free_rtip(rtip, &resource);
+	return 1;
+    }
+    struct soltab *solid = prepare_brep_solid(rtip, brep);
+    if (!solid) {
+	std::printf("FAIL: adjacent-span sphere preparation\n");
+	delete brep;
+	free_rtip(rtip, &resource);
+	return 1;
+    }
+
+    const bool span_tree_valid = prepared_span_tree_available(solid) &&
+	prepared_face_tree_available(solid);
+    if (!span_tree_valid)
+	std::printf("FAIL: adjacent-span hierarchy preparation\n");
+
+    const point_t forward_origin = {
+	ADJACENT_SPAN_RADIAL_OFFSET, 0.0,
+	-ADJACENT_SPAN_ORIGIN_DISTANCE
+    };
+    const vect_t forward_direction = {0.0, 0.0, 1.0};
+    const point_t reverse_origin = {
+	ADJACENT_SPAN_RADIAL_OFFSET, 0.0,
+	ADJACENT_SPAN_ORIGIN_DISTANCE
+    };
+    const vect_t reverse_direction = {0.0, 0.0, -1.0};
+    const double distance_tolerance = std::max(1.0e-9,
+	(double)rtip->rti_tol.dist);
+
+    struct rt_brep_shot_trace forward_trace = {};
+    const int forward_hits = shoot_trace(solid, rtip, &resource,
+	forward_origin, forward_direction, forward_trace);
+    const bool forward_valid = trace_has_expected_events(forward_trace,
+	forward_hits, distance_tolerance);
+    if (!forward_valid) {
+	std::printf("FAIL: forward adjacent-span stream\n");
+	report_trace(forward_trace);
+    }
+
+    struct rt_brep_shot_trace reverse_trace = {};
+    const int reverse_hits = shoot_trace(solid, rtip, &resource,
+	reverse_origin, reverse_direction, reverse_trace);
+    const bool reverse_valid = trace_has_expected_events(reverse_trace,
+	reverse_hits, distance_tolerance);
+    if (!reverse_valid) {
+	std::printf("FAIL: reverse adjacent-span stream\n");
+	report_trace(reverse_trace);
+    }
+
+    bool hierarchy_valid = false;
+    ON_Brep *hierarchy_brep = make_hierarchy_span_sphere();
+    if (hierarchy_brep && hierarchy_brep->IsSolid()) {
+	struct soltab *hierarchy_solid = prepare_brep_solid(rtip, hierarchy_brep);
+	if (hierarchy_solid) {
+	    const bool hierarchy_tree_valid =
+		prepared_span_tree_available(hierarchy_solid) &&
+		prepared_face_tree_available(hierarchy_solid) &&
+		prepared_face_bvh_available(hierarchy_solid);
+	    struct rt_brep_shot_trace hierarchy_forward_trace = {};
+	    const int hierarchy_forward_hits = shoot_trace(hierarchy_solid, rtip,
+		&resource, forward_origin, forward_direction,
+		hierarchy_forward_trace);
+	    struct rt_brep_shot_trace hierarchy_reverse_trace = {};
+	    const int hierarchy_reverse_hits = shoot_trace(hierarchy_solid, rtip,
+		&resource, reverse_origin, reverse_direction,
+		hierarchy_reverse_trace);
+	    const bool hierarchy_pruned = hierarchy_tree_valid &&
+		trace_has_hierarchy_pruning(hierarchy_solid,
+		    hierarchy_forward_trace, hierarchy_forward_hits,
+		    distance_tolerance) && trace_has_hierarchy_pruning(hierarchy_solid,
+		    hierarchy_reverse_trace, hierarchy_reverse_hits,
+		    distance_tolerance);
+	    bool spatial_fallback_valid = false;
+	    if (disable_prepared_face_bvh(hierarchy_solid)) {
+		const bool spatial_fallback_tree_valid =
+		    prepared_span_tree_available(hierarchy_solid) &&
+		    prepared_face_tree_available(hierarchy_solid) &&
+		    !prepared_face_bvh_available(hierarchy_solid);
+		struct rt_brep_shot_trace fallback_forward_trace = {};
+		const int fallback_forward_hits = shoot_trace(hierarchy_solid,
+		    rtip, &resource, forward_origin, forward_direction,
+		    fallback_forward_trace);
+		struct rt_brep_shot_trace fallback_reverse_trace = {};
+		const int fallback_reverse_hits = shoot_trace(hierarchy_solid,
+		    rtip, &resource, reverse_origin, reverse_direction,
+		    fallback_reverse_trace);
+		spatial_fallback_valid = spatial_fallback_tree_valid &&
+		    trace_has_hierarchy_pruning(hierarchy_solid, fallback_forward_trace,
+		    fallback_forward_hits, distance_tolerance) &&
+		    trace_has_hierarchy_pruning(hierarchy_solid, fallback_reverse_trace,
+		    fallback_reverse_hits, distance_tolerance);
+	    }
+	    bool face_fallback_valid = false;
+	    if (disable_prepared_face_tree(hierarchy_solid)) {
+		const bool face_fallback_tree_valid =
+		    prepared_span_tree_available(hierarchy_solid) &&
+		    !prepared_face_tree_available(hierarchy_solid) &&
+		    !prepared_face_bvh_available(hierarchy_solid);
+		struct rt_brep_shot_trace fallback_forward_trace = {};
+		const int fallback_forward_hits = shoot_trace(hierarchy_solid,
+		    rtip, &resource, forward_origin, forward_direction,
+		    fallback_forward_trace);
+		struct rt_brep_shot_trace fallback_reverse_trace = {};
+		const int fallback_reverse_hits = shoot_trace(hierarchy_solid,
+		    rtip, &resource, reverse_origin, reverse_direction,
+		    fallback_reverse_trace);
+		face_fallback_valid = face_fallback_tree_valid &&
+		    trace_has_hierarchy_pruning(hierarchy_solid, fallback_forward_trace,
+		    fallback_forward_hits, distance_tolerance) &&
+		    trace_has_hierarchy_pruning(hierarchy_solid, fallback_reverse_trace,
+		    fallback_reverse_hits, distance_tolerance);
+	    }
+	    bool span_fallback_valid = false;
+	    if (disable_prepared_span_tree(hierarchy_solid)) {
+		struct rt_brep_shot_trace fallback_forward_trace = {};
+		const int fallback_forward_hits = shoot_trace(hierarchy_solid,
+		    rtip, &resource, forward_origin, forward_direction,
+		    fallback_forward_trace);
+		struct rt_brep_shot_trace fallback_reverse_trace = {};
+		const int fallback_reverse_hits = shoot_trace(hierarchy_solid,
+		    rtip, &resource, reverse_origin, reverse_direction,
+		    fallback_reverse_trace);
+		span_fallback_valid = trace_has_expected_events(fallback_forward_trace,
+		    fallback_forward_hits, distance_tolerance) &&
+		    trace_has_expected_events(fallback_reverse_trace,
+		    fallback_reverse_hits, distance_tolerance);
+	    }
+	    hierarchy_valid = hierarchy_pruned && spatial_fallback_valid &&
+		face_fallback_valid && span_fallback_valid;
+	    free_solid(hierarchy_solid);
+	} else {
+	    delete hierarchy_brep;
+	}
+    } else {
+	delete hierarchy_brep;
+    }
+
+    free_solid(solid);
+    free_rtip(rtip, &resource);
+    return span_tree_valid && forward_valid && reverse_valid &&
+	hierarchy_valid ? 0 : 1;
+}
+
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C++
+ * c-basic-offset: 4
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */

src/librt/tests/brep_span_junction.cpp

@@ -0,0 +1,616 @@
+/*             B R E P _ S P A N _ J U N C T I O N . C P P
+ * BRL-CAD
+ *
+ * Internal Bezier-span junction BREP ray regression.
+ */
+
+#include "common.h"
+
+#include <cerrno>
+#include <climits>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+
+#include "bu/app.h"
+#include "bu/malloc.h"
+#include "raytrace.h"
+#include "rt/geom.h"
+#include "brep.h"
+#include "brep/cobb.h"
+#include "../librt_private.h"
+
+
+static const int SPAN_JUNCTION_SURFACE_TREE_DEPTH = 0;
+static const int SPAN_JUNCTION_DEFAULT_FACE_INDEX = 5;
+static const double SPAN_JUNCTION_RADIUS = 10.0;
+static const double SPAN_JUNCTION_RAY_CLEARANCE = 20.0;
+static const double SPAN_JUNCTION_OBLIQUE_TANGENT_FRACTION = 0.25;
+static const double SPAN_JUNCTION_BOUNDARY_TRANSVERSE_FRACTION = 0.25;
+static const size_t SPAN_JUNCTION_STREAM_EVENTS = 2;
+static const size_t SPAN_JUNCTION_JUNCTION_BOXES = 4;
+static const size_t SPAN_JUNCTION_ADJACENT_BOXES = 2;
+static const size_t SPAN_JUNCTION_SINGLETON_BOXES = 1;
+static const double SPAN_JUNCTION_ORIGIN_RADIUS =
+    SPAN_JUNCTION_RADIUS + SPAN_JUNCTION_RAY_CLEARANCE;
+static const double SPAN_JUNCTION_ENTRY_DISTANCE =
+    SPAN_JUNCTION_RAY_CLEARANCE;
+static const double SPAN_JUNCTION_EXIT_DISTANCE =
+    SPAN_JUNCTION_ORIGIN_RADIUS + SPAN_JUNCTION_RADIUS;
+
+
+static struct soltab *
+prepare_solid(struct rt_i *rtip, struct rt_db_internal *intern, int type)
+{
+    struct soltab *stp = (struct soltab *)bu_calloc(1, sizeof(struct soltab),
+	"span junction ray soltab");
+    stp->l.magic = RT_SOLTAB_MAGIC;
+    stp->l2.magic = RT_SOLTAB2_MAGIC;
+    stp->st_rtip = rtip;
+    struct directory *directory = (struct directory *)bu_calloc(1,
+	sizeof(struct directory), "span junction ray directory");
+    directory->d_magic = RT_DIR_MAGIC;
+    directory->d_namep = (char *)"span_junction_ray.s";
+    stp->st_dp = directory;
+    stp->st_id = type;
+    stp->st_meth = &OBJ[type];
+
+    if (OBJ[type].ft_prep(stp, intern, rtip)) {
+	if (stp->st_specific && stp->st_meth && stp->st_meth->ft_free)
+	    stp->st_meth->ft_free(stp);
+	bu_free((void *)stp->st_dp, "span junction ray directory");
+	bu_free(stp, "span junction ray soltab");
+	return NULL;
+    }
+    return stp;
+}
+
+
+static void
+free_solid(struct soltab *stp)
+{
+    if (!stp)
+	return;
+    if (stp->st_meth && stp->st_meth->ft_free)
+	stp->st_meth->ft_free(stp);
+    bu_free((void *)stp->st_dp, "span junction ray directory");
+    bu_free(stp, "span junction ray soltab");
+}
+
+
+static void
+free_rtip(struct rt_i *rtip, struct resource *resource)
+{
+    if (!rtip || !resource)
+	return;
+    rt_clean_resource_basic(rtip, resource);
+    BU_PTBL_SET(&rtip->rti_resources, 0, NULL);
+    rt_i_destroy(rtip);
+}
+
+
+static ON_Brep *
+make_span_junction_sphere()
+{
+    const ON_3dPoint center(0.0, 0.0, 0.0);
+    ON_Brep *brep = ON_Brep_CobbSphereSewn(SPAN_JUNCTION_RADIUS, center);
+    if (!brep)
+	return NULL;
+    for (int face_index = 0; face_index < brep->m_F.Count(); ++face_index) {
+	const int surface_index = brep->m_F[face_index].m_si;
+	ON_NurbsSurface *surface = surface_index >= 0 &&
+	    surface_index < brep->m_S.Count() ? ON_NurbsSurface::Cast(
+	brep->m_S[surface_index]) : NULL;
+	if (!surface) {
+	    delete brep;
+	    return NULL;
+	}
+	for (int direction = 0; direction < 2; ++direction) {
+	    const ON_Interval domain = surface->Domain(direction);
+	    if (!domain.IsIncreasing() ||
+		!surface->InsertKnot(direction, domain.Mid())) {
+		delete brep;
+		return NULL;
+	    }
+	}
+    }
+	return brep;
+}
+
+
+static ON_NurbsSurface *
+span_junction_face_surface(ON_Brep *brep, int face_index)
+{
+    if (!brep || face_index < 0 || face_index >= brep->m_F.Count())
+	return NULL;
+    const int surface_index = brep->m_F[face_index].m_si;
+    return surface_index >= 0 && surface_index < brep->m_S.Count() ?
+	ON_NurbsSurface::Cast(brep->m_S[surface_index]) : NULL;
+}
+
+
+static bool
+span_junction_radial(ON_Brep *brep, int face_index, ON_3dVector &radial)
+{
+    const ON_3dPoint center(0.0, 0.0, 0.0);
+    ON_NurbsSurface *surface = span_junction_face_surface(brep, face_index);
+    if (!surface)
+	return false;
+    const ON_3dPoint junction = surface->PointAt(surface->Domain(0).Mid(),
+	surface->Domain(1).Mid());
+    radial = junction - center;
+    return junction.IsValid() && radial.Unitize();
+}
+
+
+static bool
+span_junction_boundary_radial(ON_Brep *brep, int face_index,
+    double transverse_fraction, ON_3dVector &radial)
+{
+    const ON_3dPoint center(0.0, 0.0, 0.0);
+	if (!(transverse_fraction > 0.0) || !(transverse_fraction < 1.0))
+	return false;
+    ON_NurbsSurface *surface = span_junction_face_surface(brep, face_index);
+    if (!surface)
+	return false;
+    const ON_Interval u_domain = surface->Domain(0);
+    const ON_Interval v_domain = surface->Domain(1);
+    if (!u_domain.IsIncreasing() || !v_domain.IsIncreasing())
+	return false;
+    const ON_3dPoint boundary = surface->PointAt(u_domain.Mid(),
+	v_domain.ParameterAt(transverse_fraction));
+    radial = boundary - center;
+    return boundary.IsValid() && radial.Unitize();
+}
+
+
+static struct soltab *
+prepare_brep_solid(struct rt_i *rtip, ON_Brep *brep)
+{
+    if (!rtip || !brep)
+	return NULL;
+    struct rt_brep_internal brep_internal = {};
+    brep_internal.magic = RT_BREP_INTERNAL_MAGIC;
+    brep_internal.brep = brep;
+    struct rt_db_internal database_internal;
+    RT_DB_INTERNAL_INIT(&database_internal);
+    database_internal.idb_major_type = DB5_MAJORTYPE_BRLCAD;
+    database_internal.idb_type = ID_BREP;
+    database_internal.idb_meth = &OBJ[ID_BREP];
+    database_internal.idb_ptr = &brep_internal;
+    return prepare_solid(rtip, &database_internal, ID_BREP);
+}
+
+
+static int
+shoot_trace(struct soltab *stp, struct rt_i *rtip, struct resource *resource,
+    const point_t origin, const vect_t direction,
+    struct rt_brep_shot_trace &trace)
+{
+    struct application application;
+    struct seg seghead;
+    struct xray ray;
+    RT_APPLICATION_INIT(&application);
+    application.a_rt_i = rtip;
+    application.a_resource = resource;
+    VMOVE(ray.r_pt, origin);
+    VMOVE(ray.r_dir, direction);
+    ray.magic = RT_RAY_MAGIC;
+    BU_LIST_INIT(&seghead.l);
+    const int hits = _rt_brep_shot_trace(stp, &ray, &application, &seghead,
+	&trace);
+
+    struct seg *segment;
+    while (BU_LIST_WHILE(segment, seg, &seghead.l)) {
+	BU_LIST_DEQUEUE(&segment->l);
+	RT_FREE_SEG(segment, resource);
+    }
+    return hits;
+}
+
+
+static void
+report_trace(const struct rt_brep_shot_trace &trace)
+{
+    std::printf("  selected/fallback/hits=%zu/%d/%zu boxes/roots/events="
+	"%zu/%zu/%zu stream=%zu/%zu/%zu/%zu/%zu stage=%d\n",
+	trace.prepared_production_selected, trace.prepared_production_fallback,
+	trace.prepared_production_hits, trace.stored_surface_boxes,
+	trace.stored_local_roots, trace.stored_physical_events,
+	trace.physical_event_regular_stream_attempts,
+	trace.physical_event_regular_stream_certified,
+	trace.physical_event_regular_stream_components,
+	trace.physical_event_regular_stream_boxes,
+	trace.physical_event_regular_stream_roots,
+	trace.physical_event_regular_stream_failure_stage);
+    for (size_t root_index = 0; root_index < trace.stored_local_roots;
+	 ++root_index) {
+	const struct rt_brep_trace_local_root &root =
+	    trace.local_roots[root_index];
+	std::printf("  root[%zu] face/span=%d/%d uv=(%.17g,%.17g) t=%.17g "
+	    "normal/direction=%.17g/%d\n", root_index, root.face_index,
+	    root.span_index, root.uv[0], root.uv[1], root.dist,
+	    root.normal_dot, root.direction);
+    }
+    for (size_t box_index = 0; box_index < trace.stored_surface_boxes;
+	 ++box_index) {
+	const struct rt_brep_trace_surface_box &box =
+	    trace.surface_boxes[box_index];
+	std::printf("  box[%zu] face/span=%d/%d uv=[%.17g,%.17g]x"
+	    "[%.17g,%.17g] t=[%.17g,%.17g] state/sign=%d/%d\n",
+	    box_index, box.face_index, box.span_index, box.uv_min[0],
+	    box.uv_max[0], box.uv_min[1], box.uv_max[1], box.t_min,
+	    box.t_max, box.disposition, box.determinant_sign);
+    }
+    for (size_t event_index = 0;
+	 event_index < trace.stored_physical_events; ++event_index) {
+	const struct rt_brep_trace_physical_event &event =
+	    trace.physical_events[event_index];
+	std::printf("  event[%zu] face/span=%d/%d t=[%.17g,%.17g] "
+	    "distance=%.17g root/boxes=%zu/%zu certificate/direction=%d/%d\n",
+	    event_index, event.face_index, event.span_index, event.t_min,
+	    event.t_max, event.dist, event.source_root, event.source_box_count,
+	    event.certificate, event.direction);
+    }
+}
+
+
+static bool
+trace_has_expected_events(const struct rt_brep_shot_trace &trace, int hits,
+    double distance_tolerance)
+{
+    if (hits != 2 || trace.prepared_production_selected != 1 ||
+	trace.prepared_production_fallback != RT_BREP_PREPARED_FALLBACK_NONE ||
+	trace.prepared_production_hits != 2 || trace.physical_event_complete != 1 ||
+	trace.physical_event_regular_stream_attempts != 1 ||
+	trace.physical_event_regular_stream_certified != 1 ||
+	trace.physical_event_regular_stream_components != 2 ||
+	trace.physical_event_regular_stream_boxes !=
+	trace.stored_surface_boxes ||
+	trace.physical_event_regular_stream_roots !=
+	trace.stored_local_roots ||
+	trace.physical_event_regular_stream_failure_stage ||
+	trace.stored_local_roots != 8 || trace.stored_physical_events != 2 ||
+	trace.physical_event_regular != 2 || trace.physical_event_unresolved ||
+	trace.physical_event_direction_mismatches ||
+	trace.physical_event_state_failures ||
+	trace.physical_event_material_segments != 1 ||
+	trace.physical_event_subminimum_contacts ||
+	trace.physical_event_tolerance_ambiguous)
+	return false;
+
+    const double expected_distance[2] = {
+	SPAN_JUNCTION_ENTRY_DISTANCE, SPAN_JUNCTION_EXIT_DISTANCE
+    };
+    const int expected_direction[2] = {
+	RT_BREP_TRACE_ENTERING, RT_BREP_TRACE_LEAVING
+    };
+    size_t source_boxes = 0;
+    for (size_t event_index = 0; event_index < 2; ++event_index) {
+	const struct rt_brep_trace_physical_event &event =
+	    trace.physical_events[event_index];
+	if (event.certificate != RT_BREP_TRACE_EVENT_REGULAR_INTERIOR ||
+		event.source_kind != RT_BREP_TRACE_EVENT_SOURCE_LOCAL_ROOT ||
+		event.source_root >= trace.stored_local_roots ||
+		event.source_box >= trace.stored_surface_boxes ||
+		event.source_box_count < 4 || !event.determinant_sign ||
+		event.direction != expected_direction[event_index] ||
+		fabs(event.dist - expected_distance[event_index]) >
+		distance_tolerance || event.dist < event.t_min ||
+		event.dist > event.t_max)
+	    return false;
+	source_boxes += event.source_box_count;
+    }
+    if (source_boxes != trace.stored_surface_boxes)
+	return false;
+    for (size_t box_index = 0; box_index < trace.stored_surface_boxes;
+	 ++box_index) {
+	const struct rt_brep_trace_surface_box &box =
+	    trace.surface_boxes[box_index];
+	if (box.disposition != RT_BREP_TRACE_BOX_RESOLVED_REGULAR_STREAM ||
+		!box.determinant_sign)
+	    return false;
+    }
+    return true;
+}
+
+
+static bool
+trace_has_span_grid_events(const struct rt_brep_shot_trace &trace, int hits,
+    size_t forward_box_count, size_t reverse_box_count, bool reverse)
+{
+    if (!forward_box_count || !reverse_box_count)
+	return false;
+    const size_t local_root_count = forward_box_count + reverse_box_count;
+    if (hits != 2 || trace.prepared_production_selected != 1 ||
+	trace.prepared_production_fallback != RT_BREP_PREPARED_FALLBACK_NONE ||
+	trace.prepared_production_hits != 2 || trace.physical_event_complete != 1 ||
+	trace.physical_event_regular_stream_attempts != 1 ||
+	trace.physical_event_regular_stream_certified != 1 ||
+	trace.physical_event_regular_stream_components !=
+	SPAN_JUNCTION_STREAM_EVENTS ||
+	trace.physical_event_regular_stream_boxes !=
+	trace.stored_surface_boxes ||
+	trace.physical_event_regular_stream_roots !=
+	trace.stored_local_roots ||
+	trace.physical_event_regular_stream_failure_stage ||
+	trace.stored_surface_boxes != local_root_count ||
+	trace.stored_local_roots != local_root_count ||
+	trace.stored_physical_events != SPAN_JUNCTION_STREAM_EVENTS ||
+	trace.physical_event_regular != SPAN_JUNCTION_STREAM_EVENTS ||
+	trace.physical_event_unresolved ||
+	trace.physical_event_direction_mismatches ||
+	trace.physical_event_state_failures ||
+	trace.physical_event_material_segments != 1 ||
+	trace.physical_event_subminimum_contacts ||
+	trace.physical_event_tolerance_ambiguous)
+	return false;
+
+    const size_t expected_boxes[SPAN_JUNCTION_STREAM_EVENTS] = {
+	reverse ? reverse_box_count : forward_box_count,
+	reverse ? forward_box_count : reverse_box_count
+    };
+    const int expected_direction[SPAN_JUNCTION_STREAM_EVENTS] = {
+	RT_BREP_TRACE_ENTERING, RT_BREP_TRACE_LEAVING
+    };
+    size_t source_boxes = 0;
+    for (size_t event_index = 0;
+	event_index < SPAN_JUNCTION_STREAM_EVENTS; ++event_index) {
+	const struct rt_brep_trace_physical_event &event =
+	    trace.physical_events[event_index];
+	if (event.certificate != RT_BREP_TRACE_EVENT_REGULAR_INTERIOR ||
+		event.source_kind != RT_BREP_TRACE_EVENT_SOURCE_LOCAL_ROOT ||
+		event.source_root >= trace.stored_local_roots ||
+		event.source_box >= trace.stored_surface_boxes ||
+		event.source_box_count != expected_boxes[event_index] ||
+		!event.determinant_sign ||
+		event.direction != expected_direction[event_index] ||
+		event.dist < event.t_min || event.dist > event.t_max)
+	    return false;
+	source_boxes += event.source_box_count;
+    }
+    if (source_boxes != trace.stored_surface_boxes)
+	return false;
+    for (size_t box_index = 0; box_index < trace.stored_surface_boxes;
+	 ++box_index) {
+	const struct rt_brep_trace_surface_box &box =
+	    trace.surface_boxes[box_index];
+	if (box.disposition != RT_BREP_TRACE_BOX_RESOLVED_REGULAR_STREAM ||
+		!box.determinant_sign)
+	    return false;
+    }
+    return true;
+}
+
+
+static bool
+trace_oblique_span_grid(struct soltab *solid, struct rt_i *rtip,
+    struct resource *resource, const ON_3dVector &radial,
+    ON_3dVector tangent, size_t forward_box_count,
+    size_t reverse_box_count)
+{
+    const double projection = tangent.x * radial.x + tangent.y * radial.y +
+	tangent.z * radial.z;
+    tangent.x -= projection * radial.x;
+    tangent.y -= projection * radial.y;
+    tangent.z -= projection * radial.z;
+    if (!tangent.Unitize())
+	return false;
+    ON_3dVector forward_direction(
+	-radial.x + SPAN_JUNCTION_OBLIQUE_TANGENT_FRACTION * tangent.x,
+	-radial.y + SPAN_JUNCTION_OBLIQUE_TANGENT_FRACTION * tangent.y,
+	-radial.z + SPAN_JUNCTION_OBLIQUE_TANGENT_FRACTION * tangent.z);
+    if (!forward_direction.Unitize())
+	return false;
+    const double radial_dot = radial.x * forward_direction.x +
+	radial.y * forward_direction.y + radial.z * forward_direction.z;
+    const double exit_distance = SPAN_JUNCTION_RAY_CLEARANCE -
+	2.0 * SPAN_JUNCTION_RADIUS * radial_dot;
+    if (!std::isfinite(exit_distance) ||
+	exit_distance <= SPAN_JUNCTION_RAY_CLEARANCE)
+	return false;
+
+    point_t forward_origin;
+    vect_t forward_ray;
+    VSET(forward_origin, SPAN_JUNCTION_RADIUS * radial.x -
+	SPAN_JUNCTION_RAY_CLEARANCE * forward_direction.x,
+	SPAN_JUNCTION_RADIUS * radial.y -
+	SPAN_JUNCTION_RAY_CLEARANCE * forward_direction.y,
+	SPAN_JUNCTION_RADIUS * radial.z -
+	SPAN_JUNCTION_RAY_CLEARANCE * forward_direction.z);
+    VSET(forward_ray, forward_direction.x, forward_direction.y,
+	forward_direction.z);
+    struct rt_brep_shot_trace forward_trace = {};
+    const int forward_hits = shoot_trace(solid, rtip, resource,
+	forward_origin, forward_ray, forward_trace);
+    const bool forward_valid = trace_has_span_grid_events(forward_trace,
+	forward_hits, forward_box_count, reverse_box_count, false);
+    if (!forward_valid)
+	report_trace(forward_trace);
+
+    point_t reverse_origin;
+    vect_t reverse_ray;
+    VSET(reverse_origin, SPAN_JUNCTION_RADIUS * radial.x +
+	exit_distance * forward_direction.x,
+	SPAN_JUNCTION_RADIUS * radial.y + exit_distance * forward_direction.y,
+	SPAN_JUNCTION_RADIUS * radial.z + exit_distance * forward_direction.z);
+    VSET(reverse_ray, -forward_direction.x, -forward_direction.y,
+	-forward_direction.z);
+    struct rt_brep_shot_trace reverse_trace = {};
+    const int reverse_hits = shoot_trace(solid, rtip, resource,
+	reverse_origin, reverse_ray, reverse_trace);
+    const bool reverse_valid = trace_has_span_grid_events(reverse_trace,
+	reverse_hits, forward_box_count, reverse_box_count, true);
+    if (!reverse_valid)
+	report_trace(reverse_trace);
+    return forward_valid && reverse_valid;
+}
+
+
+static ON_3dVector
+span_junction_tangent_axis(const ON_3dVector &radial)
+{
+    ON_3dVector tangent(1.0, 0.0, 0.0);
+    if (fabs(radial.y) < fabs(radial.x) &&
+	fabs(radial.y) <= fabs(radial.z))
+	tangent = ON_3dVector(0.0, 1.0, 0.0);
+    else if (fabs(radial.z) < fabs(radial.x) &&
+	fabs(radial.z) < fabs(radial.y))
+	tangent = ON_3dVector(0.0, 0.0, 1.0);
+    return tangent;
+}
+
+
+static bool
+trace_mixed_span_junction(struct soltab *solid, struct rt_i *rtip,
+    struct resource *resource, const ON_3dVector &radial)
+{
+    return trace_oblique_span_grid(solid, rtip, resource, radial,
+	span_junction_tangent_axis(radial), SPAN_JUNCTION_JUNCTION_BOXES,
+	SPAN_JUNCTION_ADJACENT_BOXES);
+}
+
+
+static bool
+trace_singleton_span_junction(struct soltab *solid, struct rt_i *rtip,
+    struct resource *resource, const ON_3dVector &radial)
+{
+    const ON_3dVector tangent(1.0, 1.0, 1.0);
+    return trace_oblique_span_grid(solid, rtip, resource, radial, tangent,
+	SPAN_JUNCTION_JUNCTION_BOXES,
+	SPAN_JUNCTION_SINGLETON_BOXES);
+}
+
+
+static bool
+trace_adjacent_singleton_span(struct soltab *solid, struct rt_i *rtip,
+    struct resource *resource, ON_Brep *brep, int face_index)
+{
+    ON_3dVector boundary_radial;
+    if (!span_junction_boundary_radial(brep, face_index,
+	SPAN_JUNCTION_BOUNDARY_TRANSVERSE_FRACTION, boundary_radial))
+	return false;
+    return trace_oblique_span_grid(solid, rtip, resource,
+	boundary_radial, span_junction_tangent_axis(boundary_radial),
+	SPAN_JUNCTION_ADJACENT_BOXES, SPAN_JUNCTION_SINGLETON_BOXES);
+}
+
+
+int
+main(int argc, char **argv)
+{
+    bu_setprogname(argv[0]);
+    int face_index = SPAN_JUNCTION_DEFAULT_FACE_INDEX;
+    if (argc == 2) {
+	char *end = NULL;
+	errno = 0;
+	const long parsed_index = strtol(argv[1], &end, 10);
+	if (errno || !end || end == argv[1] || *end ||
+		parsed_index < 0 || parsed_index > INT_MAX)
+	    return 1;
+	face_index = (int)parsed_index;
+    } else if (argc != 1) {
+	return 1;
+    }
+    struct rt_i *rtip = rt_dirbuild_inmem(NULL, 0, NULL, 0);
+    if (!rtip) {
+	std::printf("FAIL: in-memory raytrace setup\n");
+	return 1;
+    }
+    rtip->rti_tol.magic = BN_TOL_MAGIC;
+    rtip->rti_tol.dist = 0.0005;
+    rtip->rti_tol.dist_sq = rtip->rti_tol.dist * rtip->rti_tol.dist;
+    rtip->rti_tol.perp = 1.0e-6;
+    rtip->rti_tol.para = 1.0 - rtip->rti_tol.perp;
+    if (!_rt_brep_set_surface_tree_depth(rtip,
+	SPAN_JUNCTION_SURFACE_TREE_DEPTH)) {
+	std::printf("FAIL: BREP surface tree setup\n");
+	rt_i_destroy(rtip);
+	return 1;
+    }
+
+    struct resource resource = {};
+    rt_init_resource(&resource, 0, rtip);
+    ON_Brep *brep = make_span_junction_sphere();
+    if (!brep || !brep->IsSolid()) {
+	std::printf("FAIL: span-junction sphere construction\n");
+	delete brep;
+	free_rtip(rtip, &resource);
+	return 1;
+    }
+    struct soltab *solid = prepare_brep_solid(rtip, brep);
+    if (!solid) {
+	std::printf("FAIL: span-junction sphere preparation\n");
+	delete brep;
+	free_rtip(rtip, &resource);
+	return 1;
+    }
+
+    ON_3dVector radial;
+    if (!span_junction_radial(brep, face_index, radial)) {
+	free_solid(solid);
+	free_rtip(rtip, &resource);
+	return 1;
+    }
+
+    point_t forward_origin;
+    vect_t forward_direction;
+    point_t reverse_origin;
+    vect_t reverse_direction;
+    VSET(forward_origin, SPAN_JUNCTION_ORIGIN_RADIUS * radial.x,
+	SPAN_JUNCTION_ORIGIN_RADIUS * radial.y,
+	SPAN_JUNCTION_ORIGIN_RADIUS * radial.z);
+    VSET(forward_direction, -radial.x, -radial.y, -radial.z);
+    VSET(reverse_origin, -SPAN_JUNCTION_ORIGIN_RADIUS * radial.x,
+	-SPAN_JUNCTION_ORIGIN_RADIUS * radial.y,
+	-SPAN_JUNCTION_ORIGIN_RADIUS * radial.z);
+    VSET(reverse_direction, radial.x, radial.y, radial.z);
+    const double distance_tolerance = std::max(1.0e-9,
+	(double)rtip->rti_tol.dist);
+
+    struct rt_brep_shot_trace forward_trace = {};
+    const int forward_hits = shoot_trace(solid, rtip, &resource,
+	forward_origin, forward_direction, forward_trace);
+    const bool forward_valid = trace_has_expected_events(forward_trace,
+	forward_hits, distance_tolerance);
+    if (!forward_valid) {
+	std::printf("FAIL: forward span-junction stream\n");
+	report_trace(forward_trace);
+    }
+
+    struct rt_brep_shot_trace reverse_trace = {};
+    const int reverse_hits = shoot_trace(solid, rtip, &resource,
+	reverse_origin, reverse_direction, reverse_trace);
+    const bool reverse_valid = trace_has_expected_events(reverse_trace,
+	reverse_hits, distance_tolerance);
+    if (!reverse_valid) {
+	std::printf("FAIL: reverse span-junction stream\n");
+	report_trace(reverse_trace);
+    }
+
+    const bool mixed_valid = face_index != SPAN_JUNCTION_DEFAULT_FACE_INDEX ||
+	trace_mixed_span_junction(solid, rtip, &resource, radial);
+    const bool singleton_valid =
+	face_index != SPAN_JUNCTION_DEFAULT_FACE_INDEX ||
+	trace_singleton_span_junction(solid, rtip, &resource, radial);
+    const bool adjacent_singleton_valid =
+	face_index != SPAN_JUNCTION_DEFAULT_FACE_INDEX ||
+	trace_adjacent_singleton_span(solid, rtip, &resource, brep, face_index);
+
+    free_solid(solid);
+    free_rtip(rtip, &resource);
+    return forward_valid && reverse_valid && mixed_valid && singleton_valid &&
+	adjacent_singleton_valid ? 0 : 1;
+}
+
+
+/*
+ * Local Variables:
+ * tab-width: 8
+ * mode: C++
+ * c-basic-offset: 4
+ * indent-tabs-mode: t
+ * c-file-style: "stroustrup"
+ * End:
+ * ex: shiftwidth=4 tabstop=8
+ */

src/librt/tree.c

@@ -28,6 +28,7 @@
 
 #include "bu/color.h"
 #include "bu/parallel.h"
+#include "bu/snooze.h"
 #include "vmath.h"
 #include "bn.h"
 #include "rt/db_attr.h"
@@ -246,11 +247,18 @@ _rt_gettree_region_end(struct db_tree_state *tsp, const struct db_full_path *pat
 	db_mater_color_region(tsp->ts_dbip, rp);
     bu_avs_free(&avs);
 
+    /* d_uses also indexes solid instances and must use the same directory-
+     * hash semaphore as _rt_find_identical_solid(). */
+    {
+	int hash = db_dirhash(dp->d_namep);
+	ACQUIRE_SEMAPHORE_TREE(hash);
+	rp->reg_instnum = dp->d_uses++;
+	RELEASE_SEMAPHORE_TREE(hash);
+    }
+
     /* enter critical section */
     bu_semaphore_acquire(RT_SEM_RESULTS);
 
-    rp->reg_instnum = dp->d_uses++;
-
     /*
      * Add the region to the linked list of regions.
      * Positions in the region bit vector are established at this time.
@@ -346,6 +354,7 @@ _rt_find_identical_solid(const matp_t mat, struct directory *dp, struct rt_i *rt
 
     hash = db_dirhash(dp->d_namep);
 
+retry:
     /* Enter the appropriate dual critical-section */
     ACQUIRE_SEMAPHORE_TREE(hash);
 
@@ -386,6 +395,16 @@ _rt_find_identical_solid(const matp_t mat, struct directory *dp, struct rt_i *rt
 	     */
 	    if (stp->st_rtip != rtip) continue;
 
+	    /* A matching soltab is linked before its potentially expensive
+	     * primitive prep begins.  Do not publish or reference its mutable
+	     * prep result until the creating worker has completed under this
+	     * same hash semaphore. */
+	    if (!stp->st_prep_complete) {
+		RELEASE_SEMAPHORE_TREE(hash);
+		bu_snooze(1000);
+		goto retry;
+	    }
+
 	    /*
 	     * stp now points to re-referenced solid.  stp->st_id is
 	     * non-zero, indicating pre-existing solid.
@@ -558,6 +577,7 @@ _rt_gettree_leaf(struct db_tree_state *tsp, const struct db_full_path *pathp, st
 	hash = db_dirhash(dp->d_namep);
 	ACQUIRE_SEMAPHORE_TREE(hash);
 	stp->st_aradius = -1;
+	stp->st_prep_complete = 1;
 	stp->st_uses--;
 	RELEASE_SEMAPHORE_TREE(hash);
 	return TREE_NULL;		/* BAD */
@@ -619,6 +639,15 @@ _rt_gettree_leaf(struct db_tree_state *tsp, const struct db_full_path *pathp, st
 	bu_vls_free(&str);
     }
 
+    /* Publish every primitive-specific pointer, bound, and path write before
+     * another reference is allowed to consume this soltab. */
+    {
+	int hash = db_dirhash(dp->d_namep);
+	ACQUIRE_SEMAPHORE_TREE(hash);
+	stp->st_prep_complete = 1;
+	RELEASE_SEMAPHORE_TREE(hash);
+    }
+
 found_it:
     BU_GET(curtree, union tree);
     RT_TREE_INIT(curtree);

src/proc-db/CMakeLists.txt

@@ -10,7 +10,7 @@ set(PROC_DB_WDB_NMG_LIBS "libwdb;libnmg;${M_LIBRARY}")
 set(PROC_DB_WDB_OPENNURBS_LIBS "libwdb;OPENNURBS::OPENNURBS;${M_LIBRARY}")
 
 brlcad_addexec(bottest bottest.c "${PROC_DB_WDB_LIBS}" NO_INSTALL)
-brlcad_addexec(brep_cobb brep_cobb.cpp "${PROC_DB_WDB_OPENNURBS_LIBS}" NO_STRICT NO_INSTALL)
+brlcad_addexec(brep_cobb brep_cobb.cpp "libbrep;${PROC_DB_WDB_OPENNURBS_LIBS}" NO_STRICT NO_INSTALL)
 brlcad_addexec(brep_cube "brep_cube.cpp;twistedcube.cpp" "${PROC_DB_WDB_OPENNURBS_LIBS}" NO_STRICT NO_INSTALL)
 target_include_directories(brep_cube BEFORE PRIVATE ${CMAKE_CURRENT_SOURCE_DIR})
 brlcad_addexec(brep_invalid "brep_invalid.cpp" "${PROC_DB_WDB_OPENNURBS_LIBS}" NO_STRICT NO_INSTALL)

src/proc-db/brep_cobb.cpp

@@ -1,293 +1,35 @@
-/*                   B R E P _ C O B B . C P P
+/*                    B R E P _ C O B B . C P P
  * BRL-CAD
  *
  * Copyright (c) 2013-2026 United States Government as represented by
  * the U.S. Army Research Laboratory.
  *
- * Based off of code from Ayam:
- *
- * This software is copyrighted (c) 1998-2012 by
- * Randolf Schultz (randolf.schultz@gmail.com).
- * All rights reserved.
- *
- * The author hereby grants permission to use, copy, modify, distribute,
- * and license this software and its documentation for any purpose, provided
- * that existing copyright notices are retained in all copies and that this
- * notice is included verbatim in any distributions. No written agreement,
- * license, or royalty fee is required for any of the authorized uses.
- * Modifications to this software may be copyrighted by their authors
- * and need not follow the licensing terms described here, provided that
- * the new terms are clearly indicated on the first page of each file where
- * they apply.
- *
- * IN NO EVENT SHALL THE AUTHOR OR DISTRIBUTORS BE LIABLE TO ANY PARTY
- * FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES
- * ARISING OUT OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY
- * DERIVATIVES THEREOF, EVEN IF THE AUTHOR HAS BEEN ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- *
- * THE AUTHOR AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES,
- * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY,
- * FITNESS FOR A PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE
- * IS PROVIDED ON AN "AS IS" BASIS, AND THE AUTHOR AND DISTRIBUTORS HAVE
- * NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR
- * MODIFICATIONS.
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * version 2.1 as published by the Free Software Foundation.
  */
 /** @file proc-db/brep_cobb.cpp
  *
  * Creates a Cobb NURBS sphere with cube topology, per
  *
- * J.E. Cobb, “Tiling the Sphere with Rational Bezier Patches,”
- *    Tech. Report TR UUCS-88-009, Computer Science Dept.,
- *    Univ. of Utah, Salt Lake City, 1988.
+ * J. E. Cobb, "Tiling the Sphere with Rational Bezier Patches,"
+ * Tech. Report TR UUCS-88-009, University of Utah, 1988.
  */
 
 #include "common.h"
 
 #include "bio.h"
+#include "brep.h"
 #include "bu/app.h"
 #include "wdb.h"
 
 
-ON_BezierSurface *
-ON_CobbSphereFace(double rotation_x, double rotation_z)
-{
-    ON_3dVector rotation_x_axis(1, 0, 0);
-    ON_3dVector rotation_z_axis(0, 0, 1);
-    ON_3dVector rotation_center(0, 0, 0);
-    double rot_x = rotation_x * ON_PI/180.0;
-    double rot_z = rotation_z * ON_PI/180.0;
-
-    double t = sqrt(3.0);
-    double d = sqrt(2.0);
-    ON_4dPoint cv;
-    ON_BezierSurface *p1 = new ON_BezierSurface(3, true, 5, 5);
-    p1->ReserveCVCapacity(100);
-
-    cv.x = 4.0*(1.0-t);
-    cv.y = 4.0*(1.0-t);
-    cv.z = 4.0*(1.0-t);
-    cv.w = 4.0*(3.0-t);
-    p1->SetCV(0,0,cv);
-
-
-    cv.x = -d;
-    cv.y = d * (t - 4.0);
-    cv.z = d * (t - 4.0);
-    cv.w = d * (3.0*t - 2.0);
-    p1->SetCV(0,1,cv);
-
-    cv.x = 0.0;
-    cv.y = 4.0*(1.0-2.0*t)/3.0;
-    cv.z = 4.0*(1.0-2.0*t)/3.0;
-    cv.w = 4.0*(5.0-t)/3.0;
-    p1->SetCV(0,2,cv);
-
-    cv.x = d;
-    cv.y = d * (t - 4.0);
-    cv.z = d * (t - 4.0);
-    cv.w = d * (3*t - 2.0);
-    p1->SetCV(0,3,cv);
-
-    cv.x = 4.0*(t - 1.0);
-    cv.y = 4.0*(1.0-t);
-    cv.z = 4.0*(1.0-t);
-    cv.w = 4.0*(3.0-t);
-    p1->SetCV(0,4,cv);
-
-    cv.x = d*(t - 4.0);
-    cv.y = -d;
-    cv.z = d*(t-4.0);
-    cv.w = d*(3.0*t-2.0);
-    p1->SetCV(1,0,cv);
-
-    cv.x = (2.0-3.0*t)/2.0;
-    cv.y = (2.0-3.0*t)/2.0;
-    cv.z = -(t+6.0)/2.0;
-    cv.w = (t+6.0)/2.0;
-    p1->SetCV(1,1,cv);
-
-    cv.x = 0.0;
-    cv.y = d*(2.0*t-7.0)/3.0;
-    cv.z = -5.0*sqrt(6.0)/3.0;
-    cv.w = d*(t+6.0)/3.0;
-    p1->SetCV(1,2,cv);
-
-    cv.x = (3.0*t-2.0)/2.0;
-    cv.y = (2.0-3.0*t)/2.0;
-    cv.z = -(t+6.0)/2.0;
-    cv.w = (t+6.0)/2.0;
-    p1->SetCV(1,3,cv);
-
-    cv.x = d*(4.0-t);
-    cv.y = -d;
-    cv.z = d*(t-4.0);
-    cv.w = d*(3.0*t-2.0);
-    p1->SetCV(1,4,cv);
-
-    cv.x = 4.0*(1.0-2*t)/3.0;
-    cv.y = 0.0;
-    cv.z = 4.0*(1.0-2*t)/3.0;
-    cv.w = 4.0*(5.0-t)/3.0;
-    p1->SetCV(2,0,cv);
-
-    cv.x = d*(2.0*t-7.0)/3.0;
-    cv.y = 0.0;
-    cv.z = -5.0*sqrt(6.0)/3.0;
-    cv.w = d*(t+6.0)/3.0;
-    p1->SetCV(2,1,cv);
-
-    cv.x = 0.0;
-    cv.y = 0.0;
-    cv.z = 4.0*(t-5.0)/3.0;
-    cv.w = 4.0*(5.0*t-1.0)/9.0;
-    p1->SetCV(2,2,cv);
-
-    cv.x = -d*(2*t-7.0)/3.0;
-    cv.y = 0.0;
-    cv.z = -5.0*sqrt(6.0)/3.0;
-    cv.w = d*(t+6.0)/3.0;
-    p1->SetCV(2,3,cv);
-
-    cv.x = -4.0*(1.0-2.0*t)/3.0;
-    cv.y = 0.0;
-    cv.z = 4.0*(1.0-2.0*t)/3.0;
-    cv.w = 4.0*(5.0-t)/3.0;
-    p1->SetCV(2,4,cv);
-
-    cv.x = d*(t-4.0);
-    cv.y = d;
-    cv.z = d*(t-4.0);
-    cv.w = d*(3.0*t-2);
-    p1->SetCV(3,0,cv);
-
-    cv.x = (2.0-3.0*t)/2.0;
-    cv.y = -(2.0-3.0*t)/2.0;
-    cv.z = -(t+6.0)/2.0;
-    cv.w = (t+6.0)/2.0;
-    p1->SetCV(3,1,cv);
-
-    cv.x = 0.0;
-    cv.y = -d*(2.0*t-7.0)/3.0;
-    cv.z = -5.0*sqrt(6.0)/3.0;
-    cv.w = d*(t+6.0)/3.0;
-    p1->SetCV(3,2,cv);
-
-    cv.x = (3.0*t-2.0)/2.0;
-    cv.y = -(2.0-3.0*t)/2.0;
-    cv.z = -(t+6.0)/2.0;
-    cv.w = (t+6.0)/2.0;
-    p1->SetCV(3,3,cv);
-
-    cv.x = d*(4.0-t);
-    cv.y = d;
-    cv.z = d*(t-4.0);
-    cv.w = d*(3.0*t-2.0);
-    p1->SetCV(3,4,cv);
-
-    cv.x = 4.0*(1.0-t);
-    cv.y = -4.0*(1.0-t);
-    cv.z = 4.0*(1.0-t);
-    cv.w = 4.0*(3.0-t);
-    p1->SetCV(4,0,cv);
-
-    cv.x = -d;
-    cv.y = -d*(t-4.0);
-    cv.z = d*(t-4.0);
-    cv.w = d*(3.0*t-2.0);
-    p1->SetCV(4,1,cv);
-
-    cv.x = 0.0;
-    cv.y = -4.0*(1.0-2.0*t)/3.0;
-    cv.z = 4.0*(1.0-2.0*t)/3.0;
-    cv.w = 4.0*(5.0-t)/3.0;
-    p1->SetCV(4,2,cv);
-
-    cv.x = d;
-    cv.y = -d*(t-4.0);
-    cv.z = d*(t-4.0);
-    cv.w = d*(3.0*t-2.0);
-    p1->SetCV(4,3,cv);
-
-    cv.x = 4.0*(t-1.0);
-    cv.y = -4.0*(1.0-t);
-    cv.z = 4.0*(1.0-t);
-    cv.w = 4.0*(3.0-t);
-    p1->SetCV(4,4,cv);
-
-    p1->Rotate(rot_x, rotation_x_axis, rotation_center);
-    p1->Rotate(rot_z, rotation_z_axis, rotation_center);
-
-    return p1;
-}
-
-/* TODO - Need to find a more compact, efficient way to
- * do this - shouldn't need 24 3d curves... */
-ON_Brep *
-Cobb_Sphere(double UNUSED(radius), ON_3dPoint *UNUSED(origin))
-{
-    ON_Brep *b = ON_Brep::New();
-
-    // Patch 1 of 6
-    ON_BezierSurface *b1 = ON_CobbSphereFace(0, 0);
-    ON_NurbsSurface *p1_nurb = ON_NurbsSurface::New();
-    b1->GetNurbForm(*p1_nurb);
-    b->NewFace(*p1_nurb);
-    delete b1;
-
-    // Patch 2 of 6
-    ON_BezierSurface *b2 = ON_CobbSphereFace(90, 0);
-    ON_NurbsSurface *p2_nurb = ON_NurbsSurface::New();
-    b2->GetNurbForm(*p2_nurb);
-    b->NewFace(*p2_nurb);
-    delete b2;
-
-    // Patch 3 of 6
-    ON_BezierSurface *b3 = ON_CobbSphereFace(180, 0);
-    ON_NurbsSurface *p3_nurb = ON_NurbsSurface::New();
-    b3->GetNurbForm(*p3_nurb);
-    b->NewFace(*p3_nurb);
-    delete b3;
-
-    // Patch 4 of 6
-    ON_BezierSurface *b4 = ON_CobbSphereFace(270, 0);
-    ON_NurbsSurface *p4_nurb = ON_NurbsSurface::New();
-    b4->GetNurbForm(*p4_nurb);
-    b->NewFace(*p4_nurb);
-    delete b4;
-
-    // Patch 5 of 6
-    ON_BezierSurface *b5 = ON_CobbSphereFace(90, 90);
-    ON_NurbsSurface *p5_nurb = ON_NurbsSurface::New();
-    b5->GetNurbForm(*p5_nurb);
-    b->NewFace(*p5_nurb);
-    delete b5;
-
-    // Patch 6 of 6
-    ON_BezierSurface *b6 = ON_CobbSphereFace(90, 270);
-    ON_NurbsSurface *p6_nurb = ON_NurbsSurface::New();
-    b6->GetNurbForm(*p6_nurb);
-    b->NewFace(*p6_nurb);
-    delete b6;
-
-
-    b->Standardize();
-    b->Compact();
-
-    return b;
-}
-
 int
-main(int argc, char** argv)
+main(int argc, char **argv)
 {
-    struct rt_wdb* outfp;
-    ON_Brep* brep;
-    ON_TextLog error_log;
-    ON_3dPoint origin(0.0,0.0,0.0);
-    const char* id_name = "B-Rep Cobb Sphere";
-    const char* geom_name = "cobb.s";
-    const char* db_name = "brep_cobb.g";
-
+    const char *id_name = "B-Rep Cobb Sphere";
+    const char *geom_name = "cobb.s";
+    const char *db_name = "brep_cobb.g";
     bu_setprogname(argv[0]);
 
     if (argc > 1) {
@@ -303,10 +45,8 @@ main(int argc, char** argv)
     }
 
     ON::Begin();
-
-    /* export brep to file */
     bu_log("Writing a Cobb unit sphere b-rep to %s...\n", db_name);
-    outfp = wdb_fopen(db_name);
+    struct rt_wdb *outfp = wdb_fopen(db_name);
     if (!outfp) {
 	bu_log("ERROR: unable to open %s for writing\n", db_name);
 	ON::End();
@@ -314,23 +54,25 @@ main(int argc, char** argv)
     }
     mk_id(outfp, id_name);
 
-
-    brep = Cobb_Sphere(1, &origin);
+    const ON_3dPoint origin(0.0, 0.0, 0.0);
+    ON_Brep *brep = ON_Brep_CobbSphereUnsewn(1.0, origin);
+    if (!brep) {
+	bu_log("ERROR: unable to construct the Cobb sphere\n");
+	db_close(outfp->dbip);
+	ON::End();
+	return 1;
+    }
     mk_brep(outfp, geom_name, (void *)brep);
 
-    //mk_comb1(outfp, "cube.r", geom_name, 1);
     unsigned char rgb[] = {50, 255, 50};
     mk_region1(outfp, "cobb.r", geom_name, "plastic", "", rgb);
 
     db_close(outfp->dbip);
     delete brep;
-
     ON::End();
-
     return 0;
 }
 
-
 /*
  * Local Variables:
  * tab-width: 8