BRL-CAD
opennurbs_ext.cpp
Go to the documentation of this file.
1 /* O P E N N U R B S _ E X T . C P P
2  * BRL-CAD
3  *
4  * Copyright (c) 2007-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @file opennurbs_ext.cpp
21  *
22  * Implementation of routines openNURBS left out.
23  *
24  */
25 
26 #include "common.h"
27 
28 #include "bio.h"
29 #include <assert.h>
30 #include <vector>
31 
32 #include "vmath.h"
33 
34 #include "bu/log.h"
35 #include "bu/malloc.h"
36 #include "bu/parallel.h"
37 #include "brep.h"
38 #include "libbrep_brep_tools.h"
39 #include "dvec.h"
40 
41 #define RANGE_HI 0.55
42 #define RANGE_LO 0.45
43 #define UNIVERSAL_SAMPLE_COUNT 1001
44 
45 #define BBOX_GROW 0.0
46 
47 /// grows 3D BBox along each axis by this factor
48 #define BBOX_GROW_3D 0.1
49 
50 /// arbitrary calculation tolerance (need to try VDIVIDE_TOL or VUNITIZE_TOL to tighten the bounds)
51 #define TOL 0.000001
52 
53 /// another arbitrary calculation tolerance (need to try VDIVIDE_TOL or VUNITIZE_TOL to tighten the bounds)
54 #define TOL2 0.00001
55 
56 //#define _OLD_SUBDIVISION_
57 
58 void
59 brep_get_plane_ray(ON_Ray& r, plane_ray& pr)
60 {
61  vect_t v1;
62  VMOVE(v1, r.m_dir);
63  fastf_t min = MAX_FASTF;
64  int index = -1;
65  for (int i = 0; i < 3; i++) {
66  // find the smallest component
67  if (fabs(v1[i]) < min) {
68  min = fabs(v1[i]);
69  index = i;
70  }
71  }
72  v1[index] += 1; // alter the smallest component
73  VCROSS(pr.n1, v1, r.m_dir); // n1 is perpendicular to v1
74  VUNITIZE(pr.n1);
75  VCROSS(pr.n2, pr.n1, r.m_dir); // n2 is perpendicular to v1 and n1
76  VUNITIZE(pr.n2);
77  pr.d1 = VDOT(pr.n1, r.m_origin);
78  pr.d2 = VDOT(pr.n2, r.m_origin);
79  TRACE1("n1:" << ON_PRINT3(pr.n1) << " n2:" << ON_PRINT3(pr.n2) << " d1:" << pr.d1 << " d2:" << pr.d2);
80 }
81 
82 
83 void
84 brep_r(const ON_Surface* surf, plane_ray& pr, pt2d_t uv, ON_3dPoint& pt, ON_3dVector& su, ON_3dVector& sv, pt2d_t R)
85 {
86  vect_t vp;
87 
88  assert(surf->Ev1Der(uv[0], uv[1], pt, su, sv));
89  VMOVE(vp, pt);
90 
91  R[0] = VDOT(pr.n1, vp) - pr.d1;
92  R[1] = VDOT(pr.n2, vp) - pr.d2;
93 }
94 
95 
96 void
97 brep_newton_iterate(plane_ray& pr, pt2d_t R, ON_3dVector& su, ON_3dVector& sv, pt2d_t uv, pt2d_t out_uv)
98 {
99  vect_t vsu, vsv;
100  VMOVE(vsu, su);
101  VMOVE(vsv, sv);
102 
103  mat2d_t jacob = { VDOT(pr.n1, vsu), VDOT(pr.n1, vsv),
104  VDOT(pr.n2, vsu), VDOT(pr.n2, vsv) };
105  mat2d_t inv_jacob;
106  if (mat2d_inverse(inv_jacob, jacob)) {
107  // check inverse validity
108  pt2d_t tmp;
109  mat2d_pt2d_mul(tmp, inv_jacob, R);
110  pt2dsub(out_uv, uv, tmp);
111  } else {
112  TRACE2("inverse failed"); // FIXME: how to handle this?
113  move(out_uv, uv);
114  }
115 }
116 
117 
118 namespace brlcad {
119 
120 inline void
121 distribute(const int count, const ON_3dVector* v, double x[], double y[], double z[])
122 {
123  for (int i = 0; i < count; i++) {
124  x[i] = v[i].x;
125  y[i] = v[i].y;
126  z[i] = v[i].z;
127  }
128 }
129 
130 
131 //--------------------------------------------------------------------------------
132 // CurveTree
133 CurveTree::CurveTree(const ON_BrepFace* face) :
134  m_face(face)
135 {
136  m_root = initialLoopBBox();
137 
138  for (int li = 0; li < face->LoopCount(); li++) {
139  bool innerLoop = (li > 0) ? true : false;
140  ON_BrepLoop* loop = face->Loop(li);
141  // for each trim
142  for (int ti = 0; ti < loop->m_ti.Count(); ti++) {
143  int adj_face_index = -99;
144  ON_BrepTrim& trim = face->Brep()->m_T[loop->m_ti[ti]];
145 
146  if (trim.m_ei != -1) { // does not lie on a portion of a singular surface side
147  ON_BrepEdge& edge = face->Brep()->m_E[trim.m_ei];
148  switch (trim.m_type) {
149  case ON_BrepTrim::unknown:
150  bu_log("ON_BrepTrim::unknown on Face:%d\n", face->m_face_index);
151  break;
152  case ON_BrepTrim::boundary:
153  //bu_log("ON_BrepTrim::boundary on Face:%d\n", face->m_face_index);
154  break;
155  case ON_BrepTrim::mated:
156  if (edge.m_ti.Count() == 2) {
157  if (face->m_face_index == face->Brep()->m_T[edge.m_ti[0]].FaceIndexOf()) {
158  adj_face_index = face->Brep()->m_T[edge.m_ti[1]].FaceIndexOf();
159  } else {
160  adj_face_index = face->Brep()->m_T[edge.m_ti[0]].FaceIndexOf();
161  }
162  } else {
163  bu_log("Mated Edge should have 2 adjacent faces, right? Face(%d) has %d trim indexes\n", face->m_face_index, edge.m_ti.Count());
164  }
165  break;
166  case ON_BrepTrim::seam:
167  if (edge.m_ti.Count() == 2) {
168  if ((face->m_face_index == face->Brep()->m_T[edge.m_ti[0]].FaceIndexOf()) && (face->m_face_index == face->Brep()->m_T[edge.m_ti[1]].FaceIndexOf())) {
169  adj_face_index = face->m_face_index;
170  } else {
171  bu_log("Seamed Edge should have 1 faces sharing the trim so trim index should be one, right? Face(%d) has %d trim indexes\n", face->m_face_index, edge.m_ti.Count());
172  bu_log("Face(%d) has %d, %d trim indexes\n", face->m_face_index, face->Brep()->m_T[edge.m_ti[0]].FaceIndexOf(), face->Brep()->m_T[edge.m_ti[1]].FaceIndexOf());
173  }
174  } else if (edge.m_ti.Count() == 1) {
175  adj_face_index = face->m_face_index;
176  } else {
177  bu_log("Seamed Edge should have 1 faces sharing the trim so trim index should be one, right? Face(%d) has %d trim indexes\n", face->m_face_index, edge.m_ti.Count());
178  }
179  break;
180  case ON_BrepTrim::singular:
181  bu_log("ON_BrepTrim::singular on Face:%d\n", face->m_face_index);
182  break;
183  case ON_BrepTrim::crvonsrf:
184  bu_log("ON_BrepTrim::crvonsrf on Face:%d\n", face->m_face_index);
185  break;
186  case ON_BrepTrim::ptonsrf:
187  bu_log("ON_BrepTrim::ptonsrf on Face:%d\n", face->m_face_index);
188  break;
189  case ON_BrepTrim::slit:
190  bu_log("ON_BrepTrim::slit on Face:%d\n", face->m_face_index);
191  break;
192  default:
193  bu_log("ON_BrepTrim::default on Face:%d\n", face->m_face_index);
194  }
195  }
196  const ON_Curve* trimCurve = trim.TrimCurveOf();
197  double min, max;
198  (void) trimCurve->GetDomain(&min, &max);
199  ON_Interval t(min, max);
200 
201  TRACE("need to subdivide");
202  // divide on param interval
203 
204  if (!trimCurve->IsLinear()) {
205  int knotcnt = trimCurve->SpanCount();
206  double *knots = new double[knotcnt + 1];
207 
208  trimCurve->GetSpanVector(knots);
209  std::list<fastf_t> splitlist;
210  for (int knot_index = 1; knot_index <= knotcnt; knot_index++) {
211  ON_Interval range(knots[knot_index - 1], knots[knot_index]);
212 
213  if (range.Length() > TOL)
214  getHVTangents(trimCurve, range, splitlist);
215  }
216  for (std::list<fastf_t>::iterator l = splitlist.begin(); l != splitlist.end(); l++) {
217  double xmax = *l;
218  if (!NEAR_EQUAL(xmax, min, TOL)) {
219  m_root->addChild(subdivideCurve(trimCurve, adj_face_index, min, xmax, innerLoop, 0));
220  }
221  min = xmax;
222  }
223  delete [] knots;
224  } else {
225  int knotcnt = trimCurve->SpanCount();
226  double *knots = new double[knotcnt + 1];
227 
228  trimCurve->GetSpanVector(knots);
229  for (int knot_index = 1; knot_index <= knotcnt; knot_index++) {
230  double xmax = knots[knot_index];
231  if (!NEAR_EQUAL(xmax, min, TOL)) {
232  m_root->addChild(subdivideCurve(trimCurve, adj_face_index, min, xmax, innerLoop, 0));
233  }
234  min = xmax;
235  }
236  delete [] knots;
237  }
238 
239  if (!NEAR_EQUAL(max, min, TOL)) {
240  m_root->addChild(subdivideCurve(trimCurve, adj_face_index, min, max, innerLoop, 0));
241  }
242  }
243  }
244  getLeaves(m_sortedX);
245  m_sortedX.sort(sortX);
246  getLeaves(m_sortedY);
247  m_sortedY.sort(sortY);
248 
249  return;
250 }
251 
252 
253 CurveTree::~CurveTree()
254 {
255  delete m_root;
256 }
257 
258 
259 BRNode*
260 CurveTree::getRootNode() const
261 {
262  return m_root;
263 }
264 
265 
266 int
267 CurveTree::depth()
268 {
269  return m_root->depth();
270 }
271 
272 
273 ON_2dPoint
274 CurveTree::getClosestPointEstimate(const ON_3dPoint& pt)
275 {
276  return m_root->getClosestPointEstimate(pt);
277 }
278 
279 
280 ON_2dPoint
281 CurveTree::getClosestPointEstimate(const ON_3dPoint& pt, ON_Interval& u, ON_Interval& v)
282 {
283  return m_root->getClosestPointEstimate(pt, u, v);
284 }
285 
286 
287 void
288 CurveTree::getLeaves(std::list<BRNode*>& out_leaves)
289 {
290  m_root->getLeaves(out_leaves);
291 }
292 
293 
294 void
295 CurveTree::getLeavesAbove(std::list<BRNode*>& out_leaves, const ON_Interval& u, const ON_Interval& v)
296 {
297  point_t bmin, bmax;
298  double dist;
299  for (std::list<BRNode*>::iterator i = m_sortedX.begin(); i != m_sortedX.end(); i++) {
300  BRNode* br = dynamic_cast<BRNode*>(*i);
301  br->GetBBox(bmin, bmax);
302 
303  dist = TOL;//0.03*DIST_PT_PT(bmin, bmax);
304  if (bmax[X]+dist < u[0])
305  continue;
306  if (bmin[X]-dist < u[1]) {
307  if (bmax[Y]+dist > v[0]) {
308  out_leaves.push_back(br);
309  }
310  }
311  }
312 }
313 
314 
315 void
316 CurveTree::getLeavesAbove(std::list<BRNode*>& out_leaves, const ON_2dPoint& pt, fastf_t tol)
317 {
318  point_t bmin, bmax;
319  for (std::list<BRNode*>::iterator i = m_sortedX.begin(); i != m_sortedX.end(); i++) {
320  BRNode* br = dynamic_cast<BRNode*>(*i);
321  br->GetBBox(bmin, bmax);
322 
323  if (bmax[X]+tol < pt.x)
324  continue;
325  if (bmin[X]-tol < pt.x) {
326  if (bmax[Y]+tol > pt.y) {
327  out_leaves.push_back(br);
328  }
329  }
330  }
331 }
332 
333 
334 void
335 CurveTree::getLeavesRight(std::list<BRNode*>& out_leaves, const ON_Interval& u, const ON_Interval& v)
336 {
337  point_t bmin, bmax;
338  double dist;
339  for (std::list<BRNode*>::iterator i = m_sortedX.begin(); i != m_sortedX.end(); i++) {
340  BRNode* br = dynamic_cast<BRNode*>(*i);
341  br->GetBBox(bmin, bmax);
342 
343  dist = TOL;//0.03*DIST_PT_PT(bmin, bmax);
344  if (bmax[Y]+dist < v[0])
345  continue;
346  if (bmin[Y]-dist < v[1]) {
347  if (bmax[X]+dist > u[0]) {
348  out_leaves.push_back(br);
349  }
350  }
351  }
352 }
353 
354 
355 void
356 CurveTree::getLeavesRight(std::list<BRNode*>& out_leaves, const ON_2dPoint& pt, fastf_t tol)
357 {
358  point_t bmin, bmax;
359  for (std::list<BRNode*>::iterator i = m_sortedX.begin(); i != m_sortedX.end(); i++) {
360  BRNode* br = dynamic_cast<BRNode*>(*i);
361  br->GetBBox(bmin, bmax);
362 
363  if (bmax[Y]+tol < pt.y)
364  continue;
365  if (bmin[Y]-tol < pt.y) {
366  if (bmax[X]+tol > pt.x) {
367  out_leaves.push_back(br);
368  }
369  }
370  }
371 }
372 
373 
374 bool
375 CurveTree::getHVTangents(const ON_Curve* curve, ON_Interval& t, std::list<fastf_t>& list)
376 {
377  double x;
378  double midpoint = (t[1]+t[0])/2.0;
379  ON_Interval left(t[0], midpoint);
380  ON_Interval right(midpoint, t[1]);
381  int status = ON_Curve_Has_Tangent(curve, t[0], t[1], TOL);
382 
383  switch (status) {
384 
385  case 1: /* 1 Vertical tangent */
386  x = ON_Curve_Get_Vertical_Tangent(curve, t[0], t[1], TOL2);
387  list.push_back(x);
388  return true;
389 
390  case 2: /* 1 Horizontal tangent */
391  x = ON_Curve_Get_Horizontal_Tangent(curve, t[0], t[1], TOL2);
392  list.push_back(x);
393  return true;
394 
395  case 3: /* Horizontal and vertical tangents present - Simple midpoint split */
396  if (left.Length() > TOL)
397  getHVTangents(curve, left, list);
398  if (right.Length() > TOL)
399  getHVTangents(curve, right, list);
400  return true;
401 
402  default:
403  return false;
404 
405  }
406 
407  return false; //Should never get here
408 }
409 
410 
411 BRNode*
412 CurveTree::curveBBox(const ON_Curve* curve, int adj_face_index, ON_Interval& t, bool isLeaf, bool innerTrim, const ON_BoundingBox& bb)
413 {
414  BRNode* node;
415  bool vdot = true;
416 
417  if (isLeaf) {
418  TRACE("creating leaf: u(" << u.Min() << ", " << u.Max() << ") v(" << v.Min() << ", " << v.Max() << ")");
419  node = new BRNode(curve, adj_face_index, bb, m_face, t, vdot, innerTrim);
420  } else {
421  node = new BRNode(bb);
422  }
423 
424  return node;
425 
426 }
427 
428 
429 BRNode*
430 CurveTree::initialLoopBBox()
431 {
432  ON_BoundingBox bb;
433  m_face->SurfaceOf()->GetBBox(bb[0], bb[1]);
434 
435  for (int i = 0; i < m_face->LoopCount(); i++) {
436  ON_BrepLoop* loop = m_face->Loop(i);
437  if (loop->m_type == ON_BrepLoop::outer) {
438  if (loop->GetBBox(bb[0], bb[1], 0)) {
439  TRACE("BBox for Loop min<" << bb[0][0] << ", " << bb[0][1] ", " << bb[0][2] << ">");
440  TRACE("BBox for Loop max<" << bb[1][0] << ", " << bb[1][1] ", " << bb[1][2] << ">");
441  }
442  break;
443  }
444  }
445  BRNode* node = new BRNode(bb);
446  return node;
447 }
448 
449 
450 BRNode*
451 CurveTree::subdivideCurve(const ON_Curve* curve, int adj_face_index, double min, double max, bool innerTrim, int divDepth)
452 {
453  ON_Interval dom = curve->Domain();
454  ON_3dPoint points[2];
455  points[0] = curve->PointAt(min);
456  points[1] = curve->PointAt(max);
457  point_t minpt, maxpt;
458  VSETALL(minpt, INFINITY);
459  VSETALL(maxpt, -INFINITY);
460  for (int i = 0; i < 2; i++)
461  VMINMAX(minpt, maxpt, ((double*)points[i]));
462  points[0]=ON_3dPoint(minpt);
463  points[1]=ON_3dPoint(maxpt);
464  ON_BoundingBox bb(points[0], points[1]);
465 
466  ON_Interval t(min, max);
467  if (isLinear(curve, min, max) || divDepth >= BREP_MAX_LN_DEPTH) {
468  double delta = (max - min)/(BREP_BB_CRV_PNT_CNT-1);
469  point_t pnts[BREP_BB_CRV_PNT_CNT];
470  ON_3dPoint pnt;
471  for (int i=0;i<BREP_BB_CRV_PNT_CNT-1;i++) {
472  pnt = curve->PointAt(min + delta*i);
473  VSET(pnts[i], pnt[0], pnt[1], pnt[2]);
474  }
475  pnt = curve->PointAt(max);
476  VSET(pnts[BREP_BB_CRV_PNT_CNT-1], pnt[0], pnt[1], pnt[2]);
477 
478  VSETALL(minpt, MAX_FASTF);
479  VSETALL(maxpt, -MAX_FASTF);
480  for (int i = 0; i < BREP_BB_CRV_PNT_CNT; i++)
481  VMINMAX(minpt, maxpt, ((double*)pnts[i]));
482 
483  VMOVE(pnt, minpt);
484  bb.Set(pnt, false);
485  VMOVE(pnt, maxpt);
486  bb.Set(pnt, true);
487  return curveBBox(curve, adj_face_index, t, true, innerTrim, bb);
488  }
489 
490  // else subdivide
491  BRNode* parent = curveBBox(curve, adj_face_index, t, false, innerTrim, bb);
492  double mid = (max+min)/2.0;
493  BRNode* l = subdivideCurve(curve, adj_face_index, min, mid, innerTrim, divDepth+1);
494  BRNode* r = subdivideCurve(curve, adj_face_index, mid, max, innerTrim, divDepth+1);
495  parent->addChild(l);
496  parent->addChild(r);
497  return parent;
498 }
499 
500 
501 /**
502  * Determine whether a given curve segment is linear
503  */
504 bool
505 CurveTree::isLinear(const ON_Curve* curve, double min, double max)
506 {
507  ON_3dVector tangent_start = curve->TangentAt(min);
508  ON_3dVector tangent_end = curve->TangentAt(max);
509  double vdot = tangent_start * tangent_end;
510  if (vdot < BREP_CURVE_FLATNESS)
511  return false;
512 
513  ON_3dPoint pmin = curve->PointAt(min);
514  ON_3dPoint pmax = curve->PointAt(max);
515 
516  const ON_Surface* surf = m_face->SurfaceOf();
517  ON_Interval u = surf->Domain(0);
518  ON_Interval v = surf->Domain(1);
519  point_t a, b;
520  VSET(a, u[0], v[0], 0.0);
521  VSET(b, u[1], v[1], 0.0);
522  double dd = DIST_PT_PT(a, b);
523  double cd = DIST_PT_PT(pmin, pmax);
524 
525  if (cd > BREP_TRIM_SUB_FACTOR*dd)
526  return false;
527 
528  double delta = (max - min)/(BREP_BB_CRV_PNT_CNT-1);
529  ON_3dPoint points[BREP_BB_CRV_PNT_CNT];
530  for (int i=0;i<BREP_BB_CRV_PNT_CNT-1;i++) {
531  points[i] = curve->PointAt(min + delta*i);
532  }
533  points[BREP_BB_CRV_PNT_CNT-1] = curve->PointAt(max);
534 
535  ON_3dVector A;
536  ON_3dVector B;
537  vdot = 1.0;
538  A = points[BREP_BB_CRV_PNT_CNT-1] - points[0];
539  A.Unitize();
540  for (int i=1;i<BREP_BB_CRV_PNT_CNT-1;i++) {
541  B = points[i] - points[0];
542  B.Unitize();
543  vdot = vdot * (A * B);
544  if (vdot < BREP_CURVE_FLATNESS)
545  return false; //already failed
546  }
547 
548  return vdot >= BREP_CURVE_FLATNESS;
549 }
550 
551 
552 //--------------------------------------------------------------------------------
553 // SurfaceTree
554 SurfaceTree::SurfaceTree()
555  : m_removeTrimmed(false),
556  ctree(NULL),
557  m_face(NULL),
558  m_root(NULL)
559 {
560 }
561 
562 
563 SurfaceTree::SurfaceTree(const ON_BrepFace* face, bool removeTrimmed, int depthLimit, double within_distance_tol)
564  : m_removeTrimmed(removeTrimmed),
565  m_face(face)
566 {
567  // build the surface bounding volume hierarchy
568  const ON_Surface* surf = face->SurfaceOf();
569  if (!surf) {
570  TRACE("ERROR: NULL surface encountered in SurfaceTree()");
571  return;
572  }
573 
574  // may be a smaller trimmed subset of surface so worth getting
575  // face boundary
576  bool bGrowBox = false;
577  ON_3dPoint min = ON_3dPoint::UnsetPoint, max = ON_3dPoint::UnsetPoint;
578  for (int li = 0; li < face->LoopCount(); li++) {
579  for (int ti = 0; ti < face->Loop(li)->TrimCount(); ti++) {
580  ON_BrepTrim *trim = face->Loop(li)->Trim(ti);
581  trim->GetBoundingBox(min, max, bGrowBox);
582  bGrowBox = true;
583  }
584  }
585  if (!bGrowBox) {
586  surf->GetBoundingBox(min, max);
587  removeTrimmed = false;
588  }
589 
590  // first, build the Curve Tree
591  if (removeTrimmed)
592  ctree = new CurveTree(m_face);
593  else
594  ctree = NULL;
595 
596  TRACE("Creating surface tree for: " << face->m_face_index);
597 
598 #ifdef _OLD_SUBDIVISION_
599  ON_Interval u = surf->Domain(0);
600  ON_Interval v = surf->Domain(1);
601 #else
602  ON_Interval dom[2] = { ON_Interval::EmptyInterval, ON_Interval::EmptyInterval };
603  for (int i = 0; i < 2; i++) {
604  dom[i] = surf->Domain(i);
605 #ifdef LOOSEN_UV
606  min[i] -= within_distance_tol;
607  max[i] += within_distance_tol;
608 #endif
609  if ((min != ON_3dPoint::UnsetPoint) && (max != ON_3dPoint::UnsetPoint)) {
610  if ((min[i] >= dom[i].m_t[0]) && (max[i] <= dom[i].m_t[1])) {
611  dom[i].Set(min[i],max[i]);
612  }
613  }
614  }
615  ON_Interval u = dom[0];
616  ON_Interval v = dom[1];
617 #endif
618  double uq = u.Length()*0.25;
619  double vq = v.Length()*0.25;
620 
621  ///////////////////////////////////////////////////////////////////////
622  // Populate initial frames array for use in tree build
623  ON_Plane frames[9];
624  surf->FrameAt(u.Min(), v.Min(), frames[0]);
625  surf->FrameAt(u.Max(), v.Min(), frames[1]);
626  surf->FrameAt(u.Max(), v.Max(), frames[2]);
627  surf->FrameAt(u.Min(), v.Max(), frames[3]);
628  surf->FrameAt(u.Mid(), v.Mid(), frames[4]);
629  surf->FrameAt(u.Mid() - uq, v.Mid() - vq, frames[5]);
630  surf->FrameAt(u.Mid() - uq, v.Mid() + vq, frames[6]);
631  surf->FrameAt(u.Mid() + uq, v.Mid() - vq, frames[7]);
632  surf->FrameAt(u.Mid() + uq, v.Mid() + vq, frames[8]);
633 
634  m_root = subdivideSurface(surf, u, v, frames, 0, depthLimit, 1, within_distance_tol);
635 
636  if (m_root) {
637  m_root->BuildBBox();
638  }
639  TRACE("u: [" << u[0] << ", " << u[1] << "]");
640  TRACE("v: [" << v[0] << ", " << v[1] << "]");
641  TRACE("m_root: " << m_root);
642  while (!f_queue.empty()) {
643  bu_free(f_queue.front(), "free subsurface frames array");
644  f_queue.pop();
645  }
646 }
647 
648 
649 SurfaceTree::~SurfaceTree()
650 {
651  delete ctree;
652  delete m_root;
653 }
654 
655 
656 BBNode*
657 SurfaceTree::getRootNode() const
658 {
659  return m_root;
660 }
661 
662 
663 int
664 SurfaceTree::depth()
665 {
666  return m_root->depth();
667 }
668 
669 
670 ON_2dPoint
671 SurfaceTree::getClosestPointEstimate(const ON_3dPoint& pt)
672 {
673  return m_root->getClosestPointEstimate(pt);
674 }
675 
676 
677 ON_2dPoint
678 SurfaceTree::getClosestPointEstimate(const ON_3dPoint& pt, ON_Interval& u, ON_Interval& v)
679 {
680  return m_root->getClosestPointEstimate(pt, u, v);
681 }
682 
683 
684 void
685 SurfaceTree::getLeaves(std::list<BBNode*>& out_leaves)
686 {
687  m_root->getLeaves(out_leaves);
688 }
689 
690 
691 const ON_Surface *
692 SurfaceTree::getSurface()
693 {
694  return m_face->SurfaceOf();
695 }
696 
697 
698 int
699 brep_getSurfacePoint(const ON_3dPoint& pt, ON_2dPoint& uv , BBNode* node) {
700  plane_ray pr;
701  const ON_Surface *surf = node->m_face->SurfaceOf();
702  double umin, umax;
703  double vmin, vmax;
704  surf->GetDomain(0, &umin, &umax);
705  surf->GetDomain(1, &vmin, &vmax);
706 
707  ON_3dVector dir = node->m_normal;
708  dir.Reverse();
709  ON_Ray ray((ON_3dPoint&)pt, dir);
710  brep_get_plane_ray(ray, pr);
711 
712  //know use this as guess to iterate to closer solution
713  pt2d_t Rcurr;
714  pt2d_t new_uv;
715  ON_3dVector su, sv;
716  bool found=false;
717  fastf_t Dlast = MAX_FASTF;
718  pt2d_t nuv;
719  nuv[0] = (node->m_u[1] + node->m_u[0])/2.0;
720  nuv[1] = (node->m_v[1] + node->m_v[0])/2.0;
721  ON_3dPoint newpt;
722  for (int i = 0; i < BREP_MAX_ITERATIONS; i++) {
723  brep_r(surf, pr, nuv, newpt, su, sv, Rcurr);
724  fastf_t d = v2mag(Rcurr);
725 
727  TRACE1("R:"<<ON_PRINT2(Rcurr));
728  found = true;
729  break;
730  } else if (d < BREP_INTERSECTION_ROOT_SETTLE) {
731  found = true;
732  }
733  brep_newton_iterate(pr, Rcurr, su, sv, nuv, new_uv);
734 
735  //Check for closed surface wrap around
736  if (surf->IsClosed(0)) {
737  CLAMP(new_uv[0], umin, umax);
738  }
739  if (surf->IsClosed(1)) {
740  CLAMP(new_uv[1], vmin, vmax);
741  }
742 
743 #ifdef HOOD
744  //push answer back to within node bounds
745  double ufluff = (node->m_u[1] - node->m_u[0])*0.01;
746  double vfluff = (node->m_v[1] - node->m_v[0])*0.01;
747 #else
748  //push answer back to within node bounds
749  double ufluff = 0.0;
750  double vfluff = 0.0;
751 #endif
752 
753  if (new_uv[0] < node->m_u[0] - ufluff)
754  new_uv[0] = node->m_u[0];
755  else if (new_uv[0] > node->m_u[1] + ufluff)
756  new_uv[0] = node->m_u[1];
757 
758  if (new_uv[1] < node->m_v[0] - vfluff)
759  new_uv[1] = node->m_v[0];
760  else if (new_uv[1] > node->m_v[1] + vfluff)
761  new_uv[1] = node->m_v[1];
762 
763 
764  surface_EvNormal(surf,new_uv[0], new_uv[1], newpt, ray.m_dir);
765  ray.m_dir.Reverse();
766  brep_get_plane_ray(ray, pr);
767 
768  if (d < Dlast) {
769  move(nuv, new_uv);
770  Dlast = d;
771  }
772  }
773  if (found) {
774  uv.x = nuv[0];
775  uv.y = nuv[1];
776  return 1;
777  }
778  return -1;
779 }
780 
781 
782 int
783 SurfaceTree::getSurfacePoint(const ON_3dPoint& pt, ON_2dPoint& uv, const ON_3dPoint& from, double tolerance) const
784 {
785  std::list<BBNode*> nodes;
786  (void)m_root->getLeavesBoundingPoint(from, nodes);
787 
788  double min_dist = MAX_FASTF;
789  ON_2dPoint curr_uv(0.0, 0.0);
790  bool found = false;
791 
792  std::list<BBNode*>::iterator i;
793  for (i = nodes.begin(); i != nodes.end(); i++) {
794  BBNode* node = (*i);
795  if (brep_getSurfacePoint(pt, curr_uv, node)) {
796  ON_3dPoint fp = m_face->SurfaceOf()->PointAt(curr_uv.x, curr_uv.y);
797  double dist = fp.DistanceTo(pt);
798  if (NEAR_ZERO(dist, BREP_SAME_POINT_TOLERANCE)) {
799  uv = curr_uv;
800  found = true;
801  return 1; //close enough to same point so no sense in looking for one closer
802  } else if (NEAR_ZERO(dist, tolerance)) {
803  if (dist < min_dist) {
804  uv = curr_uv;
805  min_dist = dist;
806  found = true; //within tolerance but may be a point closer so keep looking
807  }
808  }
809  }
810  }
811  if (found) {
812  return 1;
813  }
814 
815  nodes.clear();
816  (void)m_root->getLeavesBoundingPoint(pt, nodes);
817  for (i = nodes.begin(); i != nodes.end(); i++) {
818  BBNode* node = (*i);
819  if (brep_getSurfacePoint(pt, curr_uv, node)) {
820  ON_3dPoint fp = m_face->SurfaceOf()->PointAt(curr_uv.x, curr_uv.y);
821  double dist = fp.DistanceTo(pt);
822  if (NEAR_ZERO(dist, BREP_SAME_POINT_TOLERANCE)) {
823  uv = curr_uv;
824  found = true;
825  return 1; //close enough to same point so no sense in looking for one closer
826  } else if (NEAR_ZERO(dist, tolerance)) {
827  if (dist < min_dist) {
828  uv = curr_uv;
829  min_dist = dist;
830  found = true; //within tolerance but may be a point closer so keep looking
831  }
832  }
833  }
834  }
835 
836  return -1;
837 }
838 
839 
840 //static int bb_cnt=0;
841 BBNode*
842 SurfaceTree::surfaceBBox(const ON_Surface *localsurf,
843  bool isLeaf,
844  ON_Plane *m_frames,
845  const ON_Interval& u,
846  const ON_Interval& v,
847  double within_distance_tol)
848 {
849  point_t min, max, buffer;
850 #ifdef _OLD_SUBDIVISION_
851  ON_BoundingBox bbox = localsurf->BoundingBox();
852 #else
853  ON_BoundingBox bbox = ON_BoundingBox::EmptyBoundingBox;
854  if (!surface_GetBoundingBox(localsurf,u,v,bbox,false)) {
855  return NULL;
856  }
857 #endif
858 
859  VSETALL(buffer, within_distance_tol);
860 
861  //bu_log("in bb%d rpp %f %f %f %f %f %f\n", bb_cnt, min[0], max[0], min[1], max[1], min[2], max[2]);
862  VMOVE(min, bbox.Min());
863  VMOVE(max, bbox.Max());
864 
865  // calculate the estimate point on the surface: i.e. use the point
866  // on the surface defined by (u.Mid(), v.Mid()) as a heuristic for
867  // finding the uv domain bounding a portion of the surface close
868  // to an arbitrary model-space point (see
869  // getClosestPointEstimate())
870  ON_3dPoint estimate;
871  ON_3dVector normal;
872  estimate = m_frames[4].origin;
873  normal = m_frames[4].zaxis;
874 
875  BBNode* node;
876  if (isLeaf) {
877  /* vect_t delta; */
878 
879  VSUB2(min, min, buffer);
880  VADD2(max, max, buffer);
881  /* VSUB2(delta, max, min);
882  VSCALE(delta, delta, BBOX_GROW_3D);
883  VSUB2(min, min, delta);
884  VADD2(max, max, delta);
885  */
886  TRACE("creating leaf: u(" << u.Min() << ", " << u.Max() <<
887  ") v(" << v.Min() << ", " << v.Max() << ")");
888  node = new BBNode(ctree, ON_BoundingBox(ON_3dPoint(min),
889  ON_3dPoint(max)),
890  m_face,
891  u, v);
892  node->prepTrims();
893 
894  } else {
895  node = new BBNode(ctree, ON_BoundingBox(ON_3dPoint(min),
896  ON_3dPoint(max)));
897  }
898 
899  node->m_estimate = estimate;
900  node->m_normal = normal;
901  node->m_face = m_face;
902  node->m_u = u;
903  node->m_v = v;
904  return node;
905 }
906 
907 
908 BBNode*
909 initialBBox(CurveTree* ctree, const ON_Surface* surf, const ON_BrepFace* face, const ON_Interval& u, const ON_Interval& v)
910 {
911 #ifdef _OLD_SUBDIVISION_
912  ON_BoundingBox bb = surf->BoundingBox();
913 #else
914  ON_BoundingBox bb = ON_BoundingBox::EmptyBoundingBox;
915  if (!surface_GetBoundingBox(surf,u,v,bb,false)) {
916  return NULL;
917  }
918 #endif
919  BBNode* node = new BBNode(ctree, bb, face, u, v, false, false);
920  ON_3dPoint estimate;
921  ON_3dVector normal;
922  if (!surface_EvNormal(surf,surf->Domain(0).Mid(), surf->Domain(1).Mid(), estimate, normal)) {
923  bu_bomb("Could not evaluate estimate point on surface");
924  }
925  node->m_estimate = estimate;
926  node->m_normal = normal;
927  node->m_face = face;
928  node->m_u = u;
929  node->m_v = v;
930  return node;
931 }
932 
933 
934 // Cache surface information as file static to ensure initialization once;
935 static const ON_Surface *prev_surf[MAX_PSW] = {NULL};
936 static ON_Interval dom[MAX_PSW][2] = {{ON_Interval::EmptyInterval, ON_Interval::EmptyInterval}};
937 static int span_cnt[MAX_PSW][2] = {{0, 0}};
938 static double *span[MAX_PSW][2] = {{NULL, NULL}};
939 bool
940 hasSplit(const ON_Surface *surf, const int dir,const ON_Interval& interval,double &split)
941 {
942  int p = bu_parallel_id();
943  if (surf == NULL) {
944  // clean up statics and return
945  prev_surf[p] = NULL;
946  for(int i=0; i<2; i++) {
947  span_cnt[p][i] = 0;
948  if (span[p][i])
949  delete [] span[p][i];
950  span[p][i] = NULL;
951  dom[p][i] = ON_Interval::EmptyInterval;
952  }
953 
954  return false;
955  }
956  if (prev_surf[p] != surf ) {
957  // load new surf info
958  for(int i=0; i<2; i++) {
959  if (span[p][i])
960  delete [] span[p][i];
961  dom[p][i] = surf->Domain(i);
962  span_cnt[p][i] = surf->SpanCount(i);
963  span[p][i] = new double[span_cnt[p][i]+1];
964  surf->GetSpanVector(i, span[p][i]);
965  }
966 
967  prev_surf[p] = surf;
968  }
969 
970  // find direction split based on setting of 'dir'
971 
972  // first, if closed in 'dir' check to see if it extends over seam and use that as split
973  if (surf->IsClosed(dir)) {
974  bool testOpen = true;
975  if (interval.Includes(dom[p][dir].m_t[0], testOpen)) { //crosses lower boundary
976  split = dom[p][dir].m_t[0];
977  return true;
978  } else if (interval.Includes(dom[p][dir].m_t[1], testOpen)) { //crosses upper boundary
979  split = dom[p][dir].m_t[1];
980  return true;
981  }
982  }
983 
984  // next lets see if we have a knots in interval, if so split on middle knot
985  if (span_cnt[p][dir] > 1) {
986  int sum = 0;
987  int cnt = 0;
988  for(int i=0; i<span_cnt[p][i]+1; i++) {
989  bool testOpen = true;
990  if (interval.Includes(span[p][dir][i], testOpen)) { //crosses lower boundary
991  sum = sum + i;
992  cnt++;
993  }
994  }
995  if (cnt > 0) {
996  split = span[p][dir][sum/cnt];
997  return true;
998  }
999  }
1000 
1001  return false;
1002 }
1003 
1004 
1005 BBNode*
1006 SurfaceTree::subdivideSurface(const ON_Surface *localsurf,
1007  const ON_Interval& u,
1008  const ON_Interval& v,
1009  ON_Plane frames[],
1010  int divDepth,
1011  int depthLimit,
1012  int prev_knot,
1013  double within_distance_tol
1014  )
1015 {
1016  BBNode* quads[4];
1017  BBNode* parent = NULL;
1018  double usplit = 0.0;
1019  double vsplit = 0.0;
1020  double width = 0.0;
1021  double height = 0.0;
1022  double ratio = 5.0;
1023  double uq = u.Length()*0.25;
1024  double vq = v.Length()*0.25;
1025  localsurf->FrameAt(u.Mid() - uq, v.Mid() - vq, frames[5]);
1026  localsurf->FrameAt(u.Mid() - uq, v.Mid() + vq, frames[6]);
1027  localsurf->FrameAt(u.Mid() + uq, v.Mid() - vq, frames[7]);
1028  localsurf->FrameAt(u.Mid() + uq, v.Mid() + vq, frames[8]);
1029  unsigned int do_u_split = 0;
1030  unsigned int do_v_split = 0;
1031  unsigned int do_both_splits = 0;
1032 
1033  usplit = u.Mid();
1034  vsplit = v.Mid();
1035 
1036 #ifndef _OLD_SUBDIVISION_
1037  if (divDepth >= depthLimit) {
1038  return surfaceBBox(localsurf, true, frames, u, v, within_distance_tol);
1039  }
1040 #endif
1041 
1042  // The non-knot case where all criteria are satisfied is the
1043  // terminating case for the recursion - handle that first
1044  if (!prev_knot) {
1045 #ifdef _OLD_SUBDIVISION_
1046  localsurf->GetSurfaceSize(&width, &height);
1047 #else
1048  width = u.Length();
1049  height = v.Length();
1050 #endif
1051  if (((width/height < ratio) && (width/height > 1.0/ratio) && isFlat(frames) && isStraight(frames))
1052  || (divDepth >= depthLimit)) { //BREP_MAX_FT_DEPTH
1053  return surfaceBBox(localsurf, true, frames, u, v, within_distance_tol);
1054  }
1055  }
1056 
1057  // Knots
1058  if (prev_knot) {
1059 #ifdef _OLD_SUBDIVISION_
1060  int spanu_cnt = localsurf->SpanCount(0);
1061  int spanv_cnt = localsurf->SpanCount(1);
1062  parent = initialBBox(ctree, localsurf, m_face, u, v);
1063  if (spanu_cnt > 1) {
1064  double *spanu = new double[spanu_cnt+1];
1065  localsurf->GetSpanVector(0, spanu);
1066  do_u_split = 1;
1067  usplit = spanu[(spanu_cnt+1)/2];
1068  delete [] spanu;
1069  }
1070  if (spanv_cnt > 1) {
1071  double *spanv = new double[spanv_cnt+1];
1072  localsurf->GetSpanVector(1, spanv);
1073  do_v_split = 1;
1074  vsplit = spanv[(spanv_cnt+1)/2];
1075  delete [] spanv;
1076  }
1077  if (do_u_split && do_v_split) {
1078  do_both_splits = 1;
1079  do_u_split = 0;
1080  do_v_split = 0;
1081  }
1082 #else
1083  int dir = 0; // U direction
1084  if (hasSplit(localsurf,dir,u,usplit)) {
1085  do_u_split = 1;
1086  }
1087  dir = 1; // V direction
1088  if (hasSplit(localsurf,dir,v,vsplit)) {
1089  if (do_u_split) {
1090  do_both_splits = 1;
1091  do_u_split = 0;
1092  } else {
1093  do_v_split = 1;
1094  }
1095  }
1096  parent = initialBBox(ctree, localsurf, m_face, u, v);
1097 #endif
1098  }
1099  // Flatness
1100  if (!prev_knot) {
1101  bool isUFlat = isFlatU(frames);
1102  bool isVFlat = isFlatV(frames);
1103 
1104  parent = (divDepth == 0) ? initialBBox(ctree, localsurf, m_face, u, v) : surfaceBBox(localsurf, false, frames, u, v, within_distance_tol);
1105 
1106  if ((!isVFlat || (width/height > ratio)) && (!isUFlat || (height/width > ratio))) {
1107  do_both_splits = 1;
1108  }
1109  if ((!isUFlat || (width/height > ratio)) && !do_both_splits) {
1110  do_u_split = 1;
1111  }
1112  if (!do_both_splits && ! do_u_split) {
1113  do_v_split = 1;
1114  }
1115  }
1116 
1117  ///////////////////////////////////
1118 
1119  if (do_both_splits) {
1120  ON_Interval firstu(u.Min(), usplit);
1121  ON_Interval secondu(usplit, u.Max());
1122  ON_Interval firstv(v.Min(), vsplit);
1123  ON_Interval secondv(vsplit, v.Max());
1124 
1125 #ifdef _OLD_SUBDIVISION_
1126  ON_Surface *q0surf = NULL;
1127  ON_Surface *q1surf = NULL;
1128  ON_Surface *q2surf = NULL;
1129  ON_Surface *q3surf = NULL;
1130 
1131  bool split = ON_Surface_Quad_Split(localsurf, u, v, usplit, vsplit, &q0surf, &q1surf, &q2surf, &q3surf);
1132  /* FIXME: this needs to be handled more gracefully */
1133  if (!split) {
1134  delete parent;
1135  return NULL;
1136  }
1137  q0surf->ClearBoundingBox();
1138  q1surf->ClearBoundingBox();
1139  q2surf->ClearBoundingBox();
1140  q3surf->ClearBoundingBox();
1141 #else
1142  const ON_Surface *q0surf = localsurf;
1143  const ON_Surface *q1surf = localsurf;
1144  const ON_Surface *q2surf = localsurf;
1145  const ON_Surface *q3surf = localsurf;
1146 
1147 #endif
1148 
1149  /*********************************************************************
1150  * In order to avoid fairly expensive re-calculation of 3d points at
1151  * uv coordinates, all values that are shared between children at
1152  * the same depth of a surface subdivision are pre-computed and
1153  * passed as parameters.
1154  *
1155  * The majority of these points are already evaluated in the process
1156  * of testing whether a subdivision has produced a leaf node. These
1157  * values are in the normals and corners arrays and have index values
1158  * corresponding to the values of the figure on the left below. There
1159  * are four other shared values that are precomputed in the sharedcorners
1160  * and sharednormals arrays; their index values in those arrays are
1161  * illustrated in the figure on the right:
1162  *
1163  *
1164  * 3-------------------2 +---------2---------+
1165  * | | | |
1166  * | 6 8 | | |
1167  * | | | |
1168  * V| 4 | 1 3
1169  * | | | |
1170  * | 5 7 | | |
1171  * | | | |
1172  * 0-------------------1 +---------0---------+
1173  * U U
1174  *
1175  * Values inherited from Values pre-prepared in
1176  * parent subdivision shared arrays
1177  *
1178  *
1179  * When the four subdivisions are made, the parent parameters
1180  * are passed to the children as follows (values from the
1181  * shared arrays are prefaced with an s):
1182  *
1183  * 3--------------S2 S2--------------2
1184  * | | | |
1185  * | | | |
1186  * V | 6 | | 8 |
1187  * | | | |
1188  * | | | |
1189  * S1--------------4 4--------------S3
1190  * U U
1191  *
1192  * Quadrant 3 Quadrant 2
1193  *
1194  * S1--------------4 4--------------S3
1195  * | | | |
1196  * | | | |
1197  * V | 5 | | 7 |
1198  * | | | |
1199  * | | | |
1200  * 0--------------S0 S0--------------1
1201  * U U
1202  *
1203  * Quadrant 0 Quadrant 1
1204  *
1205  *
1206  **********************************************************************/
1207 
1208  ON_Plane sharedframes[4];
1209  localsurf->FrameAt(usplit, v.Min(), sharedframes[0]);
1210  localsurf->FrameAt(u.Min(), vsplit, sharedframes[1]);
1211  localsurf->FrameAt(usplit, v.Max(), sharedframes[2]);
1212  localsurf->FrameAt(u.Max(), vsplit, sharedframes[3]);
1213  // When splitting via knots, we don't know what point frames[4] is until
1214  // the knot is selected
1215  if (prev_knot) localsurf->FrameAt(usplit, vsplit, frames[4]);
1216 
1217  ON_Plane *newframes;
1218  if (!f_queue.empty()) {newframes = f_queue.front(); f_queue.pop();} else {newframes = (ON_Plane *)bu_malloc(9*sizeof(ON_Plane), "new frames");}
1219  newframes[0] = frames[0];
1220  newframes[1] = sharedframes[0];
1221  newframes[2] = frames[4];
1222  newframes[3] = sharedframes[1];
1223  newframes[4] = frames[5];
1224  quads[0] = subdivideSurface(q0surf, firstu, firstv, newframes, divDepth+1, depthLimit, prev_knot, within_distance_tol);
1225 #ifdef _OLD_SUBDIVISION_
1226  delete q0surf;
1227 #endif
1228  newframes[0] = sharedframes[0];
1229  newframes[1] = frames[1];
1230  newframes[2] = sharedframes[3];
1231  newframes[3] = frames[4];
1232  newframes[4] = frames[7];
1233  quads[1] = subdivideSurface(q1surf, secondu, firstv, newframes, divDepth+1, depthLimit, prev_knot, within_distance_tol);
1234 #ifdef _OLD_SUBDIVISION_
1235  delete q1surf;
1236 #endif
1237  newframes[0] = frames[4];
1238  newframes[1] = sharedframes[3];
1239  newframes[2] = frames[2];
1240  newframes[3] = sharedframes[2];
1241  newframes[4] = frames[8];
1242  quads[2] = subdivideSurface(q2surf, secondu, secondv, newframes, divDepth+1, depthLimit, prev_knot, within_distance_tol);
1243 #ifdef _OLD_SUBDIVISION_
1244  delete q2surf;
1245 #endif
1246  newframes[0] = sharedframes[1];
1247  newframes[1] = frames[4];
1248  newframes[2] = sharedframes[2];
1249  newframes[3] = frames[3];
1250  newframes[4] = frames[6];
1251  quads[3] = subdivideSurface(q3surf, firstu, secondv, newframes, divDepth+1, depthLimit, prev_knot, within_distance_tol);
1252 #ifdef _OLD_SUBDIVISION_
1253  delete q3surf;
1254 #endif
1255  memset(newframes, 0, 9 * sizeof(ON_Plane *));
1256  f_queue.push(newframes);
1257 
1258  parent->m_trimmed = true;
1259  parent->m_checkTrim = false;
1260 
1261  for (int i = 0; i < 4; i++) {
1262  if (!quads[i]) {
1263  continue;
1264  }
1265  if (!(quads[i]->m_trimmed)) {
1266  parent->m_trimmed = false;
1267  }
1268  if (quads[i]->m_checkTrim) {
1269  parent->m_checkTrim = true;
1270  }
1271  }
1272  if (m_removeTrimmed) {
1273  for (int i = 0; i < 4; i++) {
1274  if (!quads[i]) {
1275  continue;
1276  }
1277  if (!(quads[i]->m_trimmed)) {
1278  parent->addChild(quads[i]);
1279  }
1280  }
1281  } else {
1282  for (int i = 0; i < 4; i++) {
1283  parent->addChild(quads[i]);
1284  }
1285  }
1286  parent->BuildBBox();
1287  return parent;
1288  }
1289  //////////////////////////////////////
1290  if (do_u_split) {
1291  bool split;
1292  ON_Interval firstu(u.Min(), usplit);
1293  ON_Interval secondu(usplit, u.Max());
1294 #ifdef _OLD_SUBDIVISION_
1295  ON_Surface *east = NULL;
1296  ON_Surface *west = NULL;
1297 
1298  int dir = 0;
1299  if (prev_knot) {
1300  split = localsurf->Split(dir, usplit, east, west);
1301  } else {
1302  split = localsurf->Split(dir, localsurf->Domain(dir).Mid(), east, west);
1303  }
1304 
1305  /* FIXME: this needs to be handled more gracefully */
1306  if (!split || !east || !west) {
1307  bu_log("DEBUG: Split failure (split:%d, east:%p, west:%p)\n", split, (void *)east, (void *)west);
1308  delete parent;
1309  return NULL;
1310  }
1311 
1312  east->ClearBoundingBox();
1313  west->ClearBoundingBox();
1314 #else
1315  const ON_Surface *east = localsurf;
1316  const ON_Surface *west = localsurf;
1317  split = true;
1318 #endif
1319 
1320  //////////////////////////////////
1321  /*********************************************************************
1322  * In order to avoid fairly expensive re-calculation of 3d points at
1323  * uv coordinates, all values that are shared between children at
1324  * the same depth of a surface subdivision are pre-computed and
1325  * passed as parameters.
1326  *
1327  * The majority of these points are already evaluated in the process
1328  * of testing whether a subdivision has produced a leaf node. These
1329  * values are in the normals and corners arrays and have index values
1330  * corresponding to the values of the figure on the left below. There
1331  * are four other shared values that are precomputed in the sharedcorners
1332  * and sharednormals arrays; their index values in those arrays are
1333  * illustrated in the figure on the right:
1334  *
1335  *
1336  * 3-------------------2 +---------2---------+
1337  * | | | |
1338  * | 6 8 | | |
1339  * | | | |
1340  * V| 4 | 1 3
1341  * | | | |
1342  * | 5 7 | | |
1343  * | | | |
1344  * 0-------------------1 +---------0---------+
1345  * U U
1346  *
1347  * Values inherited from Values pre-prepared in
1348  * parent subdivision shared arrays
1349  *
1350  *
1351  * When the four subdivisions are made, the parent parameters
1352  * are passed to the children as follows (values from the
1353  * shared arrays are prefaced with an s):
1354  *
1355  * 3--------------S1 S1--------------2
1356  * | | | |
1357  * | | | |
1358  * V | 4 | | 5 |
1359  * | | | |
1360  * | | | |
1361  * 0--------------S0 S0--------------1
1362  * U U
1363  *
1364  * East West
1365  *
1366  *
1367  **********************************************************************/
1368  /* When the four subdivisions are made, the parent parameters
1369  * are passed to the children as follows (values from the
1370  * shared arrays are prefaced with an s):
1371  *
1372  * 3--------------S2 S2--------------2
1373  * | | | |
1374  * | | | |
1375  * V | 6 | | 8 |
1376  * | | | |
1377  * | | | |
1378  * S1--------------4 4--------------S3
1379  * U U
1380  *
1381  * Quadrant 3 Quadrant 2
1382  *
1383  * S1--------------4 4--------------S3
1384  * | | | |
1385  * | | | |
1386  * V | 5 | | 7 |
1387  * | | | |
1388  * | | | |
1389  * 0--------------S0 S0--------------1
1390  * U U
1391  *
1392  * Quadrant 0 Quadrant 1
1393  *
1394  *
1395  **********************************************************************/
1396 
1397  ON_Plane sharedframes[4];
1398  localsurf->FrameAt(usplit, v.Min(), sharedframes[0]);
1399  localsurf->FrameAt(usplit, v.Max(), sharedframes[1]);
1400 
1401  ON_Plane *newframes;
1402  if (!f_queue.empty()) {newframes = f_queue.front(); f_queue.pop();} else {newframes = (ON_Plane *)bu_malloc(9*sizeof(ON_Plane), "new frames");}
1403  newframes[0] = frames[0];
1404  newframes[1] = sharedframes[0];
1405  newframes[2] = sharedframes[1];
1406  newframes[3] = frames[3];
1407 
1408  if (prev_knot) {
1409  localsurf->FrameAt(firstu.Mid(), v.Mid(), newframes[4]);
1410  quads[0] = subdivideSurface(east, firstu, v, newframes, divDepth+1, depthLimit, prev_knot, within_distance_tol);
1411  } else {
1412  ON_Interval first(0, 0.5);
1413  localsurf->FrameAt(u.Mid() - uq, v.Mid(), newframes[4]);
1414  quads[0] = subdivideSurface(east, u.ParameterAt(first), v, newframes, divDepth + 1, depthLimit, prev_knot, within_distance_tol);
1415  }
1416 #ifdef _OLD_SUBDIVISION_
1417  delete east;
1418 #endif
1419 
1420  newframes[0] = sharedframes[0];
1421  newframes[1] = frames[1];
1422  newframes[2] = frames[2];
1423  newframes[3] = sharedframes[1];
1424 
1425  if (prev_knot) {
1426  localsurf->FrameAt(secondu.Mid(), v.Mid(), newframes[4]);
1427  quads[1] = subdivideSurface(west, secondu, v, newframes, divDepth+1, depthLimit, prev_knot, within_distance_tol);
1428  } else {
1429  ON_Interval second(0.5, 1.0);
1430  localsurf->FrameAt(u.Mid() + uq, v.Mid(), newframes[4]);
1431  quads[1] = subdivideSurface(west, u.ParameterAt(second), v, newframes, divDepth + 1, depthLimit, prev_knot, within_distance_tol);
1432  }
1433 #ifdef _OLD_SUBDIVISION_
1434  delete west;
1435 #endif
1436 
1437  memset(newframes, 0, 9 * sizeof(ON_Plane *));
1438  f_queue.push(newframes);
1439 
1440  parent->m_trimmed = true;
1441  parent->m_checkTrim = false;
1442 
1443  for (int i = 0; i < 2; i++) {
1444  if (!quads[i]) {
1445  continue;
1446  }
1447  if (!(quads[i]->m_trimmed)) {
1448  parent->m_trimmed = false;
1449  }
1450  if (quads[i]->m_checkTrim) {
1451  parent->m_checkTrim = true;
1452  }
1453  }
1454  if (m_removeTrimmed) {
1455  for (int i = 0; i < 2; i++) {
1456  if (!quads[i]) {
1457  continue;
1458  }
1459  if (!(quads[i]->m_trimmed)) {
1460  parent->addChild(quads[i]);
1461  }
1462  }
1463  } else {
1464  for (int i = 0; i < 2; i++) {
1465  parent->addChild(quads[i]);
1466  }
1467  }
1468  parent->BuildBBox();
1469  return parent;
1470  }
1471  if (do_v_split || !prev_knot) {
1472  bool split;
1473  ON_Interval firstv(v.Min(), vsplit);
1474  ON_Interval secondv(vsplit, v.Max());
1475 
1476 #ifdef _OLD_SUBDIVISION_
1477  //////////////////////////////////////
1478  ON_Surface *north = NULL;
1479  ON_Surface *south = NULL;
1480 
1481  int dir = 1;
1482  if (prev_knot) {
1483  split = localsurf->Split(dir, vsplit, south, north);
1484  } else {
1485  split = localsurf->Split(dir, localsurf->Domain(dir).Mid(), south, north);
1486  }
1487 
1488  /* FIXME: this needs to be handled more gracefully */
1489  if (!split || !south || !north) {
1490  bu_log("DEBUG: Split failure (split:%d, surf1:%p, surf2:%p)\n", split, (void *)south, (void *)north);
1491  delete parent;
1492  return NULL;
1493  }
1494 
1495  south->ClearBoundingBox();
1496  north->ClearBoundingBox();
1497 #else
1498  const ON_Surface *north = localsurf;
1499  const ON_Surface *south = localsurf;
1500  split = true;
1501 #endif
1502  //////////////////////////////////
1503  /*********************************************************************
1504  * In order to avoid fairly expensive re-calculation of 3d points at
1505  * uv coordinates, all values that are shared between children at
1506  * the same depth of a surface subdivision are pre-computed and
1507  * passed as parameters.
1508  *
1509  * The majority of these points are already evaluated in the process
1510  * of testing whether a subdivision has produced a leaf node. These
1511  * values are in the normals and corners arrays and have index values
1512  * corresponding to the values of the figure on the left below. There
1513  * are four other shared values that are precomputed in the sharedcorners
1514  * and sharednormals arrays; their index values in those arrays are
1515  * illustrated in the figure on the right:
1516  *
1517  *
1518  * 3-------------------2 +---------2---------+
1519  * | | | |
1520  * | 6 8 | | |
1521  * | | | |
1522  * V| 4 | 1 3
1523  * | | | |
1524  * | 5 7 | | |
1525  * | | | |
1526  * 0-------------------1 +---------0---------+
1527  * U U
1528  *
1529  * Values inherited from Values pre-prepared in
1530  * parent subdivision shared arrays
1531  *
1532  *
1533  * When the four subdivisions are made, the parent parameters
1534  * are passed to the children as follows (values from the
1535  * shared arrays are prefaced with an s):
1536  *
1537  * 3--------------------2
1538  * | |
1539  * | |
1540  * V | 5 |
1541  * | |
1542  * | |
1543  * S0-------------------S1
1544  * U
1545  *
1546  * North
1547  *
1548  * S0-------------------S1
1549  * | |
1550  * | |
1551  * V | 4 |
1552  * | |
1553  * | |
1554  * 0--------------------1
1555  * U
1556  *
1557  * South
1558  *
1559  *
1560  **********************************************************************/
1561 
1562  ON_Plane sharedframes[2];
1563  localsurf->FrameAt(u.Min(), vsplit, sharedframes[0]);
1564  localsurf->FrameAt(u.Max(), vsplit, sharedframes[1]);
1565 
1566  ON_Plane *newframes;
1567  if (!f_queue.empty()) {newframes = f_queue.front(); f_queue.pop();} else {newframes = (ON_Plane *)bu_malloc(9*sizeof(ON_Plane), "new frames");}
1568  newframes[0] = frames[0];
1569  newframes[1] = frames[1];
1570  newframes[2] = sharedframes[1];
1571  newframes[3] = sharedframes[0];
1572  if (prev_knot) {
1573  localsurf->FrameAt(u.Mid(), firstv.Mid(), newframes[4]);
1574  quads[0] = subdivideSurface(south, u, firstv, newframes, divDepth+1, depthLimit, prev_knot, within_distance_tol);
1575  } else {
1576  ON_Interval first(0, 0.5);
1577  localsurf->FrameAt(u.Mid(), v.Mid() - vq, newframes[4]);
1578  quads[0] = subdivideSurface(south, u, v.ParameterAt(first), newframes, divDepth + 1, depthLimit, prev_knot, within_distance_tol);
1579  }
1580 #ifdef _OLD_SUBDIVISION_
1581  delete south;
1582 #endif
1583 
1584  newframes[0] = sharedframes[0];
1585  newframes[1] = sharedframes[1];
1586  newframes[2] = frames[2];
1587  newframes[3] = frames[3];
1588 
1589  if (prev_knot) {
1590  localsurf->FrameAt(u.Mid(), secondv.Mid(), newframes[4]);
1591  quads[1] = subdivideSurface(north, u, secondv, newframes, divDepth+1, depthLimit, prev_knot, within_distance_tol);
1592  } else {
1593  ON_Interval second(0.5, 1.0);
1594  localsurf->FrameAt(u.Mid(), v.Mid() + vq, newframes[4]);
1595  quads[1] = subdivideSurface(north, u, v.ParameterAt(second), newframes, divDepth + 1, depthLimit, prev_knot, within_distance_tol);
1596  }
1597 #ifdef _OLD_SUBDIVISION_
1598  delete north;
1599 #endif
1600 
1601  memset(newframes, 0, 9 * sizeof(ON_Plane *));
1602  f_queue.push(newframes);
1603 
1604  parent->m_trimmed = true;
1605  parent->m_checkTrim = false;
1606 
1607  for (int i = 0; i < 2; i++) {
1608  if (!quads[i]) {
1609  continue;
1610  }
1611  if (!(quads[i]->m_trimmed)) {
1612  parent->m_trimmed = false;
1613  }
1614  if (quads[i]->m_checkTrim) {
1615  parent->m_checkTrim = true;
1616  }
1617  }
1618  if (m_removeTrimmed) {
1619  for (int i = 0; i < 2; i++) {
1620  if (!quads[i]) {
1621  continue;
1622  }
1623  if (!(quads[i]->m_trimmed)) {
1624  parent->addChild(quads[i]);
1625  }
1626  }
1627  } else {
1628  for (int i = 0; i < 2; i++) {
1629  parent->addChild(quads[i]);
1630  }
1631  }
1632 
1633  parent->BuildBBox();
1634  return parent;
1635  }
1636 
1637  if (!do_both_splits && !do_u_split && !do_v_split) {
1638  ((ON_Surface *)localsurf)->ClearBoundingBox();
1639  delete parent;
1640  return subdivideSurface(localsurf, u, v, frames, 0, depthLimit, 0, within_distance_tol);
1641  }
1642 
1643  // Should never get here
1644  return NULL;
1645 }
1646 
1647 
1648 bool
1649 SurfaceTree::isFlat(ON_Plane *frames)
1650 {
1651  return ON_Surface_IsFlat(frames, BREP_SURFACE_FLATNESS);
1652 }
1653 
1654 
1655 bool
1656 SurfaceTree::isStraight(ON_Plane *frames)
1657 {
1658  return ON_Surface_IsStraight(frames, BREP_SURFACE_FLATNESS);
1659 }
1660 
1661 
1662 bool
1663 SurfaceTree::isFlatU(ON_Plane *frames)
1664 {
1665  return ON_Surface_IsFlat_U(frames, BREP_SURFACE_FLATNESS);
1666 }
1667 
1668 
1669 bool
1670 SurfaceTree::isFlatV(ON_Plane *frames)
1671 {
1672  return ON_Surface_IsFlat_V(frames, BREP_SURFACE_FLATNESS);
1673 }
1674 
1675 
1676 //--------------------------------------------------------------------------------
1677 // get_closest_point implementation
1678 
1679 typedef struct _gcp_data {
1680  const ON_Surface* surf;
1681  ON_3dPoint pt;
1682 
1683  // for newton iteration
1684  ON_3dPoint S;
1685  ON_3dVector du;
1686  ON_3dVector dv;
1687  ON_3dVector duu;
1688  ON_3dVector dvv;
1689  ON_3dVector duv;
1690 } GCPData;
1691 
1692 
1693 bool
1694 gcp_gradient(pt2d_t out_grad, GCPData& data, pt2d_t uv)
1695 {
1696  bool evaluated = data.surf->Ev2Der(uv[0],
1697  uv[1],
1698  data.S,
1699  data.du,
1700  data.dv,
1701  data.duu,
1702  data.duv,
1703  data.dvv); // calc S(u, v) dS/du dS/dv d2S/du2 d2S/dv2 d2S/dudv
1704  if (!evaluated) return false;
1705  out_grad[0] = 2 * (data.du * (data.S - data.pt));
1706  out_grad[1] = 2 * (data.dv * (data.S - data.pt));
1707  return true;
1708 }
1709 
1710 
1711 bool
1712 gcp_newton_iteration(pt2d_t out_uv, GCPData& data, pt2d_t grad, pt2d_t in_uv)
1713 {
1714  ON_3dVector delta = data.S - data.pt;
1715  double g1du = 2 * (data.duu * delta) + 2 * (data.du * data.du);
1716  double g2du = 2 * (data.duv * delta) + 2 * (data.dv * data.du);
1717  double g1dv = g2du;
1718  double g2dv = 2 * (data.dvv * delta) + 2 * (data.dv * data.dv);
1719  mat2d_t jacob = { g1du, g1dv,
1720  g2du, g2dv };
1721  mat2d_t inv_jacob;
1722  if (mat2d_inverse(inv_jacob, jacob)) {
1723  pt2d_t tmp;
1724  mat2d_pt2d_mul(tmp, inv_jacob, grad);
1725  pt2dsub(out_uv, in_uv, tmp);
1726  return true;
1727  } else {
1728  std::cerr << "inverse failed!" << std::endl; // XXX fix the error handling
1729  return false;
1730  }
1731 }
1732 
1733 
1734 bool
1735 get_closest_point(ON_2dPoint& outpt,
1736  ON_BrepFace* face,
1737  const ON_3dPoint& point,
1738  SurfaceTree* tree,
1739  double tolerance)
1740 {
1741  int try_count = 0;
1742  bool delete_tree = false;
1743  bool found = false;
1744  double d_last = DBL_MAX;
1745  pt2d_t curr_grad = {0.0, 0.0};
1746  pt2d_t new_uv = {0.0, 0.0};
1747  GCPData data;
1748  data.surf = face->SurfaceOf();
1749  data.pt = point;
1750 
1751  TRACE("get_closest_point: " << PT(point));
1752 
1753  // get initial estimate
1754  SurfaceTree* a_tree = tree;
1755  if (a_tree == NULL) {
1756  a_tree = new SurfaceTree(face);
1757  delete_tree = true;
1758  }
1759  ON_Interval u, v;
1760  ON_2dPoint est = a_tree->getClosestPointEstimate(point, u, v);
1761  pt2d_t uv = {est[0], est[1]};
1762 
1763  // do the newton iterations
1764  // terminates on 1 of 3 conditions:
1765  // 1. if the gradient falls below an epsilon (preferred :-)
1766  // 2. if the gradient diverges
1767  // 3. iterated MAX_FCP_ITERATIONS
1768 try_again:
1769  int diverge_count = 0;
1770  for (int i = 0; i < BREP_MAX_FCP_ITERATIONS; i++) {
1771  assert(gcp_gradient(curr_grad, data, uv));
1772 
1773  ON_3dPoint p = data.surf->PointAt(uv[0], uv[1]);
1774  double d = p.DistanceTo(point);
1775  TRACE("dist: " << d);
1776 
1777  if (NEAR_EQUAL(d, d_last, tolerance)) {
1778  found = true; break;
1779  } else if (d > d_last) {
1780  TRACE("diverged!");
1781  diverge_count++;
1782  }
1783  if (gcp_newton_iteration(new_uv, data, curr_grad, uv)) {
1784  move(uv, new_uv);
1785  } else {
1786  // iteration failed, nudge diagonally
1787  uv[0] += VUNITIZE_TOL;
1788  uv[1] += VUNITIZE_TOL;
1789  }
1790  d_last = d;
1791  }
1792  if (found) {
1793  // check to see if we've left the surface domain
1794  double l, h;
1795  data.surf->GetDomain(0, &l, &h);
1796  CLAMP(uv[0], l, h); // make sure in range
1797 
1798  data.surf->GetDomain(1, &l, &h);
1799  CLAMP(uv[1], l, h);
1800 
1801  outpt[0] = uv[0];
1802  outpt[1] = uv[1];
1803  } else {
1804  TRACE("FAILED TO FIND CLOSEST POINT!");
1805  // XXX: try the mid point of the domain -- HACK! but it seems to work!?
1806  if (try_count == 0) {
1807  uv[0] = u.Mid();
1808  uv[1] = v.Mid();
1809  ++try_count;
1810  goto try_again;
1811  }
1812  }
1813 
1814  if (delete_tree)
1815  delete a_tree;
1816  return found;
1817 }
1818 
1819 
1820 bool
1821 sortX(BRNode* first, BRNode* second)
1822 {
1823  point_t first_min, second_min;
1824  point_t first_max, second_max;
1825 
1826  first->GetBBox(first_min, first_max);
1827  second->GetBBox(second_min, second_max);
1828 
1829  if (first_min[X] < second_min[X])
1830  return true;
1831  else
1832  return false;
1833 }
1834 
1835 
1836 bool
1837 sortY(BRNode* first, BRNode* second)
1838 {
1839  point_t first_min, second_min;
1840  point_t first_max, second_max;
1841 
1842  first->GetBBox(first_min, first_max);
1843  second->GetBBox(second_min, second_max);
1844 
1845  if (first_min[Y] < second_min[Y])
1846  return true;
1847  else
1848  return false;
1849 }
1850 
1851 
1852 } /* end namespace */
1853 
1854 // Local Variables:
1855 // tab-width: 8
1856 // mode: C++
1857 // c-basic-offset: 4
1858 // indent-tabs-mode: t
1859 // c-file-style: "stroustrup"
1860 // End:
1861 // ex: shiftwidth=4 tabstop=8
#define TOL2
another arbitrary calculation tolerance (need to try VDIVIDE_TOL or VUNITIZE_TOL to tighten the bound...
bool sortX(BRNode *first, BRNode *second)
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
bool ON_Surface_IsFlat_U(ON_Plane *frames, double f_tol)
Perform flatness test of surface in U only.
bool hasSplit(const ON_Surface *surf, const int dir, const ON_Interval &interval, double &split)
#define BREP_MAX_ITERATIONS
Definition: brep.h:63
#define VSET(a, b, c, d)
Definition: color.c:53
#define VSETALL(a, s)
Definition: color.c:54
bool gcp_newton_iteration(pt2d_t out_uv, GCPData &data, pt2d_t grad, pt2d_t in_uv)
bool isFlat(const ON_2dPoint &p1, const ON_2dPoint &m, const ON_2dPoint &p2, double flatness)
Header file for the BRL-CAD common definitions.
void brep_r(const ON_Surface *surf, plane_ray &pr, pt2d_t uv, ON_3dPoint &pt, ON_3dVector &su, ON_3dVector &sv, pt2d_t R)
BBNode * initialBBox(CurveTree *ctree, const ON_Surface *surf, const ON_BrepFace *face, const ON_Interval &u, const ON_Interval &v)
#define MAX_FASTF
Definition: defines.h:340
ustring width
bool get_closest_point(ON_2dPoint &outpt, ON_BrepFace *face, const ON_3dPoint &point, SurfaceTree *tree, double tolerance)
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
Definition: color.c:49
COMPLEX data[64]
Definition: fftest.c:34
void * memset(void *s, int c, size_t n)
int brep_getSurfacePoint(const ON_3dPoint &pt, ON_2dPoint &uv, BBNode *node)
#define TOL
arbitrary calculation tolerance (need to try VDIVIDE_TOL or VUNITIZE_TOL to tighten the bounds) ...
#define BREP_INTERSECTION_ROOT_EPSILON
Definition: brep.h:66
void brep_get_plane_ray(ON_Ray &r, plane_ray &pr)
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
bool gcp_gradient(pt2d_t out_grad, GCPData &data, pt2d_t uv)
ON_BOOL32 surface_EvNormal(const ON_Surface *surf, double s, double t, ON_3dPoint &point, ON_3dVector &normal, int side, int *hint)
Coord * point
Definition: chull3d.cpp:52
int bu_parallel_id(void)
Definition: parallel.c:152
ON_BOOL32 surface_GetBoundingBox(const ON_Surface *surf, const ON_Interval &u_interval, const ON_Interval &v_interval, ON_BoundingBox &bbox, ON_BOOL32 bGrowBox)
const ON_Surface * surf
double ON_Curve_Get_Horizontal_Tangent(const ON_Curve *curve, double min, double max, double zero_tol)
Search for a horizontal tangent on the curve between two curve parameters.
void brep_newton_iterate(plane_ray &pr, pt2d_t R, ON_3dVector &su, ON_3dVector &sv, pt2d_t uv, pt2d_t out_uv)
double ON_Curve_Get_Vertical_Tangent(const ON_Curve *curve, double min, double max, double zero_tol)
Search for a vertical tangent on the curve between two curve parameters.
#define MAX_PSW
Definition: parallel.h:48
#define BREP_INTERSECTION_ROOT_SETTLE
Definition: brep.h:69
bool ON_Surface_IsFlat_V(ON_Plane *frames, double f_tol)
Perform flatness test of surface in V only.
bool ON_Surface_IsStraight(ON_Plane *frames, double s_tol)
Perform straightness test of surface.
int ON_Curve_Has_Tangent(const ON_Curve *curve, double ct_min, double ct_max, double t_tol)
Test whether a curve interval contains one or more horizontal or vertical tangents.
#define R
Definition: msr.c:55
bool ON_Surface_IsFlat(ON_Plane *frames, double f_tol)
Perform flatness test of surface.
#define A
Definition: msr.c:51
void distribute(const int count, const ON_3dVector *v, double x[], double y[], double z[])
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
bool ON_Surface_Quad_Split(const ON_Surface *surf, const ON_Interval &u, const ON_Interval &v, double upt, double vpt, ON_Surface **q0, ON_Surface **q1, ON_Surface **q2, ON_Surface **q3)
Create four sub-surfaces from a parent surface.
bool sortY(BRNode *first, BRNode *second)
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
HIDDEN const point_t delta
Definition: sh_prj.c:618
double fastf_t
Definition: defines.h:300
Definition: color.c:50
struct brlcad::_gcp_data GCPData