saved work added here as a patch in trunk - this branch is no longer needed (svn branch delete)

svn:revision:63044
svn:branch:surfacetree
svn:account:starseeker

src/libbrep/CMakeLists.txt

@@ -27,6 +27,9 @@ set(LIBBREP_SOURCES
   px_event.cpp
   PullbackCurve.cpp
   ssx_event.cpp
+  libbrep_curvetree.cpp
+  libbrep_surfacetree.cpp
+  libbrep_util.cpp
   )
 
 set(libbrep_ignored_files
@@ -34,6 +37,9 @@ set(libbrep_ignored_files
   opennurbs_fit.h
   opennurbs_fit.cpp
   PullbackCurve.h
+  libbrep_curvetree.h
+  libbrep_surfacetree.h
+  libbrep_util.h
   )
 CMAKEFILES(${libbrep_ignored_files})
 

src/libbrep/libbrep_curvetree.cpp

@@ -0,0 +1,622 @@
+/*          L I B B R E P _ C U R V E T R E E . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2013 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.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file libbrep_curvetree.cpp
+ *
+ * Brief description
+ *
+ */
+
+#include "libbrep_curvetree.h"
+
+#include "vmath.h"
+#include <iostream>
+
+#include <vector>
+#include <map>
+#include <set>
+#include <queue>
+#include <limits>
+
+// Return 0 if there are no tangent points in the interval, 1
+// if there is a horizontal tangent, two if there is a vertical
+// tangent, 3 if a normal split is in order.
+int ON_Curve_Has_Tangent(const ON_Curve* curve, double min, double max) {
+
+    bool tanx1, tanx2, x_changed;
+    bool tany1, tany2, y_changed;
+    bool slopex, slopey;
+    double xdelta, ydelta;
+    ON_3dVector tangent1, tangent2;
+    ON_3dPoint p1, p2;
+    ON_Interval t(min, max);
+
+    tangent1 = curve->TangentAt(t[0]);
+    tangent2 = curve->TangentAt(t[1]);
+
+    tanx1 = (tangent1[X] < 0.0);
+    tanx2 = (tangent2[X] < 0.0);
+    tany1 = (tangent1[Y] < 0.0);
+    tany2 = (tangent2[Y] < 0.0);
+
+    x_changed =(tanx1 != tanx2);
+    y_changed =(tany1 != tany2);
+
+    if (x_changed && y_changed) return 3; //horz & vert
+    if (x_changed) return 1;//need to get vertical tangent
+    if (y_changed) return 2;//need to find horizontal tangent
+
+    p1 = curve->PointAt(t[0]);
+    p2 = curve->PointAt(t[1]);
+
+    xdelta = (p2[X] - p1[X]);
+    slopex = (xdelta < 0.0);
+    ydelta = (p2[Y] - p1[Y]);
+    slopey = (ydelta < 0.0);
+
+    // If we have no slope change
+    // in x or y, we have a tangent line
+    if (NEAR_ZERO(xdelta, TOL) || NEAR_ZERO(ydelta, TOL)) return 0;
+
+    if ((slopex != tanx1) || (slopey != tany1)) return 3;
+
+    return 0;
+}
+
+int ON_Curve_HV_Split(const ON_Curve *curve, double min, double max, double *mid) {
+    // Don't split if min ~= max - that's infinite loop territory
+    if (NEAR_ZERO(max - min, ON_ZERO_TOLERANCE)) return 0;
+    // Go HV tangent hunting
+    int tan_result = ON_Curve_Has_Tangent(curve, min, max);
+    if (tan_result == 3) return 1;
+    if (tan_result == 0) return 0;
+    if (tan_result == 1 || tan_result == 2) {
+	// If we have a tangent point at the min or max, that's OK
+	ON_3dVector hv_tangent = curve->TangentAt(max);
+	if (NEAR_ZERO(hv_tangent.x, TOL2) || NEAR_ZERO(hv_tangent.y, TOL2)) tan_result = 0;
+	hv_tangent = curve->TangentAt(min);
+	if (NEAR_ZERO(hv_tangent.x, TOL2) || NEAR_ZERO(hv_tangent.y, TOL2)) tan_result = 0;
+	// The tangent point is in a bad spot - find it and subdivide on it
+	if(tan_result) {
+	    bool tanmin;
+	    (tan_result == 1) ? (tanmin = (hv_tangent[X] < 0.0)) : (tanmin = (hv_tangent[Y] < 0.0));
+	    double hv_mid;
+	    double hv_min = min;
+	    double hv_max = max;
+	    // find the tangent point using binary search
+	    while (fabs(hv_min - hv_max) > TOL2) {
+		hv_mid = (hv_max + hv_min)/2.0;
+		hv_tangent = curve->TangentAt(hv_mid);
+		if ((tan_result == 1 && NEAR_ZERO(hv_tangent[X], TOL2))
+			||(tan_result == 2 && NEAR_ZERO(hv_tangent[Y], TOL2))) {
+		    (*mid) = hv_mid;
+		    return 1;
+		}
+		if ((tan_result == 1 && ((hv_tangent[X] < 0.0) == tanmin)) ||
+			(tan_result == 2 && ((hv_tangent[Y] < 0.0) == tanmin) )) {
+		    hv_min = hv_mid;
+		} else {
+		    hv_max = hv_mid;
+		}
+	    }
+	    return 0;
+	} else {
+          // The tangent was at the min and/or max, we don't need to split on account
+          // of this case
+          return 0;
+        }
+    }
+    // Shouldn't get here - something went drastically wrong if we did
+    return -1;
+}
+
+double ON_Curve_Solve_For_V(const ON_Curve *curve, double u, double min, double max, double tol) {
+
+    ON_3dVector Tan_start, Tan_end;
+    ON_3dPoint p;
+    double guess, dT;
+    double Ta = min;
+    double Tb = max;
+    ON_3dPoint A = curve->PointAt(Ta);
+    ON_3dPoint B = curve->PointAt(Tb);
+    double dU = fabs(A.x - B.x);
+    if (dU <= tol) {
+       if (A.y <= B.y) {return A.y;} else {return B.y;};
+    }
+    Tan_start = curve->TangentAt(Ta);
+    Tan_end = curve->TangentAt(Tb);
+    dT = Tb - Ta;
+
+    /* Use quick binary subdivision until derivatives at end points in 'u' are within 5 percent */
+    while (!(fabs(dU) <= tol) && !(fabs(dT) <= tol)) {
+        guess = Ta + dT/2;
+        p = curve->PointAt(guess);
+
+	/* if we hit 'u' exactly, done deal */
+	if (fabs(p.x-u) <= SMALL_FASTF) return p.y;
+
+        if (p.x > u) {
+            /* v is behind us, back up the end */
+            Tb = guess;
+            B = p;
+            Tan_end = curve->TangentAt(Tb);
+        } else {
+            /* v is in front, move start forward */
+            Ta = guess;
+            A = p;
+            Tan_start = curve->TangentAt(Ta);
+        }
+        dT = Tb - Ta;
+        dU = B.x - A.x;
+    }
+
+    dU = fabs(B.x - A.x);
+    if (dU <= tol) {  //vertical
+       if (A.y <= B.y) {return A.y;} else {return B.y;};
+    }
+
+    guess = Ta + (u - A.x) * dT/dU;
+    p = curve->PointAt(guess);
+
+    int cnt=0;
+    while ((cnt < 1000) && (!(fabs(p.x-u) <= tol))) {
+        if (p.x < u) {
+            Ta = guess;
+            A = p;
+        } else {
+            Tb = guess;
+            B = p;
+        }
+        dU = fabs(B.x - A.x);
+        if (dU <= tol) {  //vertical
+	    if (A.y <= B.y) {return A.y;} else {return B.y;};
+	}
+
+        dT = Tb - Ta;
+        guess =Ta + (u - A[X]) * dT/dU;
+        p = curve->PointAt(guess);
+        cnt++;
+    }
+    if (cnt > 999) {
+        std::cout << "ON_Curve_Solve_For_V(): estimate of 'v' given a trim curve and 'u' did not converge within iteration bound(" << cnt << ")\n";
+    }
+    return p.y;
+}
+
+int ON_Curve_Relative_Size(const ON_Curve *curve, const ON_Surface *srf, double tmin, double tmax, double factor) {
+    int retval = 0;
+    ON_Interval u = srf->Domain(0);
+    ON_Interval v = srf->Domain(1);
+    ON_3dPoint a(u[0], v[0], 0.0);
+    ON_3dPoint b(u[1], v[1], 0.0);
+    double ab = a.DistanceTo(b);
+    double cd = curve->PointAt(tmin).DistanceTo(curve->PointAt(tmax));
+    if (cd > factor*ab) retval = 1;
+    return retval;
+}
+
+long int ON_CurveTree::New_Trim_Node(int parent, int trim_index, int tmin, int tmax) {
+    // Parent node of new node
+    nodes.push_back(parent);
+    // Curve parameter information
+    nodes.push_back(tmin);
+    nodes.push_back(tmax);
+    // As yet, the node has no children
+    nodes.push_back(0);
+    nodes.push_back(0);
+    // knots/HV tangents info node
+    nodes.push_back(0);
+    // Populate uv_bbox_values
+    nodes.push_back(uv_bbox_values.size());
+    for (int i = 0; i < 2; ++i) uv_bbox_values.push_back(std::numeric_limits<double>::infinity());
+    for (int i = 2; i < 4; ++i) uv_bbox_values.push_back(-1 * std::numeric_limits<double>::infinity());
+    // Trim index in brep
+    nodes.push_back(trim_index);
+
+    return nodes.size() - TRIM_NODE_STEP;
+}
+
+unsigned int ON_CurveTree::Depth(long int node_id) {
+    if (node_id < first_trim_node) return 0;
+    int depth = 1;
+    int parent = nodes[node_id];
+    while (parent != 0 && parent >= first_trim_node) {
+        depth++;
+        parent = nodes[parent];
+    }
+    return depth;
+}
+
+int ON_Curve_Flat(const ON_Curve *curve, double min, double max, double tol){
+    int retval = 1;
+    // If min ~= max, call it flat
+    if (NEAR_ZERO(max - min, ON_ZERO_TOLERANCE)) return retval;
+    ON_3dVector tangent_start = curve->TangentAt(min);
+    ON_3dVector tangent_end = curve->TangentAt(max);
+    if (fabs(tangent_start * tangent_end) < tol) retval = 0;
+    return retval;
+}
+
+void ON_CurveTree::Subdivide_Trim_Node(long int first_node_id, const ON_BrepTrim *trim) {
+
+    // Initialize
+    ON_3dPoint pmin, pmax;
+    std::set<long int> leaf_nodes;
+    const ON_Curve *curve = trim->TrimCurveOf();
+    const ON_Surface *srf = trim->SurfaceOf();
+    double c_min, c_max;
+    int knotcnt = curve->SpanCount();
+    double *knots = (double *)malloc(sizeof(double) * (knotcnt + 1));
+    (void)curve->GetSpanVector(knots);
+    (void)curve->GetDomain(&c_min, &c_max);
+
+    // Populate the root node t_coords
+    t_coords.push_back(c_min);
+    t_coords.push_back(c_max);
+    nodes.at(first_node_id + 1) = t_coords.size() - 2;
+    nodes.at(first_node_id + 2) = t_coords.size() - 1;
+
+    // Use a queue to process nodes
+    std::queue<long int> trim_node_queue;
+
+    // Prime the queue with the root node
+    trim_node_queue.push(first_node_id);
+
+    // Begin the core of the tree building process. Pop a node
+    // off the queue, test it to see if it must be split further,
+    // and proceed accordingly.
+    while (!trim_node_queue.empty()) {
+	int split = 0;
+	long int node = trim_node_queue.front();
+        trim_node_queue.pop();
+        double node_tmin = t_coords.at(nodes.at(node+1));
+        double node_tmax = t_coords.at(nodes.at(node+2));
+        //std::cout << "tmin: " << node_tmin << ", tmax: " << node_tmax << "\n";
+	ON_Interval node_ival(node_tmin, node_tmax);
+
+	// Determine the depth of this node
+	int node_depth = this->Depth(node);
+	//std::cout << "Node depth: " << node_depth << "\n";
+
+        // Use a variable for the midpoint, since it may be overridden
+        // if we have knots or tangent points to split on.
+	double t_mid = node_ival.Mid();
+
+        // If we have knots in the interval, we need to split
+	split += ON_Interval_Find_Split_Knot(&node_ival, knotcnt, knots, &t_mid);
+        if (split) nodes.at(node+5) = 1;
+	if (split) {
+	    //std::cout << "Knot split:" << node << "," << t_mid << "\n";
+	    //if(NEAR_ZERO(t_mid-node_tmin, TOL2) || NEAR_ZERO(t_mid-node_tmax, TOL2)) std::cout << "Arrgh! Knots\n";
+	}
+        // If we aren't already splitting, check for horizontal and vertical tangents
+	if (!split) {
+	    split += ON_Curve_HV_Split(curve, node_tmin, node_tmax, &t_mid);
+	    if (split) nodes.at(node+5) = 1;
+	    if (split) {
+		//std::cout << "HV split:" << node << "," << t_mid << "\n";
+		//if(NEAR_ZERO(t_mid-node_tmin, TOL2) || NEAR_ZERO(t_mid-node_tmax, TOL2)) std::cout << "Arrgh! HV\n";
+	    }
+	}
+
+        // If we aren't already splitting, check flatness
+        if (!split) {
+	    split += !ON_Curve_Flat(curve, node_tmin, node_tmax, BREP_CURVE_FLATNESS);
+	    if (split) {
+		//std::cout << "Flatness:" << node << "\n";
+		//if(NEAR_ZERO(t_mid-node_tmin, TOL2) || NEAR_ZERO(t_mid-node_tmax, TOL2)) std::cout << "Arrgh! Flat\n";
+	    }
+	}
+
+	// We want subdivided curves to be small relative to their
+	// parent surface.
+	if (!split) {
+	    split += ON_Curve_Relative_Size(curve, srf, node_tmin, node_tmax, BREP_TRIM_SUB_FACTOR);
+	    if (split) {
+		// std::cout << "Relative:" << node << "\n";
+		//if(NEAR_ZERO(t_mid-node_tmin, TOL2) || NEAR_ZERO(t_mid-node_tmax, TOL2)) std::cout << "Arrgh! Relative\n";
+	    }
+	}
+
+        // If we DO need to split, proceed
+        if (split) {
+	    long int left_child, right_child;
+	    // Push the new midpoint onto tcoords
+	    t_coords.push_back(t_mid);
+	    left_child = this->New_Trim_Node(node, nodes.at(node+7), nodes.at(node+1), t_coords.size() - 1);
+	    nodes.at(node+3) = left_child;
+	    trim_node_queue.push(left_child);
+	    right_child = this->New_Trim_Node(node, nodes.at(node+7), t_coords.size() - 1, nodes.at(node+2));
+	    nodes.at(node+4) = right_child;
+	    trim_node_queue.push(right_child);
+        } else {
+            double pmin_u, pmin_v, pmax_u, pmax_v;
+            leaf_nodes.insert(node);
+            pmin = curve->PointAt(node_tmin);
+            pmax = curve->PointAt(node_tmax);
+            (pmin.x > pmax.x) ? (pmin_u = pmax.x) : (pmin_u = pmin.x);
+            (pmax.x > pmin.x) ? (pmax_u = pmax.x) : (pmax_u = pmin.x);
+            (pmin.y > pmax.y) ? (pmin_v = pmax.y) : (pmin_v = pmin.y);
+            (pmax.y > pmin.y) ? (pmax_v = pmax.y) : (pmax_v = pmin.y);
+            uv_bbox_values.at(nodes.at(node+6)) = pmin_u;
+            uv_bbox_values.at(nodes.at(node+6)+1) = pmin_v;
+            uv_bbox_values.at(nodes.at(node+6)+2) = pmax_u;
+            uv_bbox_values.at(nodes.at(node+6)+3) = pmax_v;
+	    //std::cout << "node " << node << " bbox:" << pmin_u << "," << pmin_v << "," << pmax_u << "," << pmax_v << "\n";
+        }
+
+    }
+    // Clean up
+    free(knots);
+
+    // Sync bboxes up to the root trim node
+    std::set<long int> node_set;
+    std::set<long int>::iterator p_it;
+    std::set<long int> *parent_nodes, *current_nodes, *temp;
+    double *p_umin, *p_vmin, *p_umax, *p_vmax;
+    double *c_umin, *c_vmin, *c_umax, *c_vmax;
+    current_nodes = &leaf_nodes;
+    parent_nodes = &node_set;
+    while (!current_nodes->empty()) {
+	for (p_it = current_nodes->begin(); p_it != current_nodes->end(); ++p_it) {
+	    if ((*p_it) > first_node_id) {
+		long int parent_id = nodes.at((*p_it));
+		p_umin = &(uv_bbox_values[nodes.at(parent_id+6)]);
+		p_vmin = &(uv_bbox_values[nodes.at(parent_id+6)+1]);
+		p_umax = &(uv_bbox_values[nodes.at(parent_id+6)+2]);
+		p_vmax = &(uv_bbox_values[nodes.at(parent_id+6)+3]);
+		//std::cout << "pnode " << parent_id << " bbox:" << *p_umin << "," << *p_vmin << "," << *p_umax << "," << *p_vmax << "\n";
+		c_umin = &(uv_bbox_values[nodes.at((*p_it)+6)]);
+		c_vmin = &(uv_bbox_values[nodes.at((*p_it)+6)+1]);
+		c_umax = &(uv_bbox_values[nodes.at((*p_it)+6)+2]);
+		c_vmax = &(uv_bbox_values[nodes.at((*p_it)+6)+3]);
+		//std::cout << "cnode " << (*p_it) << " bbox:" << *c_umin << "," << *c_vmin << "," << *c_umax << "," << *c_vmax << "\n";
+		if (*c_umin < *p_umin) (*p_umin) = (*c_umin);
+		if (*c_vmin < *p_vmin) (*p_vmin) = (*c_vmin);
+		if (*c_umax > *p_umax) (*p_umax) = (*c_umax);
+		if (*c_vmax > *p_vmax) (*p_vmax) = (*c_vmax);
+		if (parent_id > first_node_id) parent_nodes->insert(parent_id);
+		//std::cout << "pnode u" << parent_id << " bbox:" << *p_umin << "," << *p_vmin << "," << *p_umax << "," << *p_vmax << "\n";
+	    }
+	}
+        current_nodes->clear();
+        temp = parent_nodes;
+        parent_nodes = current_nodes;
+        current_nodes = temp;
+    }
+}
+
+// When testing a curve for trimming contributions, the first stage is to check
+// the curve tree node to see whether it contains a segment from below the knot
+// and HV tangent breakdowns of the curve.  If it does not, we need to walk down
+// the tree until we have the set of nodes that represent corresponding sub-segments
+// that are below the knot and HV tangent breakdowns.
+//
+// Once we have an appropriate segment set, do the pnpoly intersection checks.
+// For anything where the point is not in the segment bounding box, the test is
+// simple - if the point is within the bounding box of the node, walk down the
+// tree until either a determination can be made or the point is found to be
+// within a leaf node.
+//
+// If the point is within a leaf node, the more expensive search against the
+// actual curve to determine 2D ray intersection status is necessary.
+bool ON_CurveTree::CurveTrim(double u, double v, long int node) {
+    bool c = false;
+    std::set<long int> test_nodes;
+    std::set<long int>::iterator tn_it;
+    std::queue<long int> q1, q2;
+    std::queue<long int> *current_queue, *next_queue, *tmp;
+    // Get the node or nodes we need into test_nodes
+    q1.push(node);
+    current_queue = &q1;
+    next_queue= &q2;
+    while (!current_queue->empty()) {
+	while (!current_queue->empty()) {
+	    long int current_node = current_queue->front();
+	    current_queue->pop();
+	    if (nodes.at(current_node+5) != 1) {
+		test_nodes.insert(current_node);
+	    } else {
+		next_queue->push(nodes.at(current_node+3));
+		next_queue->push(nodes.at(current_node+4));
+	    }
+	}
+	tmp = current_queue;
+	current_queue = next_queue;
+	next_queue = tmp;
+    }
+    // Check node bounding boxes vs uv point
+    double *c_umin, *c_vmin, *c_umax, *c_vmax;
+    for(tn_it = test_nodes.begin(); tn_it != test_nodes.end(); ++tn_it) {
+	c_umin = &uv_bbox_values[nodes.at((*tn_it)+6)];
+	c_vmin = &uv_bbox_values[nodes.at((*tn_it)+6)+1];
+	c_umax = &uv_bbox_values[nodes.at((*tn_it)+6)+2];
+	c_vmax = &uv_bbox_values[nodes.at((*tn_it)+6)+3];
+        // Are we inside the u interval of the bbox?
+        if ((u <= *c_umax) && (u > *c_umin)) {
+           // Are we below the bbox?  if so, a v+ ray is guaranteed to intersect
+           if (v <= *c_vmin) {
+              c = !c;
+           } else {
+               // We're not below it - are we inside it?  If not, no intersection possible,
+               // otherwise, need closer look
+	       if (v < *c_vmax) {
+                   // Decend the tree until either we're outside
+                   // the bboxes and can make a determination, or
+                   // we find the leaf node containing the uv point
+                   // and do the full curve test.
+		   current_queue->push((*tn_it));
+		   while (!current_queue->empty()) {
+		       while (!current_queue->empty()) {
+			   long int current_node = current_queue->front();
+			   current_queue->pop();
+			   c_umin = &uv_bbox_values[nodes.at(current_node+6)];
+			   c_vmin = &uv_bbox_values[nodes.at(current_node+6)+1];
+			   c_umax = &uv_bbox_values[nodes.at(current_node+6)+2];
+			   c_vmax = &uv_bbox_values[nodes.at(current_node+6)+3];
+			   if ((u <= *c_umax) && (u >= *c_umin)) {
+			       if (v <= *c_vmin) {
+				   c = !c;
+			       } else {
+				   if (v < *c_vmax) {
+                                       // If we are at a leaf, it's time to do a
+                                       // close check.  Otherwise, queue the children
+				       if (nodes.at(current_node + 3)) {
+					   next_queue->push(nodes.at(current_node + 3));
+					   next_queue->push(nodes.at(current_node + 4));
+				       } else {
+					   // We're close to a leaf
+					   const ON_Curve *curve = this->ctree_face->Brep()->m_C2[(int)nodes.at(current_node + 7)];
+					   double v_curve = ON_Curve_Solve_For_V(curve, u, nodes.at(current_node + 1), nodes.at(current_node + 2), TOL3);
+					   if (v <= v_curve) c = !c;
+				       }
+				   }
+			       }
+			   }
+		       }
+		       tmp = current_queue;
+		       current_queue = next_queue;
+		       next_queue = tmp;
+		   }
+	       }
+           }
+        }
+    }
+    return c;
+}
+
+// Remember - unlike other
+// trimming loops, points inside the outer loop are untrimmed and
+// points outside it ARE trimmed.  Flip test results accordingly.
+bool ON_CurveTree::IsTrimmed(double u, double v, std::vector<long int> *subset, bool pnpoly_state) {
+    //std::cout << "u :" << u << ", v :" << v << "\n";
+    int last_loop = 0;
+    std::set<long int> loops_to_check;
+    long int node_pos = 0;
+    // First stage - identify which trimming loop(s) we care about, based on loop UV bounding boxes
+    while (last_loop != 1) {
+	double *p_umin = &uv_bbox_values[nodes.at(node_pos+2)];
+	double *p_vmin = &uv_bbox_values[nodes.at(node_pos+2)+1];
+	double *p_umax = &uv_bbox_values[nodes.at(node_pos+2)+2];
+	double *p_vmax = &uv_bbox_values[nodes.at(node_pos+2)+3];
+        //std::cout << "loop " << node_pos << " bbox:" << (*p_umin) << "," << (*p_vmin) << "," << (*p_umax) << "," << (*p_vmax) << "\n";
+        // Is uv point inside bbox?
+        if (!( u < (*p_umin) || v < (*p_vmin) || u > (*p_umax) || v > (*p_vmax) )) loops_to_check.insert(node_pos);
+	// Jump the starting position to the next loop.
+        node_pos = nodes.at(node_pos);
+	if (node_pos == 0) last_loop = 1;
+    }
+    std::set<long int>::iterator l_it;
+    if (loops_to_check.size() > 1) {
+	std::cout << "Loops possibly containing UV pt (" << u << "," << v << "):";
+	for(l_it = loops_to_check.begin(); l_it != loops_to_check.end(); ++l_it) {
+	    std::cout << " " << (*l_it) << " ";
+	}
+    std::cout << "\n";
+    for(l_it = loops_to_check.begin(); l_it != loops_to_check.end(); ++l_it) {
+       long int current_node_pos = (*l_it) + 3;
+       long int final_node_pos = (*l_it) + nodes.at((*l_it)+1) + 3;
+       while (current_node_pos <= final_node_pos) {
+           //std::cout << "loop " << (*l_it) << ", trim " << nodes.at(current_node_pos) << "\n";
+           current_node_pos++;
+       }
+    }
+    std::cout << "\n";
+}
+}
+
+ON_CurveTree::ON_CurveTree(const ON_BrepFace* face)
+{
+    // Save the face pointer associated with this tree, in case we
+    // need to access data from the face
+    ctree_face = face;
+
+    // We can find out in advance how many node elements we will have due to
+    // trimming loops - populate those elements up front.
+    long int elements = 0;
+    for (int loop_index = 0; loop_index < face->LoopCount(); loop_index++) {
+	elements += face->Loop(loop_index)->TrimCount();
+    }
+    elements += (face->LoopCount()) * NON_TRIM_LOOP_STEP;
+    // Import to initialize to zero - used later as stopping criteria
+    nodes.assign(elements, 0);
+
+    // Any node elements after those already in place will pertain to trim
+    // nodes - to be able to easily distinguish when an index is referring
+    // to a loop or a trim, record the index which marks the cross-over point
+    // in the vector.  The test is then a simple < or >= check.
+    first_trim_node = nodes.size();
+
+    // Now that we have the setup, add actual information and handle the trims
+    long int node_pos = 0;
+    for (int loop_index = 0; loop_index < face->LoopCount(); loop_index++) {
+	ON_BrepLoop *loop = face->Loop(loop_index);
+	nodes.at(node_pos + 1) = loop->TrimCount();
+	// Populate uv_bbox_values
+	nodes.at(node_pos + 2) = uv_bbox_values.size();
+	for (int i = 0; i < 2; ++i) uv_bbox_values.push_back(std::numeric_limits<double>::infinity());
+	for (int i = 2; i < 4; ++i) uv_bbox_values.push_back(-1 * std::numeric_limits<double>::infinity());
+
+	for (int trim_index = 0; trim_index < loop->TrimCount(); trim_index++) {
+            int curve_index = loop->Trim(trim_index)->TrimCurveIndexOf();
+            long int current_trim_node = node_pos+NON_TRIM_LOOP_STEP+trim_index;
+            // Add "leaf" trim node from loop
+            long int trim_root = this->New_Trim_Node(node_pos, curve_index, 0, 0);
+	    nodes.at(current_trim_node) = trim_root;
+            // Build KDtree below "leaf" trim
+            this->Subdivide_Trim_Node(trim_root, face->Loop(loop_index)->Trim(trim_index));
+	}
+        // Now that we've handled the trims, build the loop bbox
+	double *p_umin = &uv_bbox_values[nodes.at(node_pos+2)];
+	double *p_vmin = &uv_bbox_values[nodes.at(node_pos+2)+1];
+	double *p_umax = &uv_bbox_values[nodes.at(node_pos+2)+2];
+	double *p_vmax = &uv_bbox_values[nodes.at(node_pos+2)+3];
+	for (int trim_index = 0; trim_index < loop->TrimCount(); trim_index++) {
+	    long int current_trim_node = nodes.at(node_pos+NON_TRIM_LOOP_STEP+trim_index);
+	    double *c_umin = &uv_bbox_values[nodes.at(current_trim_node+6)];
+	    double *c_vmin = &uv_bbox_values[nodes.at(current_trim_node+6)+1];
+	    double *c_umax = &uv_bbox_values[nodes.at(current_trim_node+6)+2];
+	    double *c_vmax = &uv_bbox_values[nodes.at(current_trim_node+6)+3];
+	    if (*c_umin < *p_umin) (*p_umin) = (*c_umin);
+            if (*c_vmin < *p_vmin) (*p_vmin) = (*c_vmin);
+            if (*c_umax > *p_umax) (*p_umax) = (*c_umax);
+            if (*c_vmax > *p_vmax) (*p_vmax) = (*c_vmax);
+	}
+        std::cout << "loop " << node_pos << " bbox:" << (*p_umin) << "," << (*p_vmin) << "," << (*p_umax) << "," << (*p_vmax) << "\n";
+
+        // Identify the size of this loop node and store it in the first entry
+        // to enable "walking" of the loops
+        // Jump the starting position to the next loop.
+        if (loop_index != face->LoopCount() - 1)
+	    nodes.at(node_pos) = node_pos + loop->TrimCount() + NON_TRIM_LOOP_STEP + 1;
+	node_pos = nodes.at(node_pos);
+    }
+}
+
+ON_CurveTree::~ON_CurveTree() {
+}
+
+// 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/libbrep/libbrep_curvetree.h

@@ -0,0 +1,180 @@
+/*            L I B B R E P _ C U R V E T R E E . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2013 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.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file libbrep_curvetree.h
+ *
+ * Brief description
+ *
+ */
+
+#if !defined(__LIBBREP_CURVETREE)
+#define __LIBBREP_CURVETREE
+
+#include "libbrep_util.h"
+
+#include <vector>
+#include <map>
+#include <set>
+#include <queue>
+
+#ifndef NURBS_EXPORT
+#  if defined(NURBS_DLL_EXPORTS) && defined(NURBS_DLL_IMPORTS)
+#    error "Only NURBS_DLL_EXPORTS or NURBS_DLL_IMPORTS can be defined, not both."
+#  elif defined(NURBS_DLL_EXPORTS)
+#    define NURBS_EXPORT __declspec(dllexport)
+#  elif defined(NURBS_DLL_IMPORTS)
+#    define NURBS_EXPORT __declspec(dllimport)
+#  else
+#    define NURBS_EXPORT
+#  endif
+#endif
+
+/* trim curve point sampling count for isLinear() check and possibly growing bounding box*/
+#define BREP_BB_CRV_PNT_CNT 10
+#define BREP_CURVE_FLATNESS 0.95
+
+/* subdivision size factors */
+#define BREP_TRIM_SUB_FACTOR 1
+
+/// arbitrary calculation tolerance (need to try VDIVIDE_TOL or VUNITIZE_TOL to tighten the bounds)
+#ifndef TOL
+#  define TOL 0.000001
+#endif
+//
+///// another arbitrary calculation tolerance (need to try VDIVIDE_TOL or VUNITIZE_TOL to tighten the bounds)
+#ifndef TOL2
+#  define TOL2 0.00001
+#endif
+
+/// another arbitrary calculation tolerance (used in getting a value of v from a curve given a u value)
+#ifndef TOL3
+#  define TOL3 0.0000001
+#endif
+
+////////////////////////////////////////////////////////////////
+//
+//   ON_CurveTree
+//
+// An ON_BrepFace has 1 or more trimming loops, which
+// are in turn made up of trimming curves.
+//
+// There are actually two "tiers" to a curve tree - the loop level,
+// which manages which loops need to be addressed in particular
+// queries, and the trim level which actually uses tree structures
+// to break curves into sufficiently small sub-segments.
+//
+// The position of each vector element in a "loop node" is as follows:
+// index + 0  : holds the index of the next loop node.  Unlike
+//              trim nodes, loop node element counts aren't fixed. The
+//              final loop node will have value of zero in this slot.
+// index + 1  : number of trim curves in this loop
+// index + 2  : index into the uv_bbox_values vector identifying the umin value associated
+//              with this node. The next three values in the array after umin will be vmin, umax
+//              and vmax.
+// index + 3 + n  : index of the nth trim root node associated with this loop (counting from index 0)
+
+#define NON_TRIM_LOOP_STEP 3
+
+
+// The position of each vector element in a "trim node" is as follows:
+//
+// index + 0  : holds the index of the parent node of this node. If this is a trim node under
+//              a loop node, this index identifies the index of the parent loop node.
+// index + 1  : index into the t_coords vector of doubles identifying tmin for this node.
+// index + 2  : index into the t_coords vector of doubles identifying tmax for this node.
+// index + 3  : index of left child node of this node.
+// index + 4  : index of right child node of this node.
+// index + 5  : information on whether or not this node has problematic points such as
+//              knots or HV tangents - 1 means problem points are present, 0 is clean,
+//              2 means clean but curve segment is part of outer loop.
+// index + 6  : index into the uv_bbox_values vector identifying the umin value associated
+//              with this node. The next three values in the array after umin will be vmin, umax
+//              and vmax.  Building these values up from the leaf nodes ensures proper bounding.
+// index + 7  : index of associated curve in BRep's m_C2 array
+
+#define TRIM_NODE_STEP 8
+
+class ON_CLASS ON_CurveTree
+{
+    public:
+
+	ON_CurveTree(const ON_BrepFace*);
+	~ON_CurveTree();
+
+	unsigned int Depth(long int);  // Depth of a given node in the tree.
+
+        // Set a bounding box to the overall UV bounding box of the curve tree
+        void BoundingBox(ON_BoundingBox *);
+
+        bool CurveTrim(double u, double v, long int node);
+
+	// One of the most important queries that can be made of a Curve tree is whether
+	// a given point is trimmed or untrimmed by the loops of the tree.  This test
+	// uses a modified PNPOLY test.  The two modifications are 1) a polygon is dynamically formed
+	// for PNPOLY using the corsest polygon that can be constructed - once a point is
+	// outside a bounding box for a curve segment, the linear approximation for that
+	// curve segment is sufficient for the test - and 2), when a point is too close
+	// to a curve segment to be sufficient for the ray intersection test, use a test
+	// that actually queries the curve.
+	bool IsTrimmed(
+		double u,
+		double v,
+		std::vector<long int> *subset = NULL,
+		bool pnpoly_state = false
+		);
+
+	// We also need to be able to tell whether a UV rectangular region is fully trimmed,
+	// fully present, or partially trimmed.  The patch may want to pass in a subset of
+	// the curves it is interested in, get back a subset of the curves that overlap its
+	// domain, or both.
+	int IsTrimmed(
+		ON_Interval *u_dom,
+		ON_Interval *v_dom,
+		std::vector<long int> *included = NULL,
+		std::vector<long int> *subset = NULL,
+		bool pnpoly_state = false
+		);
+
+        // We store the face pointer so we have access to face info at
+        // need - for example, the outer loop pointer is face->OuterLoop()
+        // and the outer loop index is face->OuterLoop()->m_loop_index
+        const ON_BrepFace *ctree_face;
+
+    private:
+        long int first_trim_node;
+        long int New_Trim_Node(int parent, int curve_index, int tmin, int tmax);
+        void Subdivide_Trim_Node(long int node_id, const ON_BrepTrim *);
+
+        std::vector<long int> nodes;
+        std::vector<double> t_coords;
+        std::vector<double> uv_bbox_values;
+
+        ON_BoundingBox ctree_bbox;
+};
+
+#endif
+
+// 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/libbrep/libbrep_surfacetree.cpp

@@ -0,0 +1,439 @@
+/*        L I B B R E P _ S U R F A C E T R E E . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2013 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.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file libbrep_surfacetree.cpp
+ *
+ * Brief description
+ *
+ */
+
+#include <vector>
+#include <map>
+#include <set>
+#include <queue>
+#include <iostream>
+
+#include "libbrep_surfacetree.h"
+
+long int ON_SurfaceTree::New_Node(int parent, int umin, int umax, int vmin, int vmax) {
+    // Parent node of new node
+    nodes.push_back(parent);
+    // UV information
+    nodes.push_back(umin);
+    nodes.push_back(umax);
+    nodes.push_back(vmin);
+    nodes.push_back(vmax);
+    // As yet, the node has no children.
+    for (int i = 0; i < 4; ++i) nodes.push_back(0);
+    ON_BoundingBox surf_bbox;
+    bboxes.Append(surf_bbox);
+    nodes.push_back(bboxes.Count() - 1);
+    nodes.push_back(0); // We haven't tested for trimming curves yet.
+    nodes.push_back(0); // No temporary nodes as yet.
+    return nodes.size() - SURF_NODE_STEP;
+}
+
+// Split a node to produce new nodes in both U and V directions
+void ON_SurfaceTree::Node_Split_UV(double u_mid, double v_mid, std::pair<int, std::vector<int> *> *node, std::queue<std::vector<int> *> *frame_arrays, std::queue<std::pair<int, std::vector<int> *> > *node_queue) {
+    std::vector<int> *parent_frame_index = node->second;
+
+   // Four new frame index arrays will be needed
+    std::vector<int> *q1f = NULL;
+    std::vector<int> *q2f = NULL;
+    std::vector<int> *q3f = NULL;
+    std::vector<int> *q4f = NULL;
+
+    // Add the midpoint U and V coordiates, which will be the
+    // only UV values not already present in the double array
+    uv_coords.push_back(u_mid);
+    int umid_id = uv_coords.size() - 1;
+    uv_coords.push_back(v_mid);
+    int vmid_id = uv_coords.size() - 1;
+
+    // Create child nodes.  Reuse frame arrays if any are available to reuse
+
+    // Create the SW node
+    long int q1 = this->New_Node(node->first, nodes[node->first+1], umid_id, nodes[node->first+3], vmid_id);
+    nodes[node->first+5] = q1;
+    if (!frame_arrays->empty()) {q1f = frame_arrays->front(); frame_arrays->pop();} else {q1f = new std::vector<int>;}
+    q1f->assign(13, -1);
+    q1f->at(0) = parent_frame_index->at(0);
+    q1f->at(1) = parent_frame_index->at(9);
+    q1f->at(2) = parent_frame_index->at(4);
+    q1f->at(3) = parent_frame_index->at(10);
+    q1f->at(4) = parent_frame_index->at(5);
+    node_queue->push(std::make_pair(q1,q1f));
+
+    // Create the SE node
+    long int q2 = this->New_Node(node->first, umid_id, nodes[node->first+2], nodes[node->first+3], vmid_id);
+    nodes[node->first+6] = q2;
+    if (!frame_arrays->empty()) {q2f = frame_arrays->front(); frame_arrays->pop();} else {q2f = new std::vector<int>;}
+    q2f->assign(13, -1);
+    q2f->at(0) = parent_frame_index->at(9);
+    q2f->at(1) = parent_frame_index->at(1);
+    q2f->at(2) = parent_frame_index->at(12);
+    q2f->at(3) = parent_frame_index->at(4);
+    q2f->at(4) = parent_frame_index->at(7);
+    node_queue->push(std::make_pair(q2,q2f));
+
+    // Create the NE node
+    long int q3 = this->New_Node(node->first, umid_id, nodes[node->first+2], vmid_id, nodes[node->first+4]);
+    nodes[node->first+7] = q3;
+    if (!frame_arrays->empty()) {q3f = frame_arrays->front(); frame_arrays->pop();} else {q3f = new std::vector<int>;}
+    q3f->assign(13, -1);
+    q3f->at(0) = parent_frame_index->at(4);
+    q3f->at(1) = parent_frame_index->at(12);
+    q3f->at(2) = parent_frame_index->at(2);
+    q3f->at(3) = parent_frame_index->at(11);
+    q3f->at(4) = parent_frame_index->at(8);
+    node_queue->push(std::make_pair(q3,q3f));
+
+    // Create the NW node
+    long int q4 = this->New_Node(node->first, nodes[node->first+1], umid_id, vmid_id, nodes[node->first+4]);
+    nodes[node->first+8] = q4;
+    if (!frame_arrays->empty()) {q4f = frame_arrays->front(); frame_arrays->pop();} else {q4f = new std::vector<int>;}
+    q4f->assign(13, -1);
+    q4f->at(0) = parent_frame_index->at(10);
+    q4f->at(1) = parent_frame_index->at(4);
+    q4f->at(2) = parent_frame_index->at(11);
+    q4f->at(3) = parent_frame_index->at(3);
+    q4f->at(4) = parent_frame_index->at(6);
+    node_queue->push(std::make_pair(q4,q4f));
+}
+
+// Create two child nodes, splitting in the U direction
+void ON_SurfaceTree::Node_Split_U(double u_mid, std::pair<int, std::vector<int> *> *node, std::queue<std::vector<int> *> *frame_arrays, std::queue<std::pair<int, std::vector<int> *> > *node_queue) {
+    std::vector<int> *parent_frame_index = node->second;
+
+    // Two new frame index arrays will be needed
+    std::vector<int> *q1f = NULL;
+    std::vector<int> *q2f = NULL;
+
+    // Add the midpoint U coordiate, which will be the
+    // only UV value not already present in the double array
+    uv_coords.push_back(u_mid);
+    int umid_id = uv_coords.size() - 1;
+
+    // Create child nodes.  Reuse frame arrays if any are available to reuse
+
+    // Create the W node
+    long int q1 = this->New_Node(node->first, nodes[node->first+1], umid_id, nodes[node->first+3], nodes[node->first+4]);
+    nodes[node->first+5] = q1;
+    if (!frame_arrays->empty()) {q1f = frame_arrays->front(); frame_arrays->pop();} else {q1f = new std::vector<int>;}
+    q1f->assign(13, -1);
+    q1f->at(0) = parent_frame_index->at(0);
+    q1f->at(1) = parent_frame_index->at(9);
+    q1f->at(2) = parent_frame_index->at(11);
+    q1f->at(3) = parent_frame_index->at(3);
+    q1f->at(10) = parent_frame_index->at(10);
+    q1f->at(12) = parent_frame_index->at(4);
+    node_queue->push(std::make_pair(q1,q1f));
+
+    // Create the E node
+    long int q2 = this->New_Node(node->first, umid_id, nodes[node->first+2], nodes[node->first+3], nodes[node->first+4]);
+    nodes[node->first+6] = q2;
+    if (!frame_arrays->empty()) {q2f = frame_arrays->front(); frame_arrays->pop();} else {q2f = new std::vector<int>;}
+    q2f->assign(13, -1);
+    q2f->at(0) = parent_frame_index->at(9);
+    q2f->at(1) = parent_frame_index->at(1);
+    q2f->at(2) = parent_frame_index->at(2);
+    q2f->at(3) = parent_frame_index->at(11);
+    q2f->at(10) = parent_frame_index->at(4);
+    q2f->at(12) = parent_frame_index->at(12);
+    node_queue->push(std::make_pair(q2,q2f));
+
+}
+
+// Create two child nodes, splitting in the V direction
+void ON_SurfaceTree::Node_Split_V(double v_mid, std::pair<int, std::vector<int> *> *node, std::queue<std::vector<int> *> *frame_arrays, std::queue<std::pair<int, std::vector<int> *> > *node_queue) {
+    std::vector<int> *parent_frame_index = node->second;
+
+    // Two new frame index arrays will be needed
+    std::vector<int> *q1f = NULL;
+    std::vector<int> *q2f = NULL;
+
+    // Add the midpoint V coordiate, which will be the
+    // only UV value not already present in the double array
+    uv_coords.push_back(v_mid);
+    int vmid_id = uv_coords.size() - 1;
+
+    // Create child nodes.  Reuse frame arrays if any are available to reuse
+
+    // Create the S node
+    long int q1 = this->New_Node(node->first, nodes[node->first+1], nodes[node->first+2], nodes[node->first+3], vmid_id);
+    nodes[node->first+7] = q1;
+    if (!frame_arrays->empty()) {q1f = frame_arrays->front(); frame_arrays->pop();} else {q1f = new std::vector<int>;}
+    q1f->assign(13, -1);
+    q1f->at(0) = parent_frame_index->at(0);
+    q1f->at(1) = parent_frame_index->at(1);
+    q1f->at(2) = parent_frame_index->at(12);
+    q1f->at(3) = parent_frame_index->at(10);
+    q1f->at(9) = parent_frame_index->at(9);
+    q1f->at(11) = parent_frame_index->at(4);
+    node_queue->push(std::make_pair(q1,q1f));
+
+    // Create the N node
+    long int q2 = this->New_Node(node->first, nodes[node->first+1], nodes[node->first+2], vmid_id, nodes[node->first+4]);
+    nodes[node->first+8] = q2;
+    if (!frame_arrays->empty()) {q2f = frame_arrays->front(); frame_arrays->pop();} else {q2f = new std::vector<int>;}
+    q2f->assign(13, -1);
+    q2f->at(0) = parent_frame_index->at(10);
+    q2f->at(1) = parent_frame_index->at(12);
+    q2f->at(2) = parent_frame_index->at(2);
+    q2f->at(3) = parent_frame_index->at(3);
+    q2f->at(9) = parent_frame_index->at(4);
+    q2f->at(11) = parent_frame_index->at(11);
+    node_queue->push(std::make_pair(q2,q2f));
+
+}
+
+unsigned int ON_SurfaceTree::Depth(long int node_id) {
+    if (node_id == 1) return 0; //root node
+    int depth = 1;
+    int parent = nodes[node_id];
+    while (parent != 1) {
+	depth++;
+	parent = nodes[parent];
+    }
+    return depth;
+}
+
+ON_SurfaceTree::ON_SurfaceTree(const ON_BrepFace* face) {
+    int uv = 0;
+    int first_node_id = 0;
+    ON_SimpleArray<ON_Plane> frames;
+
+    // We need a curve tree
+    m_ctree = new ON_CurveTree(face);
+
+    // Initialize
+    std::set<long int> leaf_nodes;
+    srftree_face = face;
+    const ON_Surface *srf = face->SurfaceOf();
+    int u_knotcnt = srf->SpanCount(0);
+    int v_knotcnt = srf->SpanCount(1);
+    double *u_knots = (double *)malloc(sizeof(double) *(u_knotcnt + 1));
+    double *v_knots = (double *)malloc(sizeof(double) *(v_knotcnt + 1));
+    (void)srf->GetSpanVector(0, u_knots);
+    (void)srf->GetSpanVector(1, v_knots);
+    ON_Interval u_domain = srf->Domain(0);
+    ON_Interval v_domain = srf->Domain(1);
+
+    nodes.push_back(0);
+    temp_nodes.push_back(1);
+    // When we start splitting surfaces later, we need
+    // some scratch surfaces to avoid lots of malloc
+    // activity.
+    ON_Surface *t1 = NULL;
+    ON_Surface *t2 = NULL;
+    ON_Surface *t3 = NULL;
+    ON_Surface *t4 = NULL;
+    ON_Surface *sub_surface = NULL;
+    ON_Surface_Create_Scratch_Surfaces(&t1, &t2, &t3, &t4);
+
+    // Initialize the tree with the root node.
+    uv_coords.push_back(u_domain.Min());
+    uv_coords.push_back(u_domain.Max());
+    uv_coords.push_back(v_domain.Min());
+    uv_coords.push_back(v_domain.Max());
+    uv = uv_coords.size() - 1;
+    first_node_id = this->New_Node(1, uv - 3, uv - 2, uv - 1, uv);
+
+    // Each node has a set of frames, but those frames are
+    // often shared with parent and sibling nodes.  In order
+    // to take at least some advantage of this, we store
+    // all frames in an array and use vectors of indicies to
+    // point to a specific frame.  When a vector is no longer
+    // needed, we queue it up in an "available" queue for reuse.
+    std::queue<std::vector<int> *> frame_arrays;
+
+    // Initialize with the root frame vector
+    std::vector<int> *root_frame_index = new std::vector<int>;
+    root_frame_index->assign(13, -1);
+
+    // Processing of nodes is handled with a queue - each node,
+    // if it splits out into child nodes, queues up its child
+    // nodes in the work queue for further processing.  The
+    // queue holds only two pieces of information - a nodes position
+    // in the vector, and a pointer to the frame vector associated
+    // with this particular node.
+    std::queue<std::pair<int, std::vector<int> *> > node_queue;
+
+    // Prime the queue with the root node.
+    node_queue.push(std::make_pair(first_node_id,root_frame_index));
+
+    // Begin the core of the tree building process.  Pop a node
+    // off the queue, test it to see if it must be split further,
+    // and proceed accordingly.
+    while (!node_queue.empty()) {
+	int node_depth = 1;
+	int split = 0;
+	int u_split = 0;
+	int v_split = 0;
+	int trimmed = 0;
+	std::pair<int, std::vector<int> *> node = node_queue.front();
+	node_queue.pop();
+	int node_id = node.first;
+	std::vector<int> *frame_index = node.second;
+	ON_Interval u_local(uv_coords.at(nodes[node_id+1]), uv_coords.at(nodes[node_id+2]));
+	ON_Interval v_local(uv_coords.at(nodes[node_id+3]), uv_coords.at(nodes[node_id+4]));
+
+        // Determine the depth of this node
+        node_depth = this->Depth(node_id);
+        //std::cout << "Node depth: " << node_depth << "\n";
+
+	// Find out the trimming status of this node.  If it is
+	// fully trimmed out of the tree, there is no need to process it any further.
+	trimmed = ON_Surface_Patch_Trimmed(srf, m_ctree, &u_local, &v_local);
+
+	//if (!trimmed == 0) {
+	if (trimmed == 0 && node_depth < 9) {
+
+	    // Because the knots may override using Mid(), we have local variables to hold the
+	    // working values.
+	    double u_mid = u_local.Mid();
+	    double v_mid = v_local.Mid();
+            m_ctree->IsTrimmed(u_mid, v_mid);
+	    // Get the local sub-surface
+
+	    bool split_success = ON_Surface_SubSurface(srf, &sub_surface, &u_local, &v_local, &t1, &t2, &t3, &t4);
+
+	    // Test for knots - if we have knots, we must split regardless of flatness or dimensions.
+	    split += ON_Surface_Knots_Split(&u_local, &v_local, u_knotcnt, v_knotcnt, u_knots, v_knots, &u_mid, &v_mid);
+            //if (split) std::cout << "Knot split: " << u_mid << "," << v_mid << "\n";
+	    if (fabs(u_local.Mid() - u_mid) > ON_ZERO_TOLERANCE) u_split++;
+	    if (fabs(v_local.Mid() - v_mid) > ON_ZERO_TOLERANCE) v_split++;
+	    // If we're proceeding, go ahead and populate the frame
+	    if (ON_Populate_Frame_Array(srf, &u_local, &v_local, u_mid, v_mid, &frames, frame_index))
+		std::cout << "Warning - Initial frame array build failure!";
+
+
+	    // If we have one or more trimming curve nodes in the UV domain, but the collective U
+	    // domain (or V domain?  which?) of the curve nodes doesn't fully cover the U domain
+	    // of this node, divide on the node edge U that is closest to the center of the domain.
+	    // TODO - implement this - simplifies the amount of information we need to store for trimming
+	    // testing down the road, since we can make a patch-wide determination of the trimming
+	    // contribution of all curve nodes that don't actually intersect the surface patch.
+	    //
+	    // will also need to allow for both trimmed and untrimmed tree builds...
+
+
+	    // Pull the NURBS surface for this patch and test the
+	    // "unrolled" width and height - too distored, and we will need to split in just U
+	    // or just V to try and keep reasonably close to "square" patches.
+	    if (split_success)
+		split += ON_Surface_Width_vs_Height(sub_surface, 5.0, &u_split, &v_split);
+
+	    // Check the flatness and straightness of this particular patch
+	    if (!split) {
+		// test flat and straight
+		split += ON_Surface_Flat_Straight(&frames, frame_index, BREP_SURFACE_FLATNESS, BREP_SURFACE_STRAIGHTNESS);
+		if (split) {
+		    u_split += ON_Surface_Flat_U(&frames, frame_index, BREP_SURFACE_FLATNESS);
+		    v_split += ON_Surface_Flat_V(&frames, frame_index, BREP_SURFACE_FLATNESS);
+		}
+	    }
+
+	    // If we DO need to split, proceed
+	    if (split) {
+		// For efficiency, we calculate four additional frame here that are shared by the children
+		if (ON_Extend_Frame_Array(srf, &u_local, &v_local, u_mid, v_mid, &frames, frame_index))
+		    std::cout << "Warning - Frame array extension failure!";
+
+                // Generate the new nodes
+		if ((u_split && v_split) || (!u_split && !v_split)) {
+		    this->Node_Split_UV(u_mid, v_mid, &node, &frame_arrays, &node_queue);
+		} else {
+		    if (u_split && !v_split) {
+			this->Node_Split_U(u_mid, &node, &frame_arrays, &node_queue);
+		    }
+		    if (!u_split && v_split) {
+			this->Node_Split_V(v_mid, &node, &frame_arrays, &node_queue);
+		    }
+		}
+	    } else {
+                // We've got a leaf node - bounding box is presumed trustworthy
+                leaf_nodes.insert(node_id);
+                ON_BoundingBox *surf_bbox = bboxes.At((int)nodes[node_id+9]);
+                sub_surface->GetBBox(surf_bbox->m_min, surf_bbox->m_max);
+            }
+	}
+	// Whether this is a leaf or has produced children, we are done with its
+	// frame vector.  Add it to the frame vector queue for re-use.
+	node.second->clear();
+	frame_arrays.push(node.second);
+    }
+
+    // Now that we've got all our nodes, build up the bboxes
+    std::set<long int> node_set;
+    std::set<long int>::iterator p_it;
+    std::set<long int> *parent_nodes, *current_nodes, *temp;
+    current_nodes = &leaf_nodes;
+    parent_nodes = &node_set;
+    while (!current_nodes->empty()) {
+          for (p_it = current_nodes->begin(); p_it != current_nodes->end(); ++p_it) {
+              ON_BoundingBox *curr_bbox = bboxes.At((int)nodes[(*p_it)+9]);
+              ON_BoundingBox *parent_bbox = bboxes.At((int)nodes[nodes[(*p_it)]+9]);
+              parent_bbox->Set(curr_bbox->m_min, true);
+              parent_bbox->Set(curr_bbox->m_max, true);
+              if (nodes[(*p_it)] > 1) parent_nodes->insert(nodes[(*p_it)]);
+          }
+	  current_nodes->clear();
+	  temp = parent_nodes;
+	  parent_nodes = current_nodes;
+	  current_nodes = temp;
+    }
+
+    // Now that we're all done with sub-surfaces, delete our working surface memory
+    delete t1;
+    delete t2;
+    delete t3;
+    delete t4;
+    free(u_knots);
+    free(v_knots);
+    while (!frame_arrays.empty()) {
+	delete frame_arrays.front();
+	frame_arrays.pop();
+    }
+    frames.Destroy();
+}
+
+
+ON_SurfaceTree::~ON_SurfaceTree()
+{
+}
+
+/*
+TODO Add the fast bbox/ray intersection test from:
+http://www.cg.cs.tu-bs.de/media/publications/fast-rayaxis-aligned-bounding-box-overlap-tests-using-ray-slopes.pdf
+http://graphics.tu-bs.de/people/eisemann/code/RaySlopeIntersection.zip
+
+Will probably have use for this too (inside-outside polygon test, does the 2d raycasting bit):
+http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
+*/
+
+
+// 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/libbrep/libbrep_surfacetree.h

@@ -0,0 +1,200 @@
+/*          L I B B R E P _ S U R F A C E T R E E . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2013 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.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file libbrep_surfacetree.h
+ *
+ * Brief description
+ *
+ */
+
+#if !defined(LIBBREP_SURFACETREE_)
+#define LIBBREP_SURFACETREE_
+
+#include <vector>
+#include <map>
+#include <set>
+#include <queue>
+
+#include "libbrep_util.h"
+#include "libbrep_curvetree.h"
+
+
+/* Maximum per-surface BVH depth */
+#define BREP_MAX_FT_DEPTH 8
+#define BREP_MAX_LN_DEPTH 20
+#define SIGN(x) ((x) >= 0 ? 1 : -1)
+/* Surface flatness parameter, Abert says between 0.8-0.9 */
+#define BREP_SURFACE_FLATNESS 0.85
+#define BREP_SURFACE_STRAIGHTNESS 0.75
+/* Max newton iterations when finding closest point */
+#define BREP_MAX_FCP_ITERATIONS 50
+/* Root finding epsilon */
+#define BREP_FCP_ROOT_EPSILON 1e-5
+
+/* subdivision size factors */
+#define BREP_SURF_SUB_FACTOR 1
+
+/**
+ * The EDGE_MISS_TOLERANCE setting is critical in a couple of ways -
+ * too small and the allowed uncertainty region near edges will be
+ * smaller than the actual uncertainty needed for accurate solid
+ * raytracing, too large and trimming will not be adequate.  May need
+ * to adapt this to the scale of the model, perhaps using bounding box
+ * size to key off of.
+ */
+/* #define BREP_EDGE_MISS_TOLERANCE 5e-2 */
+#define BREP_EDGE_MISS_TOLERANCE 5e-3
+#define BREP_SAME_POINT_TOLERANCE 1e-6
+
+
+/// arbitrary calculation tolerance (need to try VDIVIDE_TOL or VUNITIZE_TOL to tighten the bounds)
+#ifndef TOL
+#  define TOL 0.000001
+#endif
+//
+///// another arbitrary calculation tolerance (need to try VDIVIDE_TOL or VUNITIZE_TOL to tighten the bounds)
+#ifndef TOL2
+#  define TOL2 0.00001
+#endif
+
+//////////////////////////////////////////////////
+
+//         Surface Tree Functionality
+
+//////////////////////////////////////////////////
+/*
+struct ON_SplitTest_Container {
+    int split = 0;
+    int u_split = 0;
+    int v_split = 0;
+    int trimmed = 0;
+    std::vector<int> *frame_indices;
+    ON_Interval u_local;
+    ON_Interval v_local;
+    ON_Surface *t1 = NULL;
+    ON_Surface *t2 = NULL;
+    ON_Surface *t3 = NULL;
+    ON_Surface *sub_surface = NULL;
+    ON_SimpleArray<ON_Plane> frames;
+}
+*/
+
+////////////////////////////////////////////////////////////////
+//
+//   ON_SurfaceTree
+//
+//  A "node" is a series of long integers stored in a vector,
+//  which point to either the index location of other nodes
+//  in the vector or the index into an asset array.
+//
+//  The first element (index 0) in the vector identifies the
+//  tree type.  The categories are as follows
+//
+//  0 - standard quad tree
+//  1 - temporary quad tree built using leaf nodes of a primary tree
+//
+//  The purpose of this mechanism is to allow the primary tree to
+//  contain references into a temporary tree, while still allowing the
+//  tree walking and other mechanisms to access the correct data. All
+//  functions doing array lookups should be supplied with an array
+//  type as a function parameter and use it to work with either nodes
+//  or temp_nodes to retrieve data.
+//
+//  The position of each vector element in a "node" is as follows:
+//
+//  index + 0  :  holds the index of the parent node of this node.
+//                The root node (1) is its own parent and is set to
+//                1.  The tree's root node (1) is never a valid child node.
+//  index + 1  :  Index into the vector of doubles identifying umin for this node.
+//  index + 2  :  Index into the vector of doubles identifying umax for this node.
+//  index + 3  :  Index into the vector of doubles identifying vmin for this node.
+//  index + 4  :  Index into the vector of doubles identifying vmax for this node.
+//                Specific values are retrieved with the offsets:
+//  index + 5  :  holds the index of the southwest child node of
+//                this node.  A value of zero indicates this node
+//                does not have this child node.
+//  index + 6  :  holds the index of the southeast child node of
+//                this node.  A value of zero indicates this node
+//                does not have this child node.
+//  index + 7  :  holds the index of the northeast child node of
+//                this node.  A value of zero indicates this node
+//                does not have this child node.
+//  index + 8  :  holds the index of the northwest child node of
+//                this node.  A value of zero indicates this node
+//                does not have this child node.
+//  index + 9  :  Index of the node's BoundingBox.
+//  index + 10 :  Trimming status of this node: 0 means untested,
+//                1 means fully trimmed away, 2 means no trimming curves
+//                intersect this node, and 3 means a test must
+//                be made to determine if a particular point within
+//                this patch is trimmed or untrimmed.  Matters only
+//                for leaf nodes - those with childeren will ask them.
+//  index + 11 :  The index of a temporary child node created
+//                using this node as a starting point.  The default
+//                value is zero, indicating no such node exists.
+//                The purpose of this mechanism is to allow a temporary
+//                deepening of a surface tree locally, and then allow
+//                the temporary resources to be cleared while leaving
+//                the primary tree intact.  If the index is non-zero
+//                it is a reference into the temp_nodes vector, rather
+//                than the primary vector.
+
+#define SURF_NODE_STEP 12
+
+class ON_CLASS  ON_SurfaceTree
+{
+    public:
+
+	ON_SurfaceTree(const ON_BrepFace*);
+	~ON_SurfaceTree();
+
+        unsigned int Depth(long int);  // Depth of a given node in the tree
+
+        // Store the face pointer, in case face information is needed
+        const ON_BrepFace *srftree_face;
+	//void LeafNodes(std::set<int_pair_t> *);
+    private:
+
+        long int New_Node(int, int, int, int, int);
+	void Node_Split_UV(double, double, std::pair<int, std::vector<int> *> *, std::queue<std::vector<int> *> *, std::queue<std::pair<int, std::vector<int> *> > *);
+	void Node_Split_U(double, std::pair<int, std::vector<int> *> *, std::queue<std::vector<int> *> *, std::queue<std::pair<int, std::vector<int> *> > *);
+	void Node_Split_V(double, std::pair<int, std::vector<int> *> *, std::queue<std::vector<int> *> *, std::queue<std::pair<int, std::vector<int> *> > *);
+
+	ON_CurveTree *m_ctree;
+	std::vector<long int> nodes;
+	std::vector<double> uv_coords;
+	ON_SimpleArray<ON_BoundingBox> bboxes;
+	// Contains for temporary data - these are intended
+	// to be cleared after completion of a task requiring
+	// local tree refinement
+	std::vector<long int> temp_nodes;
+	std::vector<double> temp_uv_coords;
+	ON_SimpleArray<ON_BoundingBox> temp_bboxes;
+};
+
+#endif
+
+// 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/libbrep/libbrep_util.cpp

@@ -0,0 +1,356 @@
+/*                L I B B R E P _ U T I L . C P P
+ * BRL-CAD
+ *
+ * Copyright (c) 2013 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.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file libbrep_util.cpp
+ *
+ * Brief description
+ *
+ */
+
+#include <vector>
+#include <iostream>
+
+#include "opennurbs.h"
+
+int ON_Surface_Patch_Trimmed(const ON_Surface *srf, ON_CurveTree *m_ctree, ON_Interval *u_val, ON_Interval *v_val)
+{
+    return 0;
+}
+
+int ON_Surface_Flat_Straight(ON_SimpleArray<ON_Plane> *frames, std::vector<int> *frame_index, double flatness_threshold, double straightness_threshold) {
+    double Fdot = 1.0;
+    double Sdot = 1.0;
+    // This uses 0 and 1 only - bug?
+    if ((Sdot = Sdot * frames->At(frame_index->at(0))->xaxis * frames->At(frame_index->at(1))->xaxis) < straightness_threshold) {
+	//std::cout << "Straightness: " << Sdot << "\n";
+	return 1;
+    }
+    for(int i=0; i<8; i++) {
+	for( int j=i+1; j<9; j++) {
+	    //std::cout << i << ": frames-At(" << frame_index->at(i) << ")" << "\n";
+	    //std::cout << j << ": frames-At(" << frame_index->at(j) << ")" << "\n";
+	    if ((Fdot = Fdot * frames->At(frame_index->at(i))->zaxis * frames->At(frame_index->at(j))->zaxis) < flatness_threshold) {
+		//std::cout << "Flatness: " << Fdot << "\n";
+		return 1;
+	    }
+	}
+    }
+    return 0;
+}
+
+// Directional flatness functions
+int ON_Surface_Flat_U(ON_SimpleArray<ON_Plane> *frames, std::vector<int> *f_ind, double flatness_threshold)
+{
+    // check surface normals in U direction
+    double Ndot = 1.0;
+    if ((Ndot=frames->At(f_ind->at(0))->zaxis * frames->At(f_ind->at(1))->zaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Ndot=Ndot * frames->At(f_ind->at(2))->zaxis * frames->At(f_ind->at(3))->zaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Ndot=Ndot * frames->At(f_ind->at(5))->zaxis * frames->At(f_ind->at(7))->zaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Ndot=Ndot * frames->At(f_ind->at(6))->zaxis * frames->At(f_ind->at(8))->zaxis) < flatness_threshold) {
+	return 1;
+    }
+
+    // check for U twist within plane
+    double Xdot = 1.0;
+    if ((Xdot=frames->At(f_ind->at(0))->xaxis * frames->At(f_ind->at(1))->xaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Xdot=Xdot * frames->At(f_ind->at(2))->xaxis * frames->At(f_ind->at(3))->xaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Xdot=Xdot * frames->At(f_ind->at(5))->xaxis * frames->At(f_ind->at(7))->xaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Xdot=Xdot * frames->At(f_ind->at(6))->xaxis * frames->At(f_ind->at(8))->xaxis) < flatness_threshold) {
+	return 1;
+    }
+
+    return 0;
+}
+
+int ON_Surface_Flat_V(ON_SimpleArray<ON_Plane> *frames, std::vector<int> *f_ind, double flatness_threshold)
+{
+    // check surface normals in V direction
+    double Ndot = 1.0;
+    if ((Ndot=frames->At(f_ind->at(0))->zaxis * frames->At(f_ind->at(3))->zaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Ndot=Ndot * frames->At(f_ind->at(1))->zaxis * frames->At(f_ind->at(2))->zaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Ndot=Ndot * frames->At(f_ind->at(5))->zaxis * frames->At(f_ind->at(6))->zaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Ndot=Ndot * frames->At(f_ind->at(7))->zaxis * frames->At(f_ind->at(8))->zaxis) < flatness_threshold) {
+	return 1;
+    }
+
+    // check for V twist within plane
+    double Xdot = 1.0;
+    if ((Xdot=frames->At(f_ind->at(0))->xaxis * frames->At(f_ind->at(3))->xaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Xdot=Xdot * frames->At(f_ind->at(1))->xaxis * frames->At(f_ind->at(2))->xaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Xdot=Xdot * frames->At(f_ind->at(5))->xaxis * frames->At(f_ind->at(6))->xaxis) < flatness_threshold) {
+	return 1;
+    } else if ((Xdot=Xdot * frames->At(f_ind->at(7))->xaxis * frames->At(f_ind->at(8))->xaxis) < flatness_threshold) {
+	return 1;
+    }
+
+    return 0;
+}
+
+
+// Check the ratio of the surface's "flattened" width to height
+int ON_Surface_Width_vs_Height(const ON_Surface *srf, double ratio, int *u_split, int *v_split)
+{
+    double width = 0;
+    double height = 0;
+    int split = 0;
+    double l_ratio = 0;
+    srf->GetSurfaceSize(&width, &height);
+    (ratio > 1.0) ? (l_ratio = ratio) : (l_ratio = 1/ratio);
+    if (width/height > l_ratio) {
+	(*u_split) += 1;
+	split++;
+    }
+     if (height/width > l_ratio) {
+	(*v_split) += 1;
+	split++;
+    }
+    return split;
+}
+
+
+// Check an individual domain and knot array for a split knot
+int ON_Interval_Find_Split_Knot(ON_Interval *val, int kcnt, double *knots, double *mid)
+{
+    int split = 0;
+    double midpt_dist = val->Length();
+
+    for (int k_ind = 1; k_ind <= kcnt; k_ind++) {
+	if (val->Includes(knots[k_ind], true)) {
+	    double d_to_midpt = fabs(val->Mid()-knots[k_ind]);
+	    if (d_to_midpt < midpt_dist) {
+		(*mid) = knots[k_ind];
+		midpt_dist = fabs(val->Mid()-knots[k_ind]);
+		split++;
+	    }
+	}
+    }
+
+    int retval = (split == 0) ? (0) : (1);
+    return retval;
+}
+
+// Figure out whether knot splitting is needed, and if so prepare values
+int ON_Surface_Knots_Split(ON_Interval *u_val, ON_Interval *v_val, int u_kcnt, int v_kcnt, double *u_knots, double *v_knots, double *u_mid, double *v_mid)
+{
+    int usplit = ON_Interval_Find_Split_Knot(u_val, u_kcnt, u_knots, u_mid);
+    int vsplit = ON_Interval_Find_Split_Knot(v_val, v_kcnt, v_knots, v_mid);
+
+    if (usplit && !vsplit) (*v_mid) = v_val->Mid();
+    if (!usplit && vsplit) (*u_mid) = u_val->Mid();
+
+    int retval = (usplit || vsplit) ? (1) : (0);
+    return retval;
+}
+
+// Prepare frames according to the standard node layout
+int ON_Populate_Frame_Array(const ON_Surface *srf, ON_Interval *u_domain, ON_Interval *v_domain, double umid, double vmid, ON_SimpleArray<ON_Plane> *frames, std::vector<int> *frame_index) {
+    ON_Plane frame;
+    if (frame_index->size() < 9) return -1;
+    double uqmin = fabs(umid - u_domain->Min())*0.25;
+    double uqmax = fabs(u_domain->Max() - umid)*0.25;
+    double vqmin = fabs(vmid - v_domain->Min())*0.25;
+    double vqmax = fabs(v_domain->Max() - vmid)*0.25;
+    if (frame_index->at(0) == -1) {
+	srf->FrameAt(u_domain->Min(), v_domain->Min(), frame);
+	frames->Append(frame);
+	frame_index->at(0) = frames->Count() - 1;
+    }
+    if (frame_index->at(1) == -1) {
+	srf->FrameAt(u_domain->Max(), v_domain->Min(), frame);
+	frames->Append(frame);
+	frame_index->at(1) = frames->Count() - 1;
+    }
+    if (frame_index->at(2) == -1) {
+	srf->FrameAt(u_domain->Max(), v_domain->Max(), frame);
+	frames->Append(frame);
+	frame_index->at(2) = frames->Count() - 1;
+    }
+    if (frame_index->at(3) == -1) {
+	srf->FrameAt(u_domain->Min(), v_domain->Max(), frame);
+	frames->Append(frame);
+	frame_index->at(3) = frames->Count() - 1;
+    }
+    if (frame_index->at(4) == -1) {
+	srf->FrameAt(umid, vmid, frame);
+	frames->Append(frame);
+	frame_index->at(4) = frames->Count() - 1;
+    }
+    if (frame_index->at(5) == -1) {
+	srf->FrameAt(umid - uqmin, vmid - vqmin, frame);
+	frames->Append(frame);
+	frame_index->at(5) = frames->Count() - 1;
+    }
+    if (frame_index->at(6) == -1) {
+	srf->FrameAt(umid - uqmin, vmid + vqmax, frame);
+	frames->Append(frame);
+	frame_index->at(6) = frames->Count() - 1;
+    }
+    if (frame_index->at(7) == -1) {
+	srf->FrameAt(umid + uqmax, vmid - vqmin, frame);
+	frames->Append(frame);
+	frame_index->at(7) = frames->Count() - 1;
+    }
+    if (frame_index->at(8) == -1) {
+	srf->FrameAt(umid + uqmax, vmid + vqmax, frame);
+	frames->Append(frame);
+	frame_index->at(8) = frames->Count() - 1;
+    }
+    return 0;
+}
+
+// Prepare frames according to the standard node layout
+int ON_Extend_Frame_Array(const ON_Surface *srf, ON_Interval *u_domain, ON_Interval *v_domain, double umid, double vmid, ON_SimpleArray<ON_Plane> *frames, std::vector<int> *frame_index) {
+    ON_Plane frame;
+    if (frame_index->size() < 13) return -1;
+    if (frame_index->at(9) == -1) {
+	srf->FrameAt(umid, v_domain->Min(), frame);
+	frames->Append(frame);
+	frame_index->at(9) = frames->Count() - 1;
+    }
+    if (frame_index->at(10) == -1) {
+	srf->FrameAt(u_domain->Min(), vmid, frame);
+	frames->Append(frame);
+	frame_index->at(10) = frames->Count() - 1;
+    }
+    if (frame_index->at(11) == -1) {
+	srf->FrameAt(umid, v_domain->Max(), frame);
+	frames->Append(frame);
+	frame_index->at(11) = frames->Count() - 1;
+    }
+    if (frame_index->at(12) == -1) {
+	srf->FrameAt(u_domain->Max(), vmid, frame);
+	frames->Append(frame);
+	frame_index->at(12) = frames->Count() - 1;
+    }
+    return 0;
+}
+
+// Check to see whether a 2D point is inside or outside of a bounding box.
+bool ON_BoundingBox_Contains_2DPoint(ON_BoundingBox *bbox, double u, double v) {
+    if (u < bbox->m_min[0] || u > bbox->m_max[0] || v < bbox->m_min[1] || v > bbox->m_max[1]) return false;
+    return true;
+}
+
+
+// For any pre-existing surface passed as one of the t* args, this is a no-op
+void ON_Surface_Create_Scratch_Surfaces(
+	ON_Surface **t1,
+	ON_Surface **t2,
+	ON_Surface **t3,
+	ON_Surface **t4)
+{
+    if (!(*t1)) {
+	ON_NurbsSurface *nt1 = ON_NurbsSurface::New();
+	(*t1)= (ON_Surface *)(nt1);
+    }
+    if (!(*t2)) {
+	ON_NurbsSurface *nt2 = ON_NurbsSurface::New();
+	(*t2)= (ON_Surface *)(nt2);
+    }
+    if (!(*t3)) {
+	ON_NurbsSurface *nt3 = ON_NurbsSurface::New();
+	(*t3)= (ON_Surface *)(nt3);
+    }
+    if (!(*t4)) {
+	ON_NurbsSurface *nt4 = ON_NurbsSurface::New();
+	(*t4)= (ON_Surface *)(nt4);
+    }
+}
+
+
+// Given a surface and UV intervals, return a NURBS surface corresponding to
+// that subset of the patch.  If t1-t3 and result are supplied externally,
+// they are re-used - if not, local versions must be created and destroyed.
+// The latter usage will have a malloc overhead penalty.
+//
+// Returns true if splits successful (or if none were needed), false if one
+// or more splits failed
+bool ON_Surface_SubSurface(
+	const ON_Surface *srf,
+	ON_Surface **result,
+	ON_Interval *u_val,
+	ON_Interval *v_val,
+	ON_Surface **t1,
+	ON_Surface **t2,
+	ON_Surface **t3,
+	ON_Surface **t4)
+{
+    bool split = true;
+    int t1_del, t2_del, t3_del;
+
+    // Make sure we have intervals with non-zero lengths
+    if ((u_val->Length() <= ON_ZERO_TOLERANCE) || (v_val->Length() <= ON_ZERO_TOLERANCE))
+	return false;
+
+    // If we have the original surface domain, just return true
+    if ((fabs(u_val->Min() - srf->Domain(0).m_t[0]) <= ON_ZERO_TOLERANCE) &&
+        (fabs(u_val->Max() - srf->Domain(0).m_t[1]) <= ON_ZERO_TOLERANCE) &&
+        (fabs(v_val->Min() - srf->Domain(1).m_t[0]) <= ON_ZERO_TOLERANCE) &&
+	(fabs(v_val->Max() - srf->Domain(1).m_t[1]) <= ON_ZERO_TOLERANCE)) {
+        (*result) = (ON_Surface *)srf;
+	return true;
+    }
+    t1_del = (*t1) ? (0) : (1);
+    t2_del = (*t2) ? (0) : (1);
+    t3_del = (*t3) ? (0) : (1);
+    ON_Surface *ssplit = (ON_Surface *)srf;
+    ON_Surface_Create_Scratch_Surfaces(t1, t2, t3, t4);
+    if (fabs(u_val->Min() - srf->Domain(0).m_t[0]) > ON_ZERO_TOLERANCE) {
+	split = ssplit->Split(0, u_val->Min(), *t1, *t2);
+	ssplit = *t2;
+    }
+    if ((fabs(u_val->Max() - srf->Domain(0).m_t[1]) > ON_ZERO_TOLERANCE) && split) {
+	split = ssplit->Split(0, u_val->Max(), *t1, *t3);
+	ssplit = *t1;
+    }
+    if ((fabs(v_val->Min() - srf->Domain(1).m_t[0]) > ON_ZERO_TOLERANCE) && split) {
+	split = ssplit->Split(1, v_val->Min(), *t2, *t3);
+        ssplit = *t3;
+    }
+    if ((fabs(v_val->Max() - srf->Domain(1).m_t[1]) > ON_ZERO_TOLERANCE) && split) {
+	split = ssplit->Split(1, v_val->Max(), *t4, *t2);
+        (*result) = *t4;
+    }
+    if (t1_del) delete *t1;
+    if (t2_del) delete *t2;
+    if (t3_del) delete *t3;
+    (*result)->SetDomain(0,u_val->Min(), u_val->Max());
+    (*result)->SetDomain(1,v_val->Min(), v_val->Max());
+    return split;
+}
+
+// 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/libbrep/libbrep_util.h

@@ -0,0 +1,159 @@
+/*      L I B B R E P _ S U R F A C E T R E E U T I L . H
+ * BRL-CAD
+ *
+ * Copyright (c) 2013 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.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this file; see the file named COPYING for more
+ * information.
+ */
+/** @file libbrep_util.h
+ *
+ * Brief description
+ *
+ */
+
+#ifndef __LIBBREP_UTILS
+#define __LIBBREP_UTILS
+
+#include <vector>
+
+#include "opennurbs.h"
+
+#ifndef NURBS_EXPORT
+#  if defined(NURBS_DLL_EXPORTS) && defined(NURBS_DLL_IMPORTS)
+#    error "Only NURBS_DLL_EXPORTS or NURBS_DLL_IMPORTS can be defined, not both."
+#  elif defined(NURBS_DLL_EXPORTS)
+#    define NURBS_EXPORT __declspec(dllexport)
+#  elif defined(NURBS_DLL_IMPORTS)
+#    define NURBS_EXPORT __declspec(dllimport)
+#  else
+#    define NURBS_EXPORT
+#  endif
+#endif
+
+// Determines the trimming status of a given patch
+extern NURBS_EXPORT int ON_Surface_Patch_Trimmed(
+	const ON_Surface *srf,
+	ON_CurveTree *c_tree,
+	ON_Interval *u_val,
+	ON_Interval *v_val);
+
+// Test a frame set from a surface patch for flatness and straightness
+extern NURBS_EXPORT int ON_Surface_Flat_Straight(
+	ON_SimpleArray<ON_Plane> *frames,
+	std::vector<int> *frame_index,
+	double flatness_threshold,
+	double straightness_threshold);
+// Directional flatness functions
+extern NURBS_EXPORT int ON_Surface_Flat_U(
+	ON_SimpleArray<ON_Plane> *frames,
+	std::vector<int> *frame_index,
+	double flatness_threshold);
+extern NURBS_EXPORT int ON_Surface_Flat_V(
+	ON_SimpleArray<ON_Plane> *frames,
+	std::vector<int> *frame_index,
+	double flatness_threshold);
+
+// Check flattened width/height ratio
+extern NURBS_EXPORT int ON_Surface_Width_vs_Height(const ON_Surface *srf, double ratio, int *u_only, int *v_only);
+
+// Check an individual domain and knot array for a split knot
+extern NURBS_EXPORT int ON_Interval_Find_Split_Knot(ON_Interval *val, int kcnt, double *knots, double *mid);
+
+// Figure out whether knot splitting is needed, and if so prepare values
+extern NURBS_EXPORT int ON_Surface_Knots_Split(
+	ON_Interval *u_val,
+	ON_Interval *v_val,
+	int u_kcnt,
+	int v_kcnt,
+	double *u_knots,
+	double *v_knots,
+	double *u_mid,
+	double *v_mid);
+
+
+// Prepare frames according to the standard node layout:
+//
+//          3---------11--------2
+//          |                   |
+//          |    6         8    |
+//          |                   |
+//        V 10        4        12
+//          |                   |
+//          |    5         7    |
+//          |                   |
+//          0---------9---------1
+//                    U
+
+extern NURBS_EXPORT int ON_Populate_Frame_Array(
+	const ON_Surface *srf,
+	ON_Interval *u_domain,
+	ON_Interval *v_domain,
+	double umid,
+	double vmid,
+	ON_SimpleArray<ON_Plane> *frames,
+	std::vector<int> *frame_index);
+
+// Extend frames according to the standard node layout for splitting
+extern NURBS_EXPORT int ON_Extend_Frame_Array(
+	const ON_Surface *srf,
+	ON_Interval *u_domain,
+	ON_Interval *v_domain,
+	double umid,
+	double vmid,
+	ON_SimpleArray<ON_Plane> *frames,
+	std::vector<int> *frame_index);
+
+// Test weather a UV point within a Bounding Box
+extern NURBS_EXPORT bool ON_BoundingBox_Contains_2DPoint(
+	ON_BoundingBox *bbox,
+	double u,
+	double v);
+
+// For any pre-existing surface passed as one of the t* args, this is a no-op
+extern NURBS_EXPORT void ON_Surface_Create_Scratch_Surfaces(
+	ON_Surface **t1,
+	ON_Surface **t2,
+	ON_Surface **t3,
+	ON_Surface **t4);
+
+
+// Given a surface and UV intervals, return a NURBS surface corresponding to
+// that subset of the patch.  If t1-t3 and result are supplied externally,
+// they are re-used - if not, local versions will be created and destroyed.
+// The latter usage will have a malloc overhead penalty.  It is always the
+// responsibility of the caller to delete the result surface.  If temporary
+// surfaces are passed in, deleting those is also the responsibility of the
+// caller.
+extern NURBS_EXPORT bool ON_Surface_SubSurface(
+	const ON_Surface *srf,
+	ON_Surface **result,
+	ON_Interval *u_val,
+	ON_Interval *v_val,
+	ON_Surface **t1,
+	ON_Surface **t2,
+	ON_Surface **t3,
+	ON_Surface **t4);
+
+
+#endif
+
+// 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/CMakeLists.txt

@@ -311,7 +311,7 @@ include_directories(
 #BRLCAD_ADDEXEC(test_bot2nurbs test_bot2nurbs.cpp "librt;libwdb;libbrep;libbu" NO_STRICT NO_INSTALL)
 #BRLCAD_ADDEXEC(test_nurbsfit test_nurbsfit.cpp "librt;libwdb;libbrep;libbu" NO_STRICT NO_INSTALL)
 #BRLCAD_ADDEXEC(test_root3-subd test_root3-subd.cpp "librt;libwdb;libbrep;libbu" NO_STRICT NO_INSTALL)
-#BRLCAD_ADDEXEC(test_surfacetree test_surfacetree.cpp "librt;libwdb;libbrep;libbu" NO_STRICT NO_INSTALL)
+BRLCAD_ADDEXEC(test_surfacetree test_surfacetree.cpp "librt;libwdb;libbrep;libbu" NO_STRICT NO_INSTALL)
 
 add_subdirectory(tests)
 

src/librt/test_surfacetree.cpp

@@ -0,0 +1,75 @@
+#include "common.h"
+#include "vmath.h"
+#include "raytrace.h"
+#include "wdb.h"
+#include "brep.h"
+
+#include "opennurbs.h"
+#include "libbrep_surfacetree.h"
+
+int
+main (int argc, char *argv[])
+{
+   struct db_i *dbip;
+   struct directory *dp;
+   struct rt_db_internal intern;
+   struct rt_brep_internal *brep_ip;
+   struct rt_wdb *wdbp;
+
+   if (argc != 3) {
+      bu_exit(1, "Usage: %s file.g object", argv[0]);
+   }
+
+   dbip = db_open(argv[1], "r+w");
+   if (dbip == DBI_NULL) {
+      bu_exit(1, "ERROR: Unable to read from %s\n", argv[1]);
+   }
+
+   if (db_dirbuild(dbip) < 0)
+      bu_exit(1, "ERROR: Unable to read from %s\n", argv[1]);
+
+   dp = db_lookup(dbip, argv[2], LOOKUP_QUIET);
+   if (dp == RT_DIR_NULL) {
+      bu_exit(1, "ERROR: Unable to look up object %s\n", argv[2]);
+   } else {
+      bu_log("Got %s\n", dp->d_namep);
+   }
+
+   RT_DB_INTERNAL_INIT(&intern)
+   if (rt_db_get_internal(&intern, dp, dbip, NULL, &rt_uniresource) < 0) {
+      bu_exit(1, "ERROR: Unable to get internal representation of %s\n", argv[2]);
+   }
+
+   if (intern.idb_minor_type != DB5_MINORTYPE_BRLCAD_BREP) {
+      bu_exit(1, "ERROR: object %s does not appear to be of type BREP\n", argv[2]);
+   } else {
+      brep_ip = (struct rt_brep_internal *)intern.idb_ptr;
+   }
+
+    ON_Brep *brep = brep_ip->brep;
+    struct bu_vls vls;
+    bu_vls_init(&vls);
+    bu_vls_printf(&vls, "surfaces:  %d\n", brep->m_S.Count());
+    bu_vls_printf(&vls, "3d curve:  %d\n", brep->m_C3.Count());
+    bu_vls_printf(&vls, "2d curves: %d\n", brep->m_C2.Count());
+    bu_vls_printf(&vls, "vertices:  %d\n", brep->m_V.Count());
+    bu_vls_printf(&vls, "edges:     %d\n", brep->m_E.Count());
+    bu_vls_printf(&vls, "trims:     %d\n", brep->m_T.Count());
+    bu_vls_printf(&vls, "loops:     %d\n", brep->m_L.Count());
+    bu_vls_printf(&vls, "faces:     %d\n", brep->m_F.Count());
+    printf("%s\n", bu_vls_addr(&vls));
+    bu_vls_free(&vls);
+
+    for(int i = 0; i < brep->m_F.Count(); i++) {
+       ON_BrepFace& face = brep->m_F[i];
+       std::cout << "Face: " << i << "\n";
+       ON_SurfaceTree *stree = new ON_SurfaceTree(&face);
+       delete stree;
+       //std::cout << "Face(old): " << i << "\n";
+       //brlcad::SurfaceTree *st = new brlcad::SurfaceTree(&face);
+    }
+
+    db_close(dbip);
+
+   return 0;
+}