BRL-CAD
intersect.cpp
Go to the documentation of this file.
1 /* I N T E R S E C T . C P P
2  * BRL-CAD
3  *
4  * Copyright (c) 2013-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 intersect.cpp
21  *
22  * Implementation of intersection routines openNURBS left out.
23  *
24  */
25 
26 #include "common.h"
27 
28 #include <assert.h>
29 #include <vector>
30 #include <algorithm>
31 
32 #include "vmath.h"
33 
34 #include "bio.h"
35 #include "bu/log.h"
36 #include "bu/malloc.h"
37 #include "brep.h"
38 #include "brep_except.h"
39 
40 // Whether to output the debug messages about b-rep intersections.
41 #define DEBUG_BREP_INTERSECT 0
42 
43 // The maximal depth for subdivision - trade-off between accuracy and
44 // performance.
45 #define NR_MAX_DEPTH 8
46 #define MAX_PCI_DEPTH NR_MAX_DEPTH
47 #define MAX_PSI_DEPTH NR_MAX_DEPTH
48 #define MAX_CCI_DEPTH NR_MAX_DEPTH
49 #define MAX_CSI_DEPTH NR_MAX_DEPTH
50 #define MAX_SSI_DEPTH NR_MAX_DEPTH
51 
52 
53 // We make the default tolerance for PSI the same as that of curve and
54 // surface intersections defined by openNURBS (see opennurbs_curve.h
55 // and opennurbs_surface.h).
56 #define NR_DEFAULT_TOLERANCE 0.001
57 #define PCI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
58 #define PSI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
59 #define CCI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
60 #define CSI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
61 #define SSI_DEFAULT_TOLERANCE NR_DEFAULT_TOLERANCE
62 
63 // Used to prevent an infinite loop in the unlikely event that we
64 // can't provide a good starting point for the Newton-Raphson
65 // Iteration.
66 #define NR_MAX_ITERATIONS 100
67 #define PCI_MAX_ITERATIONS NR_MAX_ITERATIONS
68 #define PSI_MAX_ITERATIONS NR_MAX_ITERATIONS
69 #define CCI_MAX_ITERATIONS NR_MAX_ITERATIONS
70 #define CSI_MAX_ITERATIONS NR_MAX_ITERATIONS
71 #define SSI_MAX_ITERATIONS NR_MAX_ITERATIONS
72 
73 class XEventProxy {
74 public:
75  XEventProxy(ON_X_EVENT::TYPE type);
76 
77  void SetAPoint(const ON_3dPoint &pt);
78  void SetAPoints(const ON_3dPoint &start, const ON_3dPoint &end);
79 
80  void SetBPoint(const ON_3dPoint &pt);
81  void SetBPoints(const ON_3dPoint &start, const ON_3dPoint &end);
82 
83  void SetAOverlapRange(const ON_Interval &interval);
84  void SetAOverlapRange(double start, double end);
85  void SetAOverlapRange(double t);
86  void SetACurveParameter(double t);
87 
88  void SetBOverlapRange(const ON_Interval &interval);
89  void SetBOverlapRange(double start, double end);
90  void SetBOverlapRange(double t);
91  void SetBCurveParameter(double t);
92  void SetBSurfaceParameter(double u, double v);
93  void SetBSurfaceParameters(double u1, double v1, double u2, double v2);
94 
95  ON_X_EVENT Event(void);
96 
97 private:
98  ON_X_EVENT event;
99 };
100 
101 
102 XEventProxy::XEventProxy(ON_X_EVENT::TYPE type)
103 {
104  event.m_type = type;
105 }
106 
107 
108 void
109 XEventProxy::SetAPoint(const ON_3dPoint &pt)
110 {
111  event.m_A[0] = event.m_A[1] = pt;
112 }
113 
114 
115 void
116 XEventProxy::SetAPoints(const ON_3dPoint &start, const ON_3dPoint &end)
117 {
118  event.m_A[0] = start;
119  event.m_A[1] = end;
120 }
121 
122 
123 void
124 XEventProxy::SetBPoint(const ON_3dPoint &pt)
125 {
126  event.m_B[0] = event.m_B[1] = pt;
127 }
128 
129 
130 void
131 XEventProxy::SetBPoints(const ON_3dPoint &start, const ON_3dPoint &end)
132 {
133  event.m_B[0] = start;
134  event.m_B[1] = end;
135 }
136 
137 
138 void
139 XEventProxy::SetAOverlapRange(const ON_Interval &interval)
140 {
141  event.m_a[0] = interval.Min();
142  event.m_a[1] = interval.Max();
143 }
144 
145 
146 void
147 XEventProxy::SetAOverlapRange(double start, double end)
148 {
149  SetAOverlapRange(ON_Interval(start, end));
150 }
151 
152 
153 void
155 {
156  SetAOverlapRange(t, t);
157 }
158 
159 
160 void
162 {
163  event.m_a[0] = event.m_a[1] = t;
164 }
165 
166 
167 void
168 XEventProxy::SetBOverlapRange(const ON_Interval &interval)
169 {
170  event.m_b[0] = event.m_b[2] = interval.Min();
171  event.m_b[1] = event.m_b[3] = interval.Max();
172 }
173 
174 
175 void
176 XEventProxy::SetBOverlapRange(double start, double end)
177 {
178  SetBOverlapRange(ON_Interval(start, end));
179 }
180 
181 
182 void
184 {
185  SetBOverlapRange(t, t);
186 }
187 
188 
189 void
191 {
192  event.m_b[0] = event.m_b[1] = t;
193 }
194 
195 
196 void
197 XEventProxy::SetBSurfaceParameters(double u1, double v1, double u2, double v2)
198 {
199  event.m_b[0] = u1;
200  event.m_b[1] = v1;
201  event.m_b[2] = u2;
202  event.m_b[3] = v2;
203 }
204 
205 
206 void
208 {
209  SetBSurfaceParameters(u, v, u, v);
210 }
211 
212 
213 ON_X_EVENT
215 {
216  return event;
217 }
218 
219 ON_Curve *
220 sub_curve(const ON_Curve *in, double a, double b)
221 {
222  ON_Interval dom = in->Domain();
223 
224  if (b < a) {
225  std::swap(a, b);
226  }
227 
228  ON_Curve *min_to_a = NULL, *a_to_max = NULL;
229  ON_Curve *a_to_b = NULL, *b_to_max = NULL;
230 
231  if (a < dom.m_t[0] || b > dom.m_t[1]) {
232  throw InvalidInterval("sub_curve() interval outside curve domain\n");
233  }
234 
235  // Note that ON_SQRT_EPSILON is the tolerance used inside
236  // ON_PolylineCurve::Split() to determine if the split parameter
237  // should snap to an end of the domain.
238  if (ON_NearZero(a - dom.m_t[0], ON_SQRT_EPSILON)) {
239  if (ON_NearZero(dom.m_t[1] - b, ON_SQRT_EPSILON)) {
240  // a == dom.m_t[0] && b == dom.m_t[1]
241  a_to_b = in->Duplicate();
242  } else {
243  // a == dom.m_t[0] && b < dom.m_t[1]
244  in->Split(b, a_to_b, b_to_max);
245  delete b_to_max;
246  }
247  if (!a_to_b) {
248  // dom.m_t[0] == a == b
249  throw InvalidInterval("sub_curve() degenerate interval\n");
250  }
251  return a_to_b;
252  }
253 
254  if (ON_NearZero(dom.m_t[1] - b, ON_SQRT_EPSILON)) {
255  // a > dom.m_t[0] && b == dom.m_t[1]
256  in->Split(a, min_to_a, a_to_b);
257  delete min_to_a;
258 
259  if (!a_to_b) {
260  // a == b == dom.m_t[1]
261  throw InvalidInterval("sub_curve() degenerate interval\n");
262  }
263  return a_to_b;
264  }
265 
266  // a > dom.m_t[0] && b < dom.m_t[1]
267  in->Split(a, min_to_a, a_to_max);
268  delete min_to_a;
269 
270  if (!a_to_max) {
271  // a == b == dom.m_t[1]
272  throw InvalidInterval("sub_curve() interval is degenerate\n");
273  }
274  a_to_max->Split(b, a_to_b, b_to_max);
275  delete a_to_max;
276  delete b_to_max;
277 
278  if (!a_to_b) {
279  // a == b
280  throw InvalidInterval("sub_curve() interval is degenerate\n");
281  }
282  return a_to_b;
283 }
284 
285 ON_Surface *
286 sub_surface(const ON_Surface *in, int dir, double a, double b)
287 {
288  // approach: call ON_Surface::Split() twice with a and b respectively.
289  // [min, max] -> [min, a] & [a, max]
290  // [a, max] -> [a, b] & [b, max]
291 
292  ON_Interval dom = in->Domain(dir);
293  ON_Interval sub(a, b);
294  sub.MakeIncreasing();
295  if (!sub.Intersection(dom)) {
296  return NULL;
297  }
298  ON_Surface *left = NULL, *right = NULL, *three = NULL;
299 
300  in->Split(dir, sub.m_t[0], left, right);
301  if (left) {
302  delete left;
303  }
304  left = NULL;
305  if (!right) {
306  right = in->Duplicate();
307  }
308 
309  right->Split(dir, sub.m_t[1], left, three);
310  if (!left) {
311  left = right->Duplicate();
312  }
313 
314  if (right) {
315  delete right;
316  }
317  if (three) {
318  delete three;
319  }
320  return left;
321 }
322 
323 
324 HIDDEN const ON_Interval
325 check_domain(const ON_Interval *in, const ON_Interval &domain, const char *name)
326 {
327  // Check the validity of the input domain.
328  // Returns the intersection of the input interval and the curve/surface's
329  // domain.
330  // If the input interval is NULL, just returns the whole domain.
331 
332  if (in) {
333  if (!in->IsIncreasing()) {
334  bu_log("Bogus %s passed in", name);
335  return domain;
336  }
337  ON_Interval dom;
338  dom.Intersection(*in, domain);
339  if (dom.IsEmptyInterval()) {
340  bu_log("%s is disjoint from the whole domain.\n", name);
341  } else {
342  return dom;
343  }
344  }
345  return domain;
346 }
347 
348 
349 HIDDEN bool
350 build_curve_root(const ON_Curve *curve, const ON_Interval *domain, Subcurve &root)
351 {
352  // Build the curve tree root within a given domain
353 
354  if (curve == NULL) {
355  return false;
356  }
357 
358  if (domain == NULL || *domain == curve->Domain()) {
359  root.m_curve = curve->Duplicate();
360  root.m_t = curve->Domain();
361  } else {
362  // Call sub_curve() to get the curve segment inside the input domain.
363  try {
364  root.m_curve = sub_curve(curve, domain->Min(), domain->Max());
365  } catch (InvalidInterval &) {
366  root.m_curve = NULL;
367  }
368  root.m_t = *domain;
369  }
370 
371  if (root.m_curve) {
372  root.SetBBox(root.m_curve->BoundingBox());
373  root.m_islinear = root.m_curve->IsLinear();
374  return true;
375  } else {
376  return false;
377  }
378 }
379 
380 
381 HIDDEN Subcurve *
383  const ON_Curve *curve,
384  const ON_Interval *domain)
385 {
386  Subcurve *root = new Subcurve;
387  if (!build_curve_root(curve, domain, *root)) {
388  delete root;
389  root = NULL;
390  }
391  if (root == NULL) {
392  throw std::bad_alloc();
393  }
394  return root;
395 }
396 
397 
398 HIDDEN bool
399 build_surface_root(const ON_Surface *surf, const ON_Interval *u_domain, const ON_Interval *v_domain, Subsurface &root)
400 {
401  // Build the surface tree root within a given domain
402 
403  if (surf == NULL) {
404  return false;
405  }
406 
407  const ON_Surface *u_splitted_surf; // the surface after u-split
408 
409  // First, do u split
410  if (u_domain == NULL || *u_domain == surf->Domain(0)) {
411  u_splitted_surf = surf;
412  root.m_u = surf->Domain(0);
413  } else {
414  // Call sub_surface() to get the sub-surface inside the input U domain.
415  u_splitted_surf = sub_surface(surf, 0, u_domain->Min(), u_domain->Max());
416  root.m_u = *u_domain;
417  }
418 
419  if (u_splitted_surf == NULL) {
420  return false;
421  }
422 
423  // Then v split
424  if (v_domain == NULL || *v_domain == surf->Domain(1)) {
425  root.m_surf = u_splitted_surf->Duplicate();
426  root.m_v = surf->Domain(1);
427  } else {
428  // // Call sub_surface() to get the sub-surface inside the input V domain.
429  root.m_surf = sub_surface(surf, 1, v_domain->Min(), v_domain->Max());
430  root.m_v = *v_domain;
431  }
432 
433  if (root.m_surf) {
434  ON_BoundingBox surf_bbox;
435  int ret = root.m_surf->GetBoundingBox(surf_bbox.m_min,
436  surf_bbox.m_max);
437  if (!ret) {
438  double corners_min[3], corners_max[3];
439  for (int i = 0; i < 2; ++i) {
440  for (int j = 0; j < 2; ++j) {
441  ON_3dPoint corner = root.m_surf->PointAt(root.m_u.m_t[i],
442  root.m_v.m_t[i]);
443  VMINMAX(corners_min, corners_max, corner);
444  }
445  }
446  surf_bbox.m_min = ON_3dPoint(corners_min);
447  surf_bbox.m_max = ON_3dPoint(corners_max);
448  }
449  root.SetBBox(surf_bbox);
450  root.m_isplanar = root.m_surf->IsPlanar();
451  return true;
452  } else {
453  return false;
454  }
455 }
456 
457 
458 HIDDEN Subsurface *
460  const ON_Surface *surf,
461  const ON_Interval *u_domain,
462  const ON_Interval *v_domain)
463 {
464  Subsurface *root = new Subsurface;
465  if (!build_surface_root(surf, u_domain, v_domain, *root)) {
466  delete root;
467  root = NULL;
468  }
469  if (root == NULL) {
470  throw std::bad_alloc();
471  }
472  return root;
473 }
474 
475 
476 HIDDEN ON_Curve *
477 curve_fitting(ON_Curve *in, double fitting_tol = ON_ZERO_TOLERANCE, bool delete_curve = true)
478 {
479  // Fit a *2D* curve into line, arc, ellipse or other conic curves.
480 
481  if (in == NULL) {
482  return NULL;
483  }
484 
485  bool changed = false;
486 
487  // First, for a polyline, eliminate some unnecessary points (if three
488  // neighbor points are collinear, the middle one can be removed)
489 
490  ON_3dPointArray points;
491  if (in->IsPolyline(&points)) {
492  int point_count = points.Count();
493  ON_3dPointArray new_points;
494  ON_3dPoint start = points[0];
495  new_points.Append(start);
496  for (int i = 2; i < point_count; i++) {
497  ON_3dVector v1 = points[i - 1] - start;
498  ON_3dVector v2 = points[i] - points[i - 1];
499  if (!ON_NearZero(ON_CrossProduct(v1, v2).Length()) || ON_DotProduct(v1, v2) < -ON_ZERO_TOLERANCE) {
500  // start, points[i-1], points[i] are not collinear,
501  // or v1 and v2 have opposite directions.
502  // No points can be eliminated
503  start = points[i - 1];
504  new_points.Append(start);
505  }
506  }
507  new_points.Append(points[point_count - 1]);
508  if (new_points.Count() != point_count) {
509  // Some points have been eliminated
510  if (delete_curve) {
511  delete in;
512  }
513  in = new ON_PolylineCurve(new_points);
514  changed = true;
515  if (DEBUG_BREP_INTERSECT) {
516  bu_log("fitting: %d => %d points.\n", point_count, new_points.Count());
517  }
518  }
519  }
520 
521  // Linear fitting
522  if (in->IsLinear(fitting_tol)) {
523  ON_LineCurve *linecurve = new ON_LineCurve(in->PointAtStart(), in->PointAtEnd());
524  linecurve->ChangeDimension(in->Dimension());
525  if (delete_curve || changed) {
526  delete in;
527  }
528  return linecurve;
529  }
530 
531  // Conic fitting (ellipse, parabola, hyperbola)
532  // Now we only have the ellipse fitting
533 
534  // It's only meaningful to fit the curve when it's a complex one
535  // For a polyline curve, the number of points should not be less than 10.
536  const int fit_minimum_knots = 10;
537  int knot_count = in->SpanCount();
538  if (knot_count < fit_minimum_knots) {
539  return in;
540  }
541 
542  double *knots = new double [knot_count + 1];
543  in->GetSpanVector(knots);
544 
545  // Sample 6 points along the curve, and call ON_GetConicEquationThrough6Points().
546  double conic[6];
547  double sample_pts[12];
548  int plotres = in->IsClosed() ? 6 : 5;
549 
550  // The sample points should be knots (which are accurate if the curve
551  // is a polyline curve).
552  for (int i = 0; i < 6; i++) {
553  double knot_t = knots[ON_Round((double)i / plotres * knot_count)];
554  ON_3dPoint pt3d = in->PointAt(knot_t);
555  sample_pts[2 * i] = pt3d.x;
556  sample_pts[2 * i + 1] = pt3d.y;
557  }
558 
559  if (ON_GetConicEquationThrough6Points(2, sample_pts, conic, NULL, NULL, NULL)) {
560  // It may be a conic.
561  ON_2dPoint ell_center;
562  ON_2dVector ell_A, ell_B;
563  double ell_a, ell_b;
564  // First, fitting an ellipse. It seems that ON_Curve::IsEllipse()
565  // doesn't work, so we use conic fitting first. If this is not ideal,
566  // an alternative solution is to use least squares fitting on all
567  // points.
568  if (ON_IsConicEquationAnEllipse(conic, ell_center, ell_A, ell_B, &ell_a, &ell_b)) {
569  // The conic equation is an ellipse. But since we only considered
570  // the 6 sampled points, now we need to check whether all the points
571  // are on the ellipse.
572 
573  ON_Plane ell_plane = ON_Plane(ON_3dPoint(ell_center), ON_3dVector(ell_A), ON_3dVector(ell_B));
574  ON_Ellipse ell(ell_plane, ell_a, ell_b);
575  int i;
576  for (i = 0; i <= knot_count; i++) {
577  ON_3dPoint pt = in->PointAt(knots[i]);
578  ON_3dPoint ell_pt;
579  ell_pt = ell.ClosestPointTo(pt);
580  if (ell_pt.DistanceTo(pt) > fitting_tol) {
581  break;
582  }
583  }
584 
585  if (i == knot_count + 1) {
586  // All points are on the ellipse
587  ON_NurbsCurve nurbscurve;
588  ell.GetNurbForm(nurbscurve);
589 
590  // It may not be a complete ellipse.
591  // We find the domain of its parameter.
592  // The params of the nurbscurve is between [0, 2*pi]
593  double t_min, t_max;
594  if (in->IsClosed()) {
595  t_max = ON_PI * 2.0;
596  t_min = 0.0;
597  } else {
598  double t_mid;
599  ell.ClosestPointTo(in->PointAtStart(), &t_min);
600  ell.ClosestPointTo(in->PointAtEnd(), &t_max);
601  ell.ClosestPointTo(in->PointAt(knots[knot_count / 2]), &t_mid);
602  if (t_min > t_max) {
603  std::swap(t_min, t_max);
604  }
605  if (!ON_Interval(t_min, t_max).Includes(t_mid)) {
606  // The arc crosses where t = 0 (2*PI).
607  // We make the curve "two rounds" of the ellipse.
608  t_min += 2.0 * ON_PI;
609  std::swap(t_min, t_max);
610  nurbscurve.Append(nurbscurve);
611  }
612  }
613 
614  delete []knots;
615  if (delete_curve || changed) {
616  delete in;
617  }
618 
619  // The domain of the nurbscurve is [0, 4*pi]
620  // t_min \in [0, 2*pi], t_max \in [0, 4*pi]
621  return sub_curve(&nurbscurve, t_min, t_max);
622  }
623  }
624  }
625 
626  // We have tried all fittings, but none of them succeed.
627  // So we just return the original curve.
628  delete []knots;
629  return in;
630 }
631 
632 
633 /**
634  * Point-point intersections (PPI)
635  *
636  * approach:
637  *
638  * - Calculate the distance of the two points.
639  *
640  * - If the distance is within the tolerance, the two points
641  * intersect. The mid-point of them is the intersection point,
642  * and half of the distance is the uncertainty radius.
643  */
644 
645 
646 bool
647 ON_Intersect(const ON_3dPoint &pointA,
648  const ON_3dPoint &pointB,
649  ON_ClassArray<ON_PX_EVENT> &x,
650  double tol)
651 {
652  if (tol <= 0.0) {
653  tol = ON_ZERO_TOLERANCE;
654  }
655 
656  if (pointA.DistanceTo(pointB) <= tol) {
657  ON_PX_EVENT event;
658  event.m_type = ON_PX_EVENT::ppx_point;
659  event.m_A = pointA;
660  event.m_B = pointB;
661  event.m_Mid = (pointA + pointB) * 0.5;
662  event.m_radius = pointA.DistanceTo(pointB) * 0.5;
663  x.Append(event);
664  return true;
665  }
666  return false;
667 }
668 
669 
670 /**
671  * Point-curve intersections (PCI)
672  *
673  * approach:
674  *
675  * - Sub-divide the curve, calculating bounding boxes.
676  *
677  * - If the bounding box intersects with the point, go deeper until
678  * we reach the maximal depth or the sub-curve is linear.
679  *
680  * - Use linear approximation to get an estimated intersection point.
681  *
682  * - Use Newton-Raphson iterations to calculate an accurate intersection
683  * point, with the estimated one as a start.
684  */
685 HIDDEN bool
686 newton_pci(double &t, const ON_3dPoint &pointA, const ON_Curve &curveB, double tol)
687 {
688  ON_3dPoint closest_point = curveB.PointAt(t);
689  double dist = closest_point.DistanceTo(pointA);
690 
691  // Use Newton-Raphson method with the starting point from linear
692  // approximation.
693  // It's similar to point projection (or point inversion, or get
694  // closest point on curve). We calculate the root of:
695  // (curve(t) - pointA) \dot tangent(t) = 0 (*)
696  // Then curve(t) may be the closest point on the curve to pointA.
697  //
698  // In some cases, the closest point (intersection point) is the
699  // endpoints of the curve, it does not satisfy (*), but that's
700  // OK because it already comes out with the linear approximation
701  // above.
702 
703  double last_t = DBL_MAX;
704  int iterations = 0;
705  while (!ON_NearZero(last_t - t)
706  && iterations++ < PCI_MAX_ITERATIONS) {
707  ON_3dVector tangent, s_deriv;
708  last_t = t;
709  curveB.Ev2Der(t, closest_point, tangent, s_deriv);
710  // F = (curve(t) - pointA) \dot curve'(t)
711  // F' = (curve(t) - pointA) \dot curve''(t) + curve'(t) \dot curve'(t)
712  double F = ON_DotProduct(closest_point - pointA, tangent);
713  double derivF = ON_DotProduct(closest_point - pointA, s_deriv) + ON_DotProduct(tangent, tangent);
714  if (!ON_NearZero(derivF)) {
715  t -= F / derivF;
716  }
717  }
718 
719  t = std::min(curveB.Domain().Max(), t);
720  t = std::max(curveB.Domain().Min(), t);
721  closest_point = curveB.PointAt(t);
722  dist = closest_point.DistanceTo(pointA);
723  return dist <= tol && !isnan(t);
724 }
725 
726 
727 bool
728 ON_Intersect(const ON_3dPoint &pointA,
729  const ON_Curve &curveB,
730  ON_ClassArray<ON_PX_EVENT> &x,
731  double tol,
732  const ON_Interval *curveB_domain,
733  Subcurve *treeB)
734 {
735  if (tol <= 0.0) {
736  tol = PCI_DEFAULT_TOLERANCE;
737  }
738 
739  check_domain(curveB_domain, curveB.Domain(), "curveB_domain");
740 
741  Subcurve *root = treeB;
742  try {
743  if (root == NULL) {
744  root = get_curve_root(&curveB, curveB_domain);
745  }
746  } catch (std::bad_alloc &e) {
747  bu_log("%s", e.what());
748  return false;
749  }
750 
751  if (!root->IsPointIn(pointA, tol)) {
752  if (treeB == NULL) {
753  delete root;
754  }
755  return false;
756  }
757 
758  std::vector<Subcurve *> candidates, next_candidates;
759  candidates.push_back(root);
760 
761  // find the sub-curves that possibly intersect the point
762  for (int i = 0; i < MAX_PCI_DEPTH; i++) {
763  next_candidates.clear();
764  for (size_t j = 0; j < candidates.size(); j++) {
765  if (candidates[j]->m_islinear) {
766  next_candidates.push_back(candidates[j]);
767  } else {
768  if (candidates[j]->Split() == 0) {
769  if (candidates[j]->m_children[0]->IsPointIn(pointA, tol)) {
770  next_candidates.push_back(candidates[j]->m_children[0]);
771  }
772  if (candidates[j]->m_children[1]->IsPointIn(pointA, tol)) {
773  next_candidates.push_back(candidates[j]->m_children[1]);
774  }
775  } else {
776  next_candidates.push_back(candidates[j]);
777  }
778  }
779  }
780  candidates = next_candidates;
781  }
782 
783  for (size_t i = 0; i < candidates.size(); i++) {
784  // use linear approximation to get an estimated intersection point
785  ON_Line line(candidates[i]->m_curve->PointAtStart(), candidates[i]->m_curve->PointAtEnd());
786  double line_t;
787  line.ClosestPointTo(pointA, &line_t);
788 
789  // make sure line_t belongs to [0, 1]
790  CLAMP(line_t, 0.0, 1.0);
791  double closest_point_t = candidates[i]->m_t.ParameterAt(line_t);
792 
793  // use Newton iterations to get an accurate intersection point
794  if (newton_pci(closest_point_t, pointA, curveB, tol)) {
795  ON_PX_EVENT event;
796  event.m_type = ON_PX_EVENT::pcx_point;
797  event.m_A = pointA;
798  event.m_B = curveB.PointAt(closest_point_t);
799  event.m_b[0] = closest_point_t;
800  event.m_Mid = (pointA + event.m_B) * 0.5;
801  event.m_radius = event.m_B.DistanceTo(pointA) * 0.5;
802  x.Append(event);
803  if (treeB == NULL) {
804  delete root;
805  }
806  return true;
807  }
808  }
809 
810  // all candidates have no intersection
811  if (treeB == NULL) {
812  delete root;
813  }
814  return false;
815 }
816 
817 
818 /**
819  * Point-surface intersections (PSI)
820  *
821  * approach:
822  *
823  * - Sub-divide the surface, calculating bounding boxes.
824  *
825  * - If the bounding box intersects with the point, go deeper until
826  * we reach the maximal depth or the sub-surface is planar.
827  *
828  * - Use Newton-Raphson iterations to compute the closest point on
829  * the surface to that point.
830  *
831  * - If the closest point is within the given tolerance, there is an
832  * intersection.
833  */
834 HIDDEN bool
835 newton_psi(double &u, double &v, const ON_3dPoint &pointA, const ON_Surface &surfB, double tol)
836 {
837  ON_3dPoint closest_point = surfB.PointAt(u, v);
838  double dist = closest_point.DistanceTo(pointA);
839 
840  // Use Newton-Raphson method to get an accurate point of PSI.
841  // We actually calculate the closest point on the surface to that point.
842  // (surface(u, v) - pointA) \dot deriv_u(u, v) = 0 (1)
843  // (surface(u, v) - pointA) \dot deriv_v(u, v) = 0 (2)
844  // Then surface(u, v) may be the closest point on the surface to pointA.
845 
846  double last_u = DBL_MAX * .5, last_v = DBL_MAX * .5;
847  int iterations = 0;
848  while (fabs(last_u - u) + fabs(last_v - v) > ON_ZERO_TOLERANCE
849  && iterations++ < PSI_MAX_ITERATIONS) {
850  ON_3dVector du, dv, duu, duv, dvv;
851  last_u = u;
852  last_v = v;
853  surfB.Ev2Der(u, v, closest_point, du, dv, duu, duv, dvv);
854  ON_Matrix f(2, 1), j(2, 2);
855  // F0 = (surface(u, v) - pointA) \dot du(u, v);
856  f[0][0] = ON_DotProduct(closest_point - pointA, du);
857  // F1 = (surface(u, v) - pointA) \dot dv(u, v);
858  f[1][0] = ON_DotProduct(closest_point - pointA, dv);
859  // dF0/du = du(u, v) \dot du(u, v) + (surface(u, v) - pointA) \dot duu(u, v)
860  j[0][0] = ON_DotProduct(du, du) + ON_DotProduct(closest_point - pointA, duu);
861  // dF0/dv = dv(u, v) \dot du(u, v) + (surface(u, v) - pointA) \dot duv(u, v)
862  // dF1/du = du(u, v) \dot dv(u, v) + (surface(u, v) - pointA) \dot dvu(u, v)
863  // dF0/dv = dF1/du
864  j[0][1] = j[1][0] = ON_DotProduct(du, dv) + ON_DotProduct(closest_point - pointA, duv);
865  // dF1/dv = dv(u, v) \dot dv(u, v) + (surface(u, v) - pointA) \dot dvv(u, v)
866  j[1][1] = ON_DotProduct(dv, dv) + ON_DotProduct(closest_point - pointA, dvv);
867 
868  if (!j.Invert(0.0)) {
869  break;
870  }
871  ON_Matrix delta;
872  delta.Multiply(j, f);
873  u -= delta[0][0];
874  v -= delta[1][0];
875  }
876 
877  u = std::min(u, surfB.Domain(0).Max());
878  u = std::max(u, surfB.Domain(0).Min());
879  v = std::min(v, surfB.Domain(1).Max());
880  v = std::max(v, surfB.Domain(1).Min());
881  closest_point = surfB.PointAt(u, v);
882  dist = closest_point.DistanceTo(pointA);
883  return dist <= tol && !isnan(u) && !isnan(v);
884 }
885 
886 
887 bool
888 ON_Intersect(const ON_3dPoint &pointA,
889  const ON_Surface &surfaceB,
890  ON_ClassArray<ON_PX_EVENT> &x,
891  double tol,
892  const ON_Interval *surfaceB_udomain,
893  const ON_Interval *surfaceB_vdomain,
894  Subsurface *treeB)
895 {
896  if (tol <= 0.0) {
897  tol = PCI_DEFAULT_TOLERANCE;
898  }
899 
900  ON_Interval u_domain, v_domain;
901  u_domain = check_domain(surfaceB_udomain, surfaceB.Domain(0), "surfaceB_udomain");
902  v_domain = check_domain(surfaceB_vdomain, surfaceB.Domain(1), "surfaceB_vdomain");
903 
904  Subsurface *root = treeB;
905  try {
906  if (root == NULL) {
907  root = get_surface_root(&surfaceB, &u_domain, &v_domain);
908  }
909  } catch (std::bad_alloc &e) {
910  bu_log("%s", e.what());
911  return false;
912  }
913 
914  if (!root->IsPointIn(pointA, tol)) {
915  if (treeB == NULL) {
916  delete root;
917  }
918  return false;
919  }
920 
921  std::vector<Subsurface *> candidates, next_candidates;
922  candidates.push_back(root);
923 
924  // find the sub-surfaces that possibly intersect the point
925  for (int i = 0; i < MAX_PSI_DEPTH; i++) {
926  next_candidates.clear();
927  for (size_t j = 0; j < candidates.size(); j++) {
928  if (candidates[j]->m_isplanar) {
929  next_candidates.push_back(candidates[j]);
930  } else {
931  if (candidates[j]->Split() == 0) {
932  for (int k = 0; k < 4; k++) {
933  if (candidates[j]->m_children[k] && candidates[j]->m_children[k]->IsPointIn(pointA, tol)) {
934  next_candidates.push_back(candidates[j]->m_children[k]);
935  }
936  }
937  } else {
938  next_candidates.push_back(candidates[j]);
939  }
940  }
941  }
942  candidates = next_candidates;
943  }
944 
945  for (size_t i = 0; i < candidates.size(); i++) {
946  // Use the mid point of the bounding box as the starting point,
947  // and use Newton iterations to get an accurate intersection.
948  double u = candidates[i]->m_u.Mid(), v = candidates[i]->m_v.Mid();
949  if (newton_psi(u, v, pointA, surfaceB, tol)) {
950  ON_PX_EVENT event;
951  event.m_type = ON_PX_EVENT::psx_point;
952  event.m_A = pointA;
953  event.m_B = surfaceB.PointAt(u, v);
954  event.m_b[0] = u;
955  event.m_b[1] = v;
956  event.m_Mid = (pointA + event.m_B) * 0.5;
957  event.m_radius = event.m_B.DistanceTo(pointA) * 0.5;
958  x.Append(event);
959  if (treeB == NULL) {
960  delete root;
961  }
962  return true;
963  }
964  }
965 
966  // all candidates have no intersection
967  if (treeB == NULL) {
968  delete root;
969  }
970  return false;
971 }
972 
973 
974 /**
975  * Curve-curve intersection (CCI)
976  *
977  * approach:
978  *
979  * - Sub-divide the curves and calculate their bounding boxes.
980  *
981  * - Calculate the intersection of the bounding boxes, and go deeper if
982  * they intersect, until we reach the maximal depth or they are linear.
983  *
984  * - Use Newton-Raphson iterations to calculate an accurate intersection
985  * point.
986  *
987  * - Detect overlaps, eliminate duplications, and merge the continuous
988  * overlap events.
989  */
990 
991 // We can only test a finite number of points to determine overlap.
992 // Here we test 16 points uniformly distributed.
993 #define CCI_OVERLAP_TEST_POINTS 16
994 
995 HIDDEN void
996 newton_cci(double &t_a, double &t_b, const ON_Curve *curveA, const ON_Curve *curveB, double isect_tol)
997 {
998  // Equations:
999  // x_a(t_a) - x_b(t_b) = 0 (1)
1000  // y_a(t_a) - y_b(t_b) = 0 (2)
1001  // z_a(t_a) - z_b(t_b) = 0 (3)
1002  // It's an over-determined system.
1003  // We use Newton-Raphson iterations to solve two equations of the three,
1004  // and use the third for checking.
1005  double last_t_a = DBL_MAX * .5, last_t_b = DBL_MAX * .5;
1006  ON_3dPoint pointA = curveA->PointAt(t_a);
1007  ON_3dPoint pointB = curveB->PointAt(t_b);
1008  if (pointA.DistanceTo(pointB) < isect_tol) {
1009  return;
1010  }
1011 
1012  int iteration = 0;
1013  while (fabs(last_t_a - t_a) + fabs(last_t_b - t_b) > ON_ZERO_TOLERANCE
1014  && iteration++ < CCI_MAX_ITERATIONS) {
1015  last_t_a = t_a, last_t_b = t_b;
1016  ON_3dVector derivA, derivB;
1017  curveA->Ev1Der(t_a, pointA, derivA);
1018  curveB->Ev1Der(t_b, pointB, derivB);
1019  ON_Matrix j(2, 2), f(2, 1);
1020 
1021  // First try to solve (1), (2).
1022  j[0][0] = derivA.x;
1023  j[0][1] = -derivB.x;
1024  j[1][0] = derivA.y;
1025  j[1][1] = -derivB.y;
1026  f[0][0] = pointA.x - pointB.x;
1027  f[1][0] = pointA.y - pointB.y;
1028  if (!j.Invert(0.0)) {
1029  // The first try failed. Try to solve (1), (3).
1030  j[0][0] = derivA.x;
1031  j[0][1] = -derivB.x;
1032  j[1][0] = derivA.z;
1033  j[1][1] = -derivB.z;
1034  f[0][0] = pointA.x - pointB.x;
1035  f[1][0] = pointA.z - pointB.z;
1036  if (!j.Invert(0.0)) {
1037  // Failed again. Try to solve (2), (3).
1038  j[0][0] = derivA.y;
1039  j[0][1] = -derivB.y;
1040  j[1][0] = derivA.z;
1041  j[1][1] = -derivB.z;
1042  f[0][0] = pointA.y - pointB.y;
1043  f[1][0] = pointA.z - pointB.z;
1044  if (!j.Invert(0.0)) {
1045  // cannot find a root
1046  continue;
1047  }
1048  }
1049  }
1050  ON_Matrix delta;
1051  delta.Multiply(j, f);
1052  t_a -= delta[0][0];
1053  t_b -= delta[1][0];
1054  }
1055 
1056  // make sure the solution is inside the domains
1057  t_a = std::min(t_a, curveA->Domain().Max());
1058  t_a = std::max(t_a, curveA->Domain().Min());
1059  t_b = std::min(t_b, curveB->Domain().Max());
1060  t_b = std::max(t_b, curveB->Domain().Min());
1061 }
1062 
1063 
1064 static int
1065 compare_by_m_a0(const ON_X_EVENT *a, const ON_X_EVENT *b)
1066 {
1067  if (a->m_a[0] < b->m_a[0]) {
1068  return -1;
1069  }
1070  if (a->m_a[0] > b->m_a[0]) {
1071  return 1;
1072  }
1073  return 0;
1074 }
1075 
1076 
1077 double
1079  double tol_3d,
1080  Subcurve *root,
1081  const ON_Curve *curve,
1082  const ON_Interval *curve_domain)
1083 {
1084  double tol_2d = tol_3d;
1085  double bbox_diagonal_len = root->m_curve->BoundingBox().Diagonal().Length();
1086  if (!ON_NearZero(bbox_diagonal_len)) {
1087  const ON_Interval dom = curve_domain ? *curve_domain : curve->Domain();
1088  tol_2d = tol_3d / (bbox_diagonal_len * dom.Length());
1089  }
1090  return tol_2d;
1091 }
1092 
1093 
1094 int
1095 ON_Intersect(const ON_Curve *curveA,
1096  const ON_Curve *curveB,
1097  ON_SimpleArray<ON_X_EVENT> &x,
1098  double isect_tol,
1099  double overlap_tol,
1100  const ON_Interval *curveA_domain,
1101  const ON_Interval *curveB_domain,
1102  Subcurve *treeA,
1103  Subcurve *treeB)
1104 {
1105  if (curveA == NULL || curveB == NULL) {
1106  return 0;
1107  }
1108 
1109  int original_count = x.Count();
1110 
1111  if (isect_tol <= 0.0) {
1112  isect_tol = CCI_DEFAULT_TOLERANCE;
1113  }
1114  if (overlap_tol < isect_tol) {
1115  overlap_tol = 2.0 * isect_tol;
1116  }
1117  double t1_tol = isect_tol;
1118  double t2_tol = isect_tol;
1119 
1120  check_domain(curveA_domain, curveA->Domain(), "curveA_domain");
1121  check_domain(curveB_domain, curveB->Domain(), "curveB_domain");
1122 
1123  // handle degenerate cases (one or both of them is a "point")
1124  bool ispoint_curveA = curveA->BoundingBox().Diagonal().Length() < ON_ZERO_TOLERANCE;
1125  bool ispoint_curveB = curveB->BoundingBox().Diagonal().Length() < ON_ZERO_TOLERANCE;
1126  ON_3dPoint startpointA = curveA->PointAtStart();
1127  ON_3dPoint startpointB = curveB->PointAtStart();
1128  ON_ClassArray<ON_PX_EVENT> px_event;
1129  if (ispoint_curveA && ispoint_curveB) {
1130  // both curves are degenerate (point-point intersection)
1131  if (ON_Intersect(startpointA, startpointB, px_event, isect_tol)) {
1132  XEventProxy event(ON_X_EVENT::ccx_overlap);
1133  event.SetAPoints(startpointA, curveA->PointAtEnd());
1134  event.SetBPoints(startpointB, curveB->PointAtEnd());
1135  event.SetAOverlapRange(curveA->Domain());
1136  event.SetBOverlapRange(curveB->Domain());
1137  x.Append(event.Event());
1138  return 1;
1139  } else {
1140  return 0;
1141  }
1142  } else if (ispoint_curveA) {
1143  // curveA is degenerate (point-curve intersection)
1144  if (ON_Intersect(startpointA, *curveB, px_event, isect_tol)) {
1145  XEventProxy event(ON_X_EVENT::ccx_overlap);
1146  event.SetAPoints(startpointA, curveA->PointAtEnd());
1147  event.SetBPoint(px_event[0].m_B);
1148  event.SetAOverlapRange(curveA->Domain());
1149  event.SetBOverlapRange(px_event[0].m_b[0]);
1150  x.Append(event.Event());
1151  return 1;
1152  } else {
1153  return 0;
1154  }
1155  } else if (ispoint_curveB) {
1156  // curveB is degenerate (point-curve intersection)
1157  if (ON_Intersect(startpointB, *curveA, px_event, isect_tol)) {
1158  XEventProxy event(ON_X_EVENT::ccx_overlap);
1159  event.SetAPoint(px_event[0].m_B);
1160  event.SetBPoints(startpointB, curveB->PointAtEnd());
1161  event.SetAOverlapRange(px_event[0].m_b[0]);
1162  event.SetBOverlapRange(curveB->Domain());
1163  x.Append(event.Event());
1164  return 1;
1165  } else {
1166  return 0;
1167  }
1168  }
1169 
1170  // generate the roots of the two curve trees if necessary
1171  Subcurve *rootA = treeA;
1172  Subcurve *rootB = treeB;
1173  try {
1174  if (rootA == NULL) {
1175  rootA = get_curve_root(curveA, curveA_domain);
1176  }
1177  if (rootB == NULL) {
1178  rootB = get_curve_root(curveB, curveB_domain);
1179  }
1180  } catch (std::bad_alloc &e) {
1181  bu_log("%s", e.what());
1182  if (treeA == NULL) {
1183  delete rootA;
1184  }
1185  return 0;
1186  }
1187 
1188  t1_tol = tolerance_2d_from_3d(isect_tol, rootA, curveA,
1189  curveA_domain);
1190  t2_tol = tolerance_2d_from_3d(isect_tol, rootB, curveB,
1191  curveB_domain);
1192 
1193  typedef std::vector<std::pair<Subcurve *, Subcurve *> > NodePairs;
1194  NodePairs candidates, next_candidates;
1195  if (rootA->Intersect(*rootB, isect_tol)) {
1196  candidates.push_back(std::make_pair(rootA, rootB));
1197  }
1198 
1199  // use sub-division and bounding box intersections first
1200  for (int h = 0; h <= MAX_CCI_DEPTH; h++) {
1201  if (candidates.empty()) {
1202  break;
1203  }
1204  next_candidates.clear();
1205  for (NodePairs::iterator i = candidates.begin(); i != candidates.end(); i++) {
1206  std::vector<Subcurve *> splittedA, splittedB;
1207  if (i->first->m_islinear || i->first->Split() != 0) {
1208  splittedA.push_back(i->first);
1209  } else {
1210  splittedA.push_back(i->first->m_children[0]);
1211  splittedA.push_back(i->first->m_children[1]);
1212  }
1213  if (i->second->m_islinear || i->second->Split() != 0) {
1214  splittedB.push_back(i->second);
1215  } else {
1216  splittedB.push_back(i->second->m_children[0]);
1217  splittedB.push_back(i->second->m_children[1]);
1218  }
1219  // intersected children go to the next step
1220  for (size_t j = 0; j < splittedA.size(); j++)
1221  for (size_t k = 0; k < splittedB.size(); k++)
1222  if (splittedA[j]->Intersect(*splittedB[k], isect_tol)) {
1223  next_candidates.push_back(std::make_pair(splittedA[j], splittedB[k]));
1224  }
1225  }
1226  candidates = next_candidates;
1227  }
1228 
1229  ON_SimpleArray<ON_X_EVENT> tmp_x;
1230  // calculate an accurate intersection point for intersected
1231  // bounding boxes
1232  for (NodePairs::iterator i = candidates.begin(); i != candidates.end(); i++) {
1233  // special handling for linear curves
1234  if (i->first->m_islinear && i->second->m_islinear) {
1235  // if both are linear we use ON_IntersectLineLine
1236  ON_Line lineA(curveA->PointAt(i->first->m_t.Min()), curveA->PointAt(i->first->m_t.Max()));
1237  ON_Line lineB(curveB->PointAt(i->second->m_t.Min()), curveB->PointAt(i->second->m_t.Max()));
1238  if (lineA.Direction().IsParallelTo(lineB.Direction())) {
1239  double min_dist = lineA.MinimumDistanceTo(lineB);
1240 
1241  if (min_dist >= isect_tol) {
1242  // min_dist may not be accurate if endpoints are collinear
1243  double d1 = lineA.from.DistanceTo(lineB.from);
1244  double d2 = lineA.from.DistanceTo(lineB.to);
1245  double d3 = lineA.to.DistanceTo(lineB.from);
1246  double d4 = lineA.to.DistanceTo(lineB.to);
1247 
1248  min_dist = std::min(min_dist, std::min(d1, std::min(d2,
1249  std::min(d3, d4))));
1250  }
1251 
1252  if (min_dist < isect_tol) {
1253  // curves lie on the same line, may be single
1254  // point intersection or overlap
1255  double t_a1, t_a2, t_b1, t_b2;
1256  lineA.ClosestPointTo(lineB.from, &t_a1);
1257  lineA.ClosestPointTo(lineB.to, &t_a2);
1258  lineB.ClosestPointTo(lineA.from, &t_b1);
1259  lineB.ClosestPointTo(lineA.to, &t_b2);
1260 
1261  // Since ClosestPointTo treats the lines as
1262  // infinite, t values aren't necessarily in [0.0, 1.0].
1263  CLAMP(t_a1, 0.0, 1.0);
1264  CLAMP(t_a2, 0.0, 1.0);
1265  CLAMP(t_b1, 0.0, 1.0);
1266  CLAMP(t_b2, 0.0, 1.0);
1267 
1268  // t values may not be increasing if lines have
1269  // opposite directions.
1270  if (t_a1 > t_a2) {
1271  std::swap(t_a1, t_a2);
1272  }
1273  if (t_b1 > t_b2) {
1274  std::swap(t_b1, t_b2);
1275  }
1276 
1277  if (ON_NearZero(t_a2 - t_a1, t1_tol)) {
1278  // point intersection
1279  XEventProxy event(ON_X_EVENT::ccx_point);
1280  event.SetAPoint(lineA.PointAt(t_a1));
1281  event.SetBPoint(lineB.PointAt(t_b1));
1282  event.SetACurveParameter(i->first->m_t.ParameterAt(t_a1));
1283  event.SetBCurveParameter(i->second->m_t.ParameterAt(t_b1));
1284  tmp_x.Append(event.Event());
1285  } else {
1286  // overlap intersection
1287  XEventProxy event(ON_X_EVENT::ccx_overlap);
1288  event.SetAPoints(lineA.PointAt(t_a1), lineA.PointAt(t_a2));
1289  event.SetBPoints(lineB.PointAt(t_b1), lineB.PointAt(t_b2));
1290  event.SetAOverlapRange(i->first->m_t.ParameterAt(t_a1),
1291  i->first->m_t.ParameterAt(t_a2));
1292  event.SetBOverlapRange(i->second->m_t.ParameterAt(t_b1),
1293  i->second->m_t.ParameterAt(t_b2));
1294  tmp_x.Append(event.Event());
1295  }
1296  }
1297  } else {
1298  // not parallel, check intersection point
1299  double t_lineA, t_lineB;
1300  double t_a, t_b;
1301  if (ON_IntersectLineLine(lineA, lineB, &t_lineA, &t_lineB, isect_tol, true)) {
1302  t_a = i->first->m_t.ParameterAt(t_lineA);
1303  t_b = i->second->m_t.ParameterAt(t_lineB);
1304 
1305  XEventProxy event(ON_X_EVENT::ccx_point);
1306  event.SetAPoint(lineA.PointAt(t_lineA));
1307  event.SetBPoint(lineB.PointAt(t_lineB));
1308  event.SetACurveParameter(t_a);
1309  event.SetBCurveParameter(t_b);
1310  tmp_x.Append(event.Event());
1311  }
1312  }
1313  } else {
1314  // They are not both linear.
1315 
1316  // Use two different start points - the two end-points of the interval
1317  // If they converge to one point, it's considered an intersection
1318  // point, otherwise it's considered an overlap event.
1319  // FIXME: Find a better mechanism to check overlapping, because this method
1320  // may miss some overlap cases. (Overlap events can also converge to one
1321  // point).
1322  double t_a1 = i->first->m_t.Min(), t_b1 = i->second->m_t.Min();
1323  newton_cci(t_a1, t_b1, curveA, curveB, isect_tol);
1324  double t_a2 = i->first->m_t.Max(), t_b2 = i->second->m_t.Max();
1325  newton_cci(t_a2, t_b2, curveA, curveB, isect_tol);
1326  if (isnan(t_a1) || isnan(t_b1)) {
1327  // the first iteration result is not sufficient
1328  std::swap(t_a1, t_a2);
1329  std::swap(t_b1, t_b2);
1330  }
1331  if (isnan(t_a1) || isnan(t_b1)) {
1332  continue;
1333  }
1334 
1335  ON_3dPoint pointA1 = curveA->PointAt(t_a1);
1336  ON_3dPoint pointB1 = curveB->PointAt(t_b1);
1337  ON_3dPoint pointA2 = curveA->PointAt(t_a2);
1338  ON_3dPoint pointB2 = curveB->PointAt(t_b2);
1339  if ((pointA1.DistanceTo(pointA2) < isect_tol && pointB1.DistanceTo(pointB2) < isect_tol)
1340  || (isnan(t_a2) || isnan(t_b2))) {
1341  // it's considered the same point
1342  ON_3dPoint pointA = curveA->PointAt(t_a1);
1343  ON_3dPoint pointB = curveB->PointAt(t_b1);
1344  double distance = pointA.DistanceTo(pointB);
1345  // check the validity of the solution
1346  if (distance < isect_tol) {
1347  XEventProxy event(ON_X_EVENT::ccx_point);
1348  event.SetAPoint(pointA);
1349  event.SetBPoint(pointB);
1350  event.SetACurveParameter(t_a1);
1351  event.SetBCurveParameter(t_b1);
1352  tmp_x.Append(event.Event());
1353  }
1354  } else {
1355  // check overlap
1356  // bu_log("Maybe overlap.\n");
1357  double distance1 = pointA1.DistanceTo(pointB1);
1358  double distance2 = pointA2.DistanceTo(pointB2);
1359 
1360  // check the validity of the solution
1361  if (distance1 < isect_tol && distance2 < isect_tol) {
1362  int j;
1363  for (j = 1; j < CCI_OVERLAP_TEST_POINTS; j++) {
1364  double strike = 1.0 / CCI_OVERLAP_TEST_POINTS;
1365  ON_3dPoint test_point = curveA->PointAt(t_a1 + (t_a2 - t_a1) * j * strike);
1366  ON_ClassArray<ON_PX_EVENT> pci_x;
1367  if (!ON_Intersect(test_point, *curveB, pci_x, overlap_tol)) {
1368  break;
1369  }
1370  }
1371  if (j == CCI_OVERLAP_TEST_POINTS) {
1372  XEventProxy event(ON_X_EVENT::ccx_overlap);
1373  if (t_a1 <= t_a2) {
1374  event.SetAPoints(pointA1, pointA2);
1375  event.SetBPoints(pointB1, pointB2);
1376  event.SetAOverlapRange(t_a1, t_a2);
1377  event.SetBOverlapRange(t_b1, t_b2);
1378  } else {
1379  event.SetAPoints(pointA2, pointA1);
1380  event.SetBPoints(pointB2, pointB1);
1381  event.SetAOverlapRange(t_a2, t_a1);
1382  event.SetBOverlapRange(t_b2, t_b1);
1383  }
1384  tmp_x.Append(event.Event());
1385  continue;
1386  }
1387  // If j != CCI_OVERLAP_TEST_POINTS, two ccx_point events should
1388  // be generated. Fall through.
1389  }
1390  if (distance1 < isect_tol) {
1391  XEventProxy event(ON_X_EVENT::ccx_point);
1392  event.SetAPoint(pointA1);
1393  event.SetBPoint(pointB1);
1394  event.SetACurveParameter(t_a1);
1395  event.SetBCurveParameter(t_b1);
1396  tmp_x.Append(event.Event());
1397  }
1398  if (distance2 < isect_tol) {
1399  XEventProxy event(ON_X_EVENT::ccx_point);
1400  event.SetAPoint(pointA2);
1401  event.SetBPoint(pointB2);
1402  event.SetACurveParameter(t_a2);
1403  event.SetBCurveParameter(t_b2);
1404  tmp_x.Append(event.Event());
1405  }
1406  }
1407  }
1408  }
1409 
1410  // use an O(n^2) naive approach to eliminate duplications
1411  ON_SimpleArray<ON_X_EVENT> points, overlap;
1412  for (int i = 0; i < tmp_x.Count(); i++) {
1413  int j;
1414  if (tmp_x[i].m_type == ON_X_EVENT::ccx_overlap) {
1415  // handled later (merge)
1416  overlap.Append(tmp_x[i]);
1417  continue;
1418  }
1419  // ccx_point
1420  for (j = 0; j < points.Count(); j++) {
1421  if (tmp_x[i].m_A[0].DistanceTo(points[j].m_A[0]) < isect_tol
1422  && tmp_x[i].m_A[1].DistanceTo(points[j].m_A[1]) < isect_tol
1423  && tmp_x[i].m_B[0].DistanceTo(points[j].m_B[0]) < isect_tol
1424  && tmp_x[i].m_B[1].DistanceTo(points[j].m_B[1]) < isect_tol) {
1425  break;
1426  }
1427  }
1428  if (j == points.Count()) {
1429  points.Append(tmp_x[i]);
1430  }
1431  }
1432 
1433  // merge the overlap events that are continuous
1434  ON_X_EVENT pending;
1435  overlap.QuickSort(compare_by_m_a0);
1436  if (overlap.Count()) {
1437  pending = overlap[0];
1438  for (int i = 1; i < overlap.Count(); i++) {
1439  if (pending.m_a[1] < overlap[i].m_a[0] - t1_tol) {
1440  // pending[j] and overlap[i] are disjoint. pending[j] was
1441  // already merged, and should be removed from the list.
1442  // Because the elements in overlap[] are sorted by m_a[0],
1443  // the remaining elements won't intersect pending[j].
1444  pending.m_A[0] = curveA->PointAt(pending.m_a[0]);
1445  pending.m_A[1] = curveA->PointAt(pending.m_a[1]);
1446  pending.m_B[0] = curveB->PointAt(pending.m_b[0]);
1447  pending.m_B[1] = curveB->PointAt(pending.m_b[1]);
1448  x.Append(pending);
1449  pending = overlap[i];
1450  } else if (overlap[i].m_a[0] < pending.m_a[1] + t1_tol) {
1451  ON_Interval interval_1(overlap[i].m_b[0], overlap[i].m_b[1]);
1452  ON_Interval interval_2(pending.m_b[0], pending.m_b[1]);
1453  interval_1.MakeIncreasing();
1454  interval_1.m_t[0] -= t2_tol;
1455  interval_1.m_t[1] += t2_tol;
1456  if (interval_1.Intersection(interval_2)) {
1457  // need to merge: pending[j] = union(pending[j], overlap[i])
1458  pending.m_a[1] = std::max(overlap[i].m_a[1], pending.m_a[1]);
1459  pending.m_b[0] = std::min(overlap[i].m_b[0], pending.m_b[0]);
1460  pending.m_b[1] = std::max(overlap[i].m_b[1], pending.m_b[1]);
1461  }
1462  }
1463  }
1464  pending.m_A[0] = curveA->PointAt(pending.m_a[0]);
1465  pending.m_A[1] = curveA->PointAt(pending.m_a[1]);
1466  pending.m_B[0] = curveB->PointAt(pending.m_b[0]);
1467  pending.m_B[1] = curveB->PointAt(pending.m_b[1]);
1468  x.Append(pending);
1469  }
1470 
1471  // the intersection points shouldn't be inside the overlapped parts.
1472  int overlap_events = x.Count();
1473  for (int i = 0; i < points.Count(); i++) {
1474  int j;
1475  for (j = original_count; j < overlap_events; j++) {
1476  if (points[i].m_a[0] > x[j].m_a[0] - t1_tol
1477  && points[i].m_a[0] < x[j].m_a[1] + t1_tol) {
1478  break;
1479  }
1480  }
1481  if (j == overlap_events) {
1482  x.Append(points[i]);
1483  }
1484  }
1485 
1486  if (treeA == NULL) {
1487  delete rootA;
1488  }
1489  if (treeB == NULL) {
1490  delete rootB;
1491  }
1492  return x.Count() - original_count;
1493 }
1494 
1495 
1496 /**
1497  * Curve-surface intersection (CSI)
1498  *
1499  * approach:
1500  *
1501  * - Sub-divide the curve and the surface, calculate bounding boxes.
1502  *
1503  * - If their bounding boxes intersect, go deeper until we reach the maximal
1504  * depth, or the sub-curve is linear and the sub-surface is planar.
1505  *
1506  * - Use two starting points within the bounding boxes, and perform Newton-
1507  * Raphson iterations to get an accurate result.
1508  *
1509  * - Detect overlaps, eliminate duplications, and merge the continuous
1510  * overlap events.
1511  */
1512 
1513 // We can only test a finite number of points to determine overlap.
1514 // Here we test 2 points uniformly distributed.
1515 #define CSI_OVERLAP_TEST_POINTS 2
1516 
1517 HIDDEN void
1518 newton_csi(double &t, double &u, double &v, const ON_Curve *curveA, const ON_Surface *surfB, double isect_tol, Subsurface *tree)
1519 {
1520  // Equations:
1521  // C_x(t) - S_x(u, v) = 0
1522  // C_y(t) - S_y(u, v) = 0
1523  // C_z(t) - S_z(u, v) = 0
1524  // We use Newton-Raphson iterations to solve these equations.
1525  double last_t = DBL_MAX / 3, last_u = DBL_MAX / 3, last_v = DBL_MAX / 3;
1526  ON_3dPoint pointA = curveA->PointAt(t);
1527  ON_3dPoint pointB = surfB->PointAt(u, v);
1528  if (pointA.DistanceTo(pointB) < isect_tol) {
1529  return;
1530  }
1531 
1532  // If curveA(t) is already on surfB... Newton iterations may not succeed
1533  // with this case. Use point-surface intersections.
1534  ON_ClassArray<ON_PX_EVENT> px_event;
1535  if (ON_Intersect(pointA, *surfB, px_event, isect_tol, 0, 0, tree)) {
1536  u = px_event[0].m_b.x;
1537  v = px_event[0].m_b.y;
1538  return;
1539  }
1540 
1541  int iteration = 0;
1542  while (fabs(last_t - t) + fabs(last_u - u) + fabs(last_v - v) > ON_ZERO_TOLERANCE
1543  && iteration++ < CSI_MAX_ITERATIONS) {
1544  last_t = t, last_u = u, last_v = v;
1545  // F[0] = C_x(t) - S_x(u, v)
1546  // F[1] = C_y(t) - S_y(u, v)
1547  // F[2] = C_z(t) - S_z(u, v)
1548  // J[i][0] = dF[i]/dt = C'(t)
1549  // J[i][1] = dF[i]/du = -dS(u, v)/du
1550  // J[i][2] = dF[i]/dv = -dS(u, v)/dv
1551  ON_3dVector derivA, derivBu, derivBv;
1552  curveA->Ev1Der(t, pointA, derivA);
1553  surfB->Ev1Der(u, v, pointB, derivBu, derivBv);
1554  ON_Matrix j(3, 3), f(3, 1);
1555  for (int i = 0; i < 3; i++) {
1556  j[i][0] = derivA[i];
1557  j[i][1] = -derivBu[i];
1558  j[i][2] = -derivBv[i];
1559  f[i][0] = pointA[i] - pointB[i];
1560  }
1561  if (!j.Invert(0.0)) {
1562  break;
1563  }
1564  ON_Matrix delta;
1565  delta.Multiply(j, f);
1566  t -= delta[0][0];
1567  u -= delta[1][0];
1568  v -= delta[2][0];
1569  }
1570 
1571  // make sure the solution is inside the domains
1572  t = std::min(t, curveA->Domain().Max());
1573  t = std::max(t, curveA->Domain().Min());
1574  u = std::min(u, surfB->Domain(0).Max());
1575  u = std::max(u, surfB->Domain(0).Min());
1576  v = std::min(v, surfB->Domain(1).Max());
1577  v = std::max(v, surfB->Domain(1).Min());
1578 }
1579 
1580 
1581 double
1583  double tol_3d,
1584  Subsurface *root,
1585  const ON_Surface *surf,
1586  const ON_Interval *surf_udomain,
1587  const ON_Interval *surf_vdomain)
1588 {
1589  double tol_2d = tol_3d;
1590 
1591  const ON_Interval surf_udom = surf_udomain ? *surf_udomain : surf->Domain(0);
1592  const ON_Interval surf_vdom = surf_vdomain ? *surf_vdomain : surf->Domain(1);
1593  double surf_ulen = surf_udom.Length();
1594  double surf_vlen = surf_vdom.Length();
1595  double uv_diagonal_len = ON_2dVector(surf_ulen, surf_vlen).Length();
1596 
1597  double bbox_diagonal_len = root->m_surf->BoundingBox().Diagonal().Length();
1598  if (!ON_NearZero(bbox_diagonal_len)) {
1599  tol_2d = tol_3d / (bbox_diagonal_len * uv_diagonal_len);
1600  }
1601  return tol_2d;
1602 }
1603 
1604 
1605 int
1606 ON_Intersect(const ON_Curve *curveA,
1607  const ON_Surface *surfaceB,
1608  ON_SimpleArray<ON_X_EVENT> &x,
1609  double isect_tol,
1610  double overlap_tol,
1611  const ON_Interval *curveA_domain,
1612  const ON_Interval *surfaceB_udomain,
1613  const ON_Interval *surfaceB_vdomain,
1614  ON_CurveArray *overlap2d,
1615  Subcurve *treeA,
1616  Subsurface *treeB)
1617 {
1618  if (curveA == NULL || surfaceB == NULL) {
1619  return 0;
1620  }
1621 
1622  int original_count = x.Count();
1623 
1624  if (isect_tol <= 0.0) {
1625  isect_tol = CSI_DEFAULT_TOLERANCE;
1626  }
1627  if (overlap_tol < isect_tol) {
1628  overlap_tol = 2.0 * isect_tol;
1629  }
1630 
1631  check_domain(curveA_domain, curveA->Domain(), "curveA_domain");
1632  check_domain(surfaceB_udomain, surfaceB->Domain(0), "surfaceB_udomain");
1633  check_domain(surfaceB_vdomain, surfaceB->Domain(1), "surfaceB_vdomain");
1634 
1635  // if the curve is degenerated to a point...
1636  if (curveA->BoundingBox().Diagonal().Length() < ON_ZERO_TOLERANCE) {
1637  // point-surface intersections
1638  ON_3dPoint pointA = curveA->PointAtStart();
1639  ON_ClassArray<ON_PX_EVENT> px_event;
1640  if (ON_Intersect(pointA, *surfaceB, px_event, isect_tol)) {
1641  XEventProxy event(ON_X_EVENT::csx_overlap);
1642  event.SetAPoints(pointA, curveA->PointAtEnd());
1643  event.SetBPoint(px_event[0].m_B);
1644  event.SetAOverlapRange(curveA->Domain());
1645  event.SetBSurfaceParameter(px_event[0].m_b[0], px_event[0].m_b[1]);
1646  x.Append(event.Event());
1647  if (overlap2d) {
1648  ON_X_EVENT xevent = event.Event();
1649  overlap2d->Append(new ON_LineCurve(xevent.m_B[0], xevent.m_B[1]));
1650  }
1651  return 1;
1652  } else {
1653  return 0;
1654  }
1655  }
1656 
1657  // generate the roots of the two curve trees if necessary
1658  Subcurve *rootA = treeA;
1659  Subsurface *rootB = treeB;
1660  try {
1661  if (rootA == NULL) {
1662  rootA = get_curve_root(curveA, curveA_domain);
1663  }
1664  if (rootB == NULL) {
1665  rootB = get_surface_root(surfaceB, surfaceB_udomain,
1666  surfaceB_vdomain);
1667  }
1668  } catch (std::bad_alloc &e) {
1669  bu_log("%s", e.what());
1670  if (treeA == NULL) {
1671  delete rootA;
1672  }
1673  return 0;
1674  }
1675 
1676  // adjust the tolerance from 3D scale to respective 2D scales
1677  double t_tol = isect_tol;
1678  double u_tol = isect_tol;
1679  double v_tol = isect_tol;
1680  double l = rootA->m_curve->BoundingBox().Diagonal().Length();
1681  double dl = curveA_domain ? curveA_domain->Length() : curveA->Domain().Length();
1682  if (!ON_NearZero(l)) {
1683  t_tol = isect_tol / l * dl;
1684  }
1685  l = rootB->m_surf->BoundingBox().Diagonal().Length();
1686  dl = surfaceB_udomain ? surfaceB_udomain->Length() : surfaceB->Domain(0).Length();
1687  if (!ON_NearZero(l)) {
1688  u_tol = isect_tol / l * dl;
1689  }
1690  dl = surfaceB_vdomain ? surfaceB_vdomain->Length() : surfaceB->Domain(1).Length();
1691  if (!ON_NearZero(l)) {
1692  v_tol = isect_tol / l * dl;
1693  }
1694 
1695  typedef std::vector<std::pair<Subcurve *, Subsurface *> > NodePairs;
1696  NodePairs candidates, next_candidates;
1697  if (rootB->Intersect(*rootA, isect_tol)) {
1698  candidates.push_back(std::make_pair(rootA, rootB));
1699  }
1700 
1701  // use sub-division and bounding box intersections first
1702  for (int h = 0; h <= MAX_CSI_DEPTH; h++) {
1703  if (candidates.empty()) {
1704  break;
1705  }
1706  next_candidates.clear();
1707  for (NodePairs::iterator i = candidates.begin(); i != candidates.end(); i++) {
1708  std::vector<Subcurve *> splittedA;
1709  std::vector<Subsurface *> splittedB;
1710  if (i->first->m_islinear && i->second->m_isplanar) {
1711  splittedA.push_back(i->first);
1712  splittedB.push_back(i->second);
1713  } else {
1714  if (i->first->Split() != 0) {
1715  splittedA.push_back(i->first);
1716  } else {
1717  splittedA.push_back(i->first->m_children[0]);
1718  splittedA.push_back(i->first->m_children[1]);
1719  }
1720  if (i->second->Split() != 0) {
1721  splittedB.push_back(i->second);
1722  } else {
1723  for (int j = 0; j < 4; j++) {
1724  splittedB.push_back(i->second->m_children[j]);
1725  }
1726  }
1727  }
1728  for (size_t j = 0; j < splittedA.size(); j++)
1729  for (size_t k = 0; k < splittedB.size(); k++)
1730  if (splittedB[k]->Intersect(*splittedA[j], isect_tol)) {
1731  next_candidates.push_back(std::make_pair(splittedA[j], splittedB[k]));
1732  }
1733  }
1734  candidates = next_candidates;
1735  }
1736 
1737  ON_SimpleArray<ON_X_EVENT> tmp_x;
1738  // calculate an accurate intersection point for intersected
1739  // bounding boxes
1740  for (NodePairs::iterator i = candidates.begin(); i != candidates.end(); i++) {
1741  if (i->first->m_islinear && i->second->m_isplanar) {
1742  // try optimized approach for line-plane intersections
1743  ON_Line line(curveA->PointAt(i->first->m_t.Min()), curveA->PointAt(i->first->m_t.Max()));
1744  ON_Plane plane;
1745  bool success = true;
1746  ON_3dPoint point1 = surfaceB->PointAt(i->second->m_u.Min(), i->second->m_v.Min());
1747  ON_3dPoint point2 = surfaceB->PointAt(i->second->m_u.Min(), i->second->m_v.Max());
1748  ON_3dPoint point3 = surfaceB->PointAt(i->second->m_u.Max(), i->second->m_v.Max());
1749  ON_3dPoint point4 = surfaceB->PointAt(i->second->m_u.Max(), i->second->m_v.Min());
1750  if (!plane.CreateFromPoints(point1, point2, point3))
1751  if (!plane.CreateFromPoints(point1, point2, point4))
1752  if (!plane.CreateFromPoints(point1, point3, point4))
1753  if (!plane.CreateFromPoints(point2, point3, point4)) {
1754  success = false;
1755  }
1756 
1757  if (success && !ON_NearZero(line.Length())) {
1758  if (line.Direction().IsPerpendicularTo(plane.Normal())) {
1759  // they are parallel (or overlap)
1760 
1761  if (plane.DistanceTo(line.from) < isect_tol) {
1762  // The line is on the surface's plane. The end-points of
1763  // the overlap must be the linecurve's end-points or
1764  // the intersection between the linecurve and the boundary
1765  // of the surface.
1766 
1767  ON_X_EVENT event;
1768 
1769  // first we check the endpoints of the line segment (PSI)
1770  ON_ClassArray<ON_PX_EVENT> px_event1, px_event2;
1771  int intersections = 0;
1772  if (ON_Intersect(line.from, *surfaceB, px_event1, isect_tol, 0, 0, treeB)) {
1773  event.m_A[intersections] = line.from;
1774  event.m_B[intersections] = px_event1[0].m_B;
1775  event.m_a[intersections] = i->first->m_t.Min();
1776  event.m_b[2 * intersections] = px_event1[0].m_b[0];
1777  event.m_b[2 * intersections + 1] = px_event1[0].m_b[1];
1778  intersections++;
1779  }
1780  if (ON_Intersect(line.to, *surfaceB, px_event2, isect_tol, 0, 0, treeB)) {
1781  event.m_A[intersections] = line.to;
1782  event.m_B[intersections] = px_event2[0].m_B;
1783  event.m_a[intersections] = i->first->m_t.Max();
1784  event.m_b[2 * intersections] = px_event2[0].m_b[0];
1785  event.m_b[2 * intersections + 1] = px_event2[0].m_b[1];
1786  intersections++;
1787  }
1788 
1789  // then we check the intersection of the line segment
1790  // and the boundaries of the surface (CCI)
1791  ON_Line boundary[4];
1792  ON_SimpleArray<double> line_t, boundary_t;
1793  ON_SimpleArray<int> boundary_index;
1794  boundary[0] = ON_Line(point1, point2);
1795  boundary[1] = ON_Line(point2, point3);
1796  boundary[2] = ON_Line(point4, point3);
1797  boundary[3] = ON_Line(point1, point4);
1798  for (int j = 0; j < 4; j++) {
1799  double t1, t2;
1800  if (ON_IntersectLineLine(line, boundary[j], &t1, &t2, ON_ZERO_TOLERANCE, true)) {
1801  int k;
1802  for (k = 0; k < line_t.Count(); k++) {
1803  // Check duplication
1804  if (ON_NearZero(t1 - line_t[k])) {
1805  break;
1806  }
1807  }
1808  if (k == line_t.Count()) {
1809  line_t.Append(t1);
1810  boundary_t.Append(t2);
1811  boundary_index.Append(j);
1812  }
1813  }
1814  }
1815  int count = line_t.Count();
1816 
1817  for (int j = 0; j < count && intersections < 2; j++) {
1818  ON_3dPoint line_pt = line.PointAt(line_t[j]);
1819 
1820  if (intersections > 0 &&
1821  event.m_A[0].DistanceTo(line_pt) < isect_tol)
1822  {
1823  // skip duplicate intersection point
1824  continue;
1825  }
1826 
1827  // convert from the boxes' domain to the whole surface's domain
1828  double surf_u = 0.0, surf_v = 0.0;
1829  switch (boundary_index[j]) {
1830  case 0:
1831  surf_u = i->second->m_u.Min();
1832  surf_v = i->second->m_v.ParameterAt(boundary_t[j]);
1833  break;
1834  case 1:
1835  surf_u = i->second->m_u.ParameterAt(boundary_t[j]);
1836  surf_v = i->second->m_v.Max();
1837  break;
1838  case 2:
1839  surf_u = i->second->m_u.Max();
1840  surf_v = i->second->m_v.ParameterAt(boundary_t[j]);
1841  break;
1842  case 3:
1843  surf_u = i->second->m_u.ParameterAt(boundary_t[j]);
1844  surf_v = i->second->m_v.Min();
1845  break;
1846  }
1847  event.m_A[intersections] = line_pt;
1848  event.m_B[intersections] = surfaceB->PointAt(surf_u, surf_v);
1849  event.m_a[intersections] = i->first->m_t.ParameterAt(line_t[j]);
1850  event.m_b[2 * intersections] = surf_u;
1851  event.m_b[2 * intersections + 1] = surf_v;
1852  ++intersections;
1853  }
1854 
1855  // generate an ON_X_EVENT
1856  if (intersections == 0) {
1857  continue;
1858  }
1859 
1860  if (intersections == 1) {
1861  event.m_type = ON_X_EVENT::csx_point;
1862  event.m_A[1] = event.m_A[0];
1863  event.m_B[1] = event.m_B[0];
1864  event.m_a[1] = event.m_a[0];
1865  event.m_b[2] = event.m_b[0];
1866  event.m_b[3] = event.m_b[1];
1867  } else {
1868  event.m_type = ON_X_EVENT::csx_overlap;
1869  if (event.m_a[0] > event.m_a[1]) {
1870  std::swap(event.m_A[0], event.m_A[1]);
1871  std::swap(event.m_B[0], event.m_B[1]);
1872  std::swap(event.m_a[0], event.m_a[1]);
1873  std::swap(event.m_b[0], event.m_b[2]);
1874  std::swap(event.m_b[1], event.m_b[3]);
1875  }
1876  }
1877  tmp_x.Append(event);
1878  continue;
1879  } else {
1880  continue;
1881  }
1882  } else {
1883  // not parallel (not that this means cos_angle != 0), check intersection point
1884  double cos_angle = fabs(ON_DotProduct(plane.Normal(), line.Direction())) / (plane.Normal().Length() * line.Direction().Length());
1885  double distance = fabs(plane.DistanceTo(line.from));
1886  double line_t = distance / cos_angle / line.Length();
1887  if (line_t > 1.0 + t_tol) {
1888  continue;
1889  }
1890  ON_3dPoint intersection = line.from + line.Direction() * line_t;
1891  ON_ClassArray<ON_PX_EVENT> px_event;
1892  if (!ON_Intersect(intersection, *(i->second->m_surf), px_event, isect_tol, 0, 0, treeB)) {
1893  continue;
1894  }
1895 
1896  XEventProxy event(ON_X_EVENT::csx_point);
1897  event.SetAPoint(intersection);
1898  event.SetBPoint(px_event[0].m_B);
1899  event.SetACurveParameter(i->first->m_t.ParameterAt(line_t));
1900  event.SetBSurfaceParameter(px_event[0].m_b.x, px_event[0].m_b.y);
1901  tmp_x.Append(event.Event());
1902  }
1903  }
1904  }
1905 
1906  // They are not both linear or planar, or the above attempt failed.
1907 
1908  // Use two different start points - the two end-points of the interval
1909  // of the Subcurve's m_t.
1910  // If they converge to one point, it's considered an intersection
1911  // point, otherwise it's considered an overlap event.
1912  // FIXME: Find a better mechanism to check overlapping, because this method
1913  // may miss some overlap cases. (Overlap events can also converge to one
1914  // point).
1915 
1916  double u1 = i->second->m_u.Mid(), v1 = i->second->m_v.Mid();
1917  double t1 = i->first->m_t.Min();
1918  newton_csi(t1, u1, v1, curveA, surfaceB, isect_tol, treeB);
1919  double u2 = i->second->m_u.Mid(), v2 = i->second->m_v.Mid();
1920  double t2 = i->first->m_t.Max();
1921  newton_csi(t2, u2, v2, curveA, surfaceB, isect_tol, treeB);
1922 
1923  if (isnan(u1) || isnan(v1) || isnan(t1)) {
1924  u1 = u2;
1925  v1 = v2;
1926  t1 = t2;
1927  }
1928 
1929  if (isnan(u1) || isnan(v1) || isnan(t1)) {
1930  continue;
1931  }
1932 
1933  ON_3dPoint pointA1 = curveA->PointAt(t1);
1934  ON_3dPoint pointB1 = surfaceB->PointAt(u1, v1);
1935  ON_3dPoint pointA2 = curveA->PointAt(t2);
1936  ON_3dPoint pointB2 = surfaceB->PointAt(u2, v2);
1937  if (pointA1.DistanceTo(pointA2) < isect_tol
1938  && pointB1.DistanceTo(pointB2) < isect_tol) {
1939  // it's considered the same point (both converge)
1940  double distance = pointA1.DistanceTo(pointB1);
1941  // check the validity of the solution
1942  if (distance < isect_tol) {
1943  XEventProxy event(ON_X_EVENT::csx_point);
1944  event.SetAPoint(pointA1);
1945  event.SetBPoint(pointB1);
1946  event.SetACurveParameter(t1);
1947  event.SetBSurfaceParameter(u1, v1);
1948  tmp_x.Append(event.Event());
1949  }
1950  } else {
1951  // check overlap
1952  // bu_log("Maybe overlap.\n");
1953  double distance1 = pointA1.DistanceTo(pointB1);
1954  double distance2 = pointA2.DistanceTo(pointB2);
1955 
1956  // check the validity of the solution
1957  if (distance1 < isect_tol && distance2 < isect_tol) {
1958  XEventProxy event(ON_X_EVENT::csx_overlap);
1959  // make sure that m_a[0] <= m_a[1]
1960  if (t1 <= t2) {
1961  event.SetAPoints(pointA1, pointA2);
1962  event.SetBPoints(pointB1, pointB2);
1963  event.SetAOverlapRange(t1, t2);
1964  event.SetBSurfaceParameters(u1, v1, u2, v2);
1965  } else {
1966  event.SetAPoints(pointA2, pointA1);
1967  event.SetBPoints(pointB2, pointB1);
1968  event.SetAOverlapRange(t2, t1);
1969  event.SetBSurfaceParameters(u2, v2, u1, v1);
1970  }
1971  int j;
1972  for (j = 1; j < CSI_OVERLAP_TEST_POINTS; j++) {
1973  double strike = 1.0 / CSI_OVERLAP_TEST_POINTS;
1974  ON_3dPoint test_point = curveA->PointAt(t1 + (t2 - t1) * j * strike);
1975  ON_ClassArray<ON_PX_EVENT> psi_x;
1976  if (!ON_Intersect(test_point, *surfaceB, psi_x, overlap_tol, 0, 0, treeB)) {
1977  break;
1978  }
1979  }
1980  if (j == CSI_OVERLAP_TEST_POINTS) {
1981  tmp_x.Append(event.Event());
1982  continue;
1983  }
1984  // If j != CSI_OVERLAP_TEST_POINTS, two csx_point events may
1985  // be generated. Fall through.
1986  }
1987  if (distance1 < isect_tol) {
1988  XEventProxy event(ON_X_EVENT::csx_point);
1989  event.SetAPoint(pointA1);
1990  event.SetBPoint(pointB1);
1991  event.SetACurveParameter(t1);
1992  event.SetBSurfaceParameter(u1, v1);
1993  tmp_x.Append(event.Event());
1994  }
1995  if (distance2 < isect_tol) {
1996  XEventProxy event(ON_X_EVENT::csx_point);
1997  event.SetAPoint(pointA2);
1998  event.SetBPoint(pointB2);
1999  event.SetACurveParameter(t2);
2000  event.SetBSurfaceParameter(u2, v2);
2001  tmp_x.Append(event.Event());
2002  }
2003  }
2004  }
2005 
2006  ON_SimpleArray<ON_X_EVENT> points, overlap;
2007  // use an O(n^2) naive approach to eliminate duplications
2008  for (int i = 0; i < tmp_x.Count(); i++) {
2009  int j;
2010  if (tmp_x[i].m_type == ON_X_EVENT::csx_overlap) {
2011  overlap.Append(tmp_x[i]);
2012  continue;
2013  }
2014  // csx_point
2015  for (j = 0; j < points.Count(); j++) {
2016  if (tmp_x[i].m_A[0].DistanceTo(points[j].m_A[0]) < isect_tol
2017  && tmp_x[i].m_A[1].DistanceTo(points[j].m_A[1]) < isect_tol
2018  && tmp_x[i].m_B[0].DistanceTo(points[j].m_B[0]) < isect_tol
2019  && tmp_x[i].m_B[1].DistanceTo(points[j].m_B[1]) < isect_tol) {
2020  break;
2021  }
2022  }
2023  if (j == points.Count()) {
2024  points.Append(tmp_x[i]);
2025  }
2026  }
2027 
2028  // merge the overlap events that are continuous
2029  ON_X_EVENT pending;
2030  ON_3dPointArray ptarrayB;
2031  overlap.QuickSort(compare_by_m_a0);
2032  if (overlap.Count()) {
2033  pending = overlap[0];
2034  ptarrayB.Append(ON_3dPoint(pending.m_b[0], pending.m_b[1], 0.0));
2035  for (int i = 1; i < overlap.Count(); i++) {
2036  if (pending.m_a[1] < overlap[i].m_a[0] - t_tol) {
2037  // pending[j] and overlap[i] are disjoint. pending[j] was
2038  // already merged, and should be removed from the list.
2039  // Because the elements in overlap[] are sorted by m_a[0],
2040  // the remaining elements won't intersect pending[j].
2041  x.Append(pending);
2042  if (overlap2d) {
2043  ptarrayB.Append(ON_3dPoint(pending.m_b[2], pending.m_b[3], 0.0));
2044  ON_PolylineCurve *polyline = new ON_PolylineCurve(ptarrayB);
2045  polyline->ChangeDimension(2);
2046  overlap2d->Append(curve_fitting(polyline));
2047  ptarrayB.Empty();
2048  ptarrayB.Append(ON_3dPoint(overlap[i].m_b[0], overlap[i].m_b[1], 0.0));
2049  }
2050  pending = overlap[i];
2051  } else if (overlap[i].m_a[0] < pending.m_a[1] + t_tol) {
2052  if (overlap2d) {
2053  ptarrayB.Append(ON_3dPoint(overlap[i].m_b[0], overlap[i].m_b[1], 0.0));
2054  }
2055  ON_Interval interval_u1(overlap[i].m_b[0], overlap[i].m_b[2]);
2056  ON_Interval interval_u2(pending.m_b[0], pending.m_b[2]);
2057  ON_Interval interval_v1(overlap[i].m_b[1], overlap[i].m_b[3]);
2058  ON_Interval interval_v2(pending.m_b[1], pending.m_b[3]);
2059  interval_u1.MakeIncreasing();
2060  interval_u1.m_t[0] -= u_tol;
2061  interval_u1.m_t[1] += u_tol;
2062  interval_v1.MakeIncreasing();
2063  interval_v1.m_t[0] -= v_tol;
2064  interval_v1.m_t[1] += v_tol;
2065  if (interval_u1.Intersection(interval_u2) && interval_v1.Intersection(interval_v2)) {
2066  // If the uv rectangle of them intersects, it's considered overlap.
2067  // need to merge: pending[j] = union(pending[j], overlap[i])
2068  if (overlap[i].m_a[1] > pending.m_a[1]) {
2069  pending.m_a[1] = overlap[i].m_a[1];
2070  pending.m_b[2] = overlap[i].m_b[2];
2071  pending.m_b[3] = overlap[i].m_b[3];
2072  pending.m_A[1] = overlap[i].m_A[1];
2073  pending.m_B[1] = overlap[i].m_B[1];
2074  }
2075  }
2076  }
2077  }
2078  x.Append(pending);
2079  if (overlap2d) {
2080  ptarrayB.Append(ON_3dPoint(pending.m_b[2], pending.m_b[3], 0.0));
2081  ON_PolylineCurve *polyline = new ON_PolylineCurve(ptarrayB);
2082  polyline->ChangeDimension(2);
2083  overlap2d->Append(curve_fitting(polyline));
2084  }
2085  }
2086 
2087  // The intersection points shouldn't be inside the overlapped parts.
2088  int overlap_events = x.Count();
2089  for (int i = 0; i < points.Count(); i++) {
2090  int j;
2091  for (j = original_count; j < overlap_events; j++) {
2092  if (points[i].m_a[0] > x[j].m_a[0] - t_tol
2093  && points[i].m_a[0] < x[j].m_a[1] + t_tol) {
2094  break;
2095  }
2096  }
2097  if (j == overlap_events) {
2098  x.Append(points[i]);
2099  if (overlap2d) {
2100  overlap2d->Append(NULL);
2101  }
2102  }
2103  }
2104 
2105  if (treeA == NULL) {
2106  delete rootA;
2107  }
2108  if (treeB == NULL) {
2109  delete rootB;
2110  }
2111  return x.Count() - original_count;
2112 }
2113 
2114 
2115 /**
2116  * Surface-surface intersections (SSI)
2117  *
2118  * approach:
2119  *
2120  * - Generate the bounding box of the two surfaces.
2121  *
2122  * - If their bounding boxes intersect:
2123  *
2124  * -- Split the two surfaces, both into four parts, and calculate the
2125  * sub-surfaces' bounding boxes.
2126  *
2127  * -- Calculate the intersection of sub-surfaces' bboxes, if they do
2128  * intersect, go deeper with splitting surfaces and smaller bboxes,
2129  * otherwise trace back.
2130  *
2131  * - After getting the intersecting bboxes, approximate the surface
2132  * inside the bbox as two triangles, and calculate the intersection
2133  * points of the triangles (both in 3d space and two surfaces' UV
2134  * space).
2135  *
2136  * - If an intersection point is accurate enough (within the tolerance),
2137  * we append it to the list. Otherwise, we use Newton-Raphson iterations
2138  * to solve an under-determined system to get an more accurate inter-
2139  * section point.
2140  *
2141  * - Fit the intersection points into polyline curves, and then to NURBS
2142  * curves. Points with distance less than max_dist are considered in
2143  * one curve. Points that cannot be fitted to a curve are regarded as
2144  * single points.
2145  *
2146  * - Linear fittings and conic fittings are done to simplify the curves'
2147  * representations.
2148  *
2149  * - Note: the overlap cases are handled at the beginning separately.
2150  *
2151  * See: Adarsh Krishnamurthy, Rahul Khardekar, Sara McMains, Kirk Haller,
2152  * and Gershon Elber. 2008.
2153  * Performing efficient NURBS modeling operations on the GPU.
2154  * In Proceedings of the 2008 ACM symposium on Solid and physical modeling
2155  * (SPM '08). ACM, New York, NY, USA, 257-268. DOI=10.1145/1364901.1364937
2156  * http://doi.acm.org/10.1145/1364901.1364937
2157  */
2158 struct Triangle {
2159  ON_3dPoint a, b, c;
2160  ON_2dPoint a_2d, b_2d, c_2d;
2161  inline void CreateFromPoints(ON_3dPoint &pA, ON_3dPoint &pB, ON_3dPoint &pC)
2162  {
2163  a = pA;
2164  b = pB;
2165  c = pC;
2166  }
2167  ON_3dPoint BarycentricCoordinate(const ON_3dPoint &pt) const
2168  {
2169  double x, y, z, pivot_ratio;
2170  if (ON_Solve2x2(a[0] - c[0], b[0] - c[0], a[1] - c[1], b[1] - c[1], pt[0] - c[0], pt[1] - c[1], &x, &y, &pivot_ratio) == 2) {
2171  return ON_3dPoint(x, y, 1 - x - y);
2172  }
2173  if (ON_Solve2x2(a[0] - c[0], b[0] - c[0], a[2] - c[2], b[2] - c[2], pt[0] - c[0], pt[2] - c[2], &x, &z, &pivot_ratio) == 2) {
2174  return ON_3dPoint(x, 1 - x - z, z);
2175  }
2176  if (ON_Solve2x2(a[1] - c[1], b[1] - c[1], a[2] - c[2], b[2] - c[2], pt[1] - c[1], pt[2] - c[2], &y, &z, &pivot_ratio) == 2) {
2177  return ON_3dPoint(1 - y - z, y, z);
2178  }
2179  return ON_3dPoint::UnsetPoint;
2180  }
2181  void GetLineSegments(ON_Line line[3]) const
2182  {
2183  line[0] = ON_Line(a, b);
2184  line[1] = ON_Line(a, c);
2185  line[2] = ON_Line(b, c);
2186  }
2187 };
2188 
2189 
2190 HIDDEN bool
2191 triangle_intersection(const struct Triangle &triA, const struct Triangle &triB, ON_3dPoint &center)
2192 {
2193  ON_Plane plane_a(triA.a, triA.b, triA.c);
2194  ON_Plane plane_b(triB.a, triB.b, triB.c);
2195  ON_Line intersect;
2196  if (plane_a.Normal().IsParallelTo(plane_b.Normal())) {
2197  if (!ON_NearZero(plane_a.DistanceTo(plane_b.origin))) {
2198  // parallel
2199  return false;
2200  }
2201  // the two triangles are in the same plane
2202  ON_3dPoint pt(0.0, 0.0, 0.0);
2203  int count = 0;
2204  ON_Line lineA[3], lineB[3];
2205  triA.GetLineSegments(lineA);
2206  triB.GetLineSegments(lineB);
2207  for (int i = 0; i < 3; i++) {
2208  for (int j = 0; j < 3; j++) {
2209  double t_a, t_b;
2210  if (ON_IntersectLineLine(lineA[i], lineB[j], &t_a, &t_b, ON_ZERO_TOLERANCE, true)) {
2211  // we don't care if they are parallel or coincident
2212  pt += lineA[i].PointAt(t_a);
2213  count++;
2214  }
2215  }
2216  }
2217  if (!count) {
2218  return false;
2219  }
2220  center = pt / count;
2221  return true;
2222  }
2223 
2224  if (!ON_Intersect(plane_a, plane_b, intersect)) {
2225  return false;
2226  }
2227  ON_3dVector line_normal = ON_CrossProduct(plane_a.Normal(), intersect.Direction());
2228 
2229  // dpi - >0: one side of the line, <0: another side, ==0: on the line
2230  double dp1 = ON_DotProduct(triA.a - intersect.from, line_normal);
2231  double dp2 = ON_DotProduct(triA.b - intersect.from, line_normal);
2232  double dp3 = ON_DotProduct(triA.c - intersect.from, line_normal);
2233 
2234  int points_on_one_side = (dp1 >= ON_ZERO_TOLERANCE) + (dp2 >= ON_ZERO_TOLERANCE) + (dp3 >= ON_ZERO_TOLERANCE);
2235  if (points_on_one_side == 0 || points_on_one_side == 3) {
2236  return false;
2237  }
2238 
2239  line_normal = ON_CrossProduct(plane_b.Normal(), intersect.Direction());
2240  double dp4 = ON_DotProduct(triB.a - intersect.from, line_normal);
2241  double dp5 = ON_DotProduct(triB.b - intersect.from, line_normal);
2242  double dp6 = ON_DotProduct(triB.c - intersect.from, line_normal);
2243  points_on_one_side = (dp4 >= ON_ZERO_TOLERANCE) + (dp5 >= ON_ZERO_TOLERANCE) + (dp6 >= ON_ZERO_TOLERANCE);
2244  if (points_on_one_side == 0 || points_on_one_side == 3) {
2245  return false;
2246  }
2247 
2248  double t[4];
2249  int count = 0;
2250  double t1, t2;
2251  // dpi*dpj <= 0 - the corresponding points are on different sides,
2252  // or at least one of them is on the intersection line
2253  // - the line segment between them is intersected by the plane-plane
2254  // intersection line
2255  if (dp1 * dp2 < ON_ZERO_TOLERANCE) {
2256  intersect.ClosestPointTo(triA.a, &t1);
2257  intersect.ClosestPointTo(triA.b, &t2);
2258  double d1 = triA.a.DistanceTo(intersect.PointAt(t1));
2259  double d2 = triA.b.DistanceTo(intersect.PointAt(t2));
2260  if (!ZERO(d1 + d2)) {
2261  t[count++] = t1 + d1 / (d1 + d2) * (t2 - t1);
2262  }
2263  }
2264  if (dp1 * dp3 < ON_ZERO_TOLERANCE) {
2265  intersect.ClosestPointTo(triA.a, &t1);
2266  intersect.ClosestPointTo(triA.c, &t2);
2267  double d1 = triA.a.DistanceTo(intersect.PointAt(t1));
2268  double d2 = triA.c.DistanceTo(intersect.PointAt(t2));
2269  if (!ZERO(d1 + d2)) {
2270  t[count++] = t1 + d1 / (d1 + d2) * (t2 - t1);
2271  }
2272  }
2273  if (dp2 * dp3 < ON_ZERO_TOLERANCE && count != 2) {
2274  intersect.ClosestPointTo(triA.b, &t1);
2275  intersect.ClosestPointTo(triA.c, &t2);
2276  double d1 = triA.b.DistanceTo(intersect.PointAt(t1));
2277  double d2 = triA.c.DistanceTo(intersect.PointAt(t2));
2278  if (!ZERO(d1 + d2)) {
2279  t[count++] = t1 + d1 / (d1 + d2) * (t2 - t1);
2280  }
2281  }
2282  if (count != 2) {
2283  return false;
2284  }
2285  if (dp4 * dp5 < ON_ZERO_TOLERANCE) {
2286  intersect.ClosestPointTo(triB.a, &t1);
2287  intersect.ClosestPointTo(triB.b, &t2);
2288  double d1 = triB.a.DistanceTo(intersect.PointAt(t1));
2289  double d2 = triB.b.DistanceTo(intersect.PointAt(t2));
2290  if (!ZERO(d1 + d2)) {
2291  t[count++] = t1 + d1 / (d1 + d2) * (t2 - t1);
2292  }
2293  }
2294  if (dp4 * dp6 < ON_ZERO_TOLERANCE) {
2295  intersect.ClosestPointTo(triB.a, &t1);
2296  intersect.ClosestPointTo(triB.c, &t2);
2297  double d1 = triB.a.DistanceTo(intersect.PointAt(t1));
2298  double d2 = triB.c.DistanceTo(intersect.PointAt(t2));
2299  if (!ZERO(d1 + d2)) {
2300  t[count++] = t1 + d1 / (d1 + d2) * (t2 - t1);
2301  }
2302  }
2303  if (dp5 * dp6 < ON_ZERO_TOLERANCE && count != 4) {
2304  intersect.ClosestPointTo(triB.b, &t1);
2305  intersect.ClosestPointTo(triB.c, &t2);
2306  double d1 = triB.b.DistanceTo(intersect.PointAt(t1));
2307  double d2 = triB.c.DistanceTo(intersect.PointAt(t2));
2308  if (!ZERO(d1 + d2)) {
2309  t[count++] = t1 + d1 / (d1 + d2) * (t2 - t1);
2310  }
2311  }
2312  if (count != 4) {
2313  return false;
2314  }
2315  if (t[0] > t[1]) {
2316  std::swap(t[1], t[0]);
2317  }
2318  if (t[2] > t[3]) {
2319  std::swap(t[3], t[2]);
2320  }
2321  double left = std::max(t[0], t[2]);
2322  double right = std::min(t[1], t[3]);
2323  if (left > right) {
2324  return false;
2325  }
2326  center = intersect.PointAt((left + right) / 2);
2327  return true;
2328 }
2329 
2330 
2331 struct PointPair {
2334  double tol;
2335  bool operator < (const PointPair &pp) const
2336  {
2337  if (ON_NearZero(distance3d - pp.distance3d, tol)) {
2338  return dist_uA + dist_vA + dist_uB + dist_vB < pp.dist_uA + pp.dist_vA + pp.dist_uB + pp.dist_vB;
2339  }
2340  return distance3d < pp.distance3d;
2341  }
2342 };
2343 
2344 
2345 HIDDEN void
2347  const ON_Surface *surfA,
2348  const ON_Surface *surfB,
2349  ON_3dPointArray &curvept,
2350  ON_2dPointArray &curve_uvA,
2351  ON_2dPointArray &curve_uvB)
2352 {
2353 
2354  for (int i = 0; i < 4; ++i) {
2355  int dir = i % 2;
2356  bool doA = i < 2;
2357 
2358  ON_BOOL32 is_closed = (doA ? surfA : surfB)->IsClosed(dir);
2359  if (!is_closed) {
2360  continue;
2361  }
2362  const ON_Interval domain = (doA ? surfA : surfB)->Domain(dir);
2363  double dmin = domain.Min();
2364  double dmax = domain.Max();
2365 
2366  int orig_count = curvept.Count();
2367  for (int j = 0; j < orig_count; ++j) {
2368  ON_2dPoint uvA = curve_uvA[j];
2369  ON_2dPoint uvB = curve_uvB[j];
2370  double &dval = (doA ? uvA : uvB)[dir];
2371 
2372  // if this uv point is NEAR the closed seam, add a new
2373  // point ON the seam with the same 3D coordinates
2374  if (ON_NearZero(dval - dmin)) {
2375  dval = dmax;
2376  } else if (ON_NearZero(dval - dmax)) {
2377  dval = dmin;
2378  } else {
2379  continue;
2380  }
2381  curvept.Append(curvept[j]);
2382  curve_uvA.Append(uvA);
2383  curve_uvB.Append(uvB);
2384  }
2385  }
2386 }
2387 
2388 
2389 HIDDEN bool
2390 newton_ssi(double &uA, double &vA, double &uB, double &vB, const ON_Surface *surfA, const ON_Surface *surfB, double isect_tol)
2391 {
2392  ON_3dPoint pointA = surfA->PointAt(uA, vA);
2393  ON_3dPoint pointB = surfB->PointAt(uB, vB);
2394  if (pointA.DistanceTo(pointB) < isect_tol) {
2395  return true;
2396  }
2397 
2398  // Equations:
2399  // x_a(uA, vA) - x_b(uB, vB) = 0
2400  // y_a(uA, vA) - y_b(uB, vB) = 0
2401  // z_a(uA, vA) - z_b(uB, vB) = 0
2402  // It's an under-determined system. We use Moore-Penrose pseudoinverse:
2403  // pinv(A) = transpose(A) * inv(A * transpose(A)) (A is a 3x4 Jacobian)
2404  // A * pinv(A) = I_3
2405  double last_uA = DBL_MAX / 4, last_vA = DBL_MAX / 4, last_uB = DBL_MAX / 4, last_vB = DBL_MAX / 4;
2406 
2407  int iteration = 0;
2408  while (fabs(last_uA - uA) + fabs(last_vA - vA) + fabs(last_uB - uB) + fabs(last_vB - vB) > ON_ZERO_TOLERANCE
2409  && iteration++ < SSI_MAX_ITERATIONS) {
2410  last_uA = uA, last_vA = vA, last_uB = uB, last_vB = vB;
2411  ON_3dVector deriv_uA, deriv_vA, deriv_uB, deriv_vB;
2412  surfA->Ev1Der(uA, vA, pointA, deriv_uA, deriv_vA);
2413  surfB->Ev1Der(uB, vB, pointB, deriv_uB, deriv_vB);
2414  ON_Matrix j(3, 4), f(3, 1);
2415  for (int i = 0; i < 3; i++) {
2416  j[i][0] = deriv_uA[i];
2417  j[i][1] = deriv_vA[i];
2418  j[i][2] = -deriv_uB[i];
2419  j[i][3] = -deriv_vB[i];
2420  f[i][0] = pointA[i] - pointB[i];
2421  }
2422  ON_Matrix jtrans = j;
2423  jtrans.Transpose();
2424  ON_Matrix j_jtrans;
2425  j_jtrans.Multiply(j, jtrans);
2426  if (!j_jtrans.Invert(0.0)) {
2427  break;
2428  }
2429  ON_Matrix pinv_j;
2430  pinv_j.Multiply(jtrans, j_jtrans);
2431  ON_Matrix delta;
2432  delta.Multiply(pinv_j, f);
2433  uA -= delta[0][0];
2434  vA -= delta[1][0];
2435  uB -= delta[2][0];
2436  vB -= delta[3][0];
2437  }
2438 
2439  // make sure the solution is inside the domains
2440  uA = std::min(uA, surfA->Domain(0).Max());
2441  uA = std::max(uA, surfA->Domain(0).Min());
2442  vA = std::min(vA, surfA->Domain(1).Max());
2443  vA = std::max(vA, surfA->Domain(1).Min());
2444  uB = std::min(uB, surfB->Domain(0).Max());
2445  uB = std::max(uB, surfB->Domain(0).Min());
2446  vB = std::min(vB, surfB->Domain(1).Max());
2447  vB = std::max(vB, surfB->Domain(1).Min());
2448 
2449  pointA = surfA->PointAt(uA, vA);
2450  pointB = surfB->PointAt(uB, vB);
2451  return (pointA.DistanceTo(pointB) < isect_tol) && !isnan(uA) && !isnan(vA) && !isnan(uB) & !isnan(vB);
2452 }
2453 
2454 
2455 HIDDEN ON_Curve *
2456 link_curves(ON_Curve *&c1, ON_Curve *&c2)
2457 {
2458  ON_NurbsCurve *nc1 = c1->NurbsCurve();
2459  ON_NurbsCurve *nc2 = c2->NurbsCurve();
2460  if (nc1 && nc2) {
2461  nc1->Append(*nc2);
2462  delete c1;
2463  delete c2;
2464  c1 = NULL;
2465  c2 = NULL;
2466  delete nc2;
2467  return nc1;
2468  } else if (nc1) {
2469  delete nc1;
2470  } else if (nc2) {
2471  delete nc2;
2472  }
2473  return NULL;
2474 }
2475 
2476 
2478  ON_Curve *m_curve3d, *m_curveA, *m_curveB;
2479 
2480  int m_from; // 0: iso-curve from A, 1: from B
2481  int m_dir; // 0: u iso-curve, 1: v iso-curve
2482  double m_fix; // the param fixed in the iso-curve
2483  // m_dir==0: m_fix==u; m_dir==1: m_fix==v
2484  double m_min; // minimum of the variable param in the iso-curve
2485  double m_max; // maximum of the variable param in the iso-curve
2486 
2488  {
2489  m_curve3d = m_curveA = m_curveB = NULL;
2490  }
2491 
2492  ON_2dPoint Get2DParam(double t)
2493  {
2494  return m_dir ? ON_2dPoint(t, m_fix) : ON_2dPoint(m_fix, t);
2495  }
2496 
2498  {
2499  m_curve3d = m_curveA = m_curveB = NULL;
2500  }
2501 
2503  {
2504  if (m_curve3d) {
2505  delete m_curve3d;
2506  }
2507  m_curve3d = NULL;
2508  if (m_curveA) {
2509  delete m_curveA;
2510  }
2511  m_curveA = NULL;
2512  if (m_curveB) {
2513  delete m_curveB;
2514  }
2515  m_curveB = NULL;
2516  }
2517 };
2518 
2519 
2521  ON_SSX_EVENT *m_event; // Use pointer to the ON_SSX_EVENT in case that
2522  // ~ON_SSX_EVENT() is called.
2523  ON_BoundingBox m_bboxA; // 2D surfA bounding box of the closed region.
2524  enum TYPE {
2526  outer = 1,
2527  inner = 2
2528  } m_type;
2529 
2530  std::vector<Overlapevent *> m_inside;
2531 
2533 
2534  Overlapevent(ON_SSX_EVENT *_e) : m_event(_e)
2535  {
2536  m_type = undefined;
2537  m_bboxA = _e->m_curveA->BoundingBox();
2538  }
2539 
2540  bool IsPointOnBoundary(const ON_2dPoint &_pt) const
2541  {
2542  ON_ClassArray<ON_PX_EVENT> px_event;
2543  return ON_Intersect(ON_3dPoint(_pt), *(m_event->m_curveA), px_event);
2544  }
2545 
2546  // Test whether a point is inside the closed region (including boundary).
2547  // Approach: Use line-curve intersections. If the line between this point
2548  // and a point outside the b-box has an odd number of intersections with
2549  // the boundary, it's considered inside, otherwise outside. (Tangent
2550  // intersections are counted too.)
2551  bool IsPointIn(const ON_2dPoint &_pt) const
2552  {
2553  if (!m_bboxA.IsPointIn(_pt)) {
2554  return false;
2555  }
2556  if (IsPointOnBoundary(_pt)) {
2557  return true;
2558  }
2559  ON_ClassArray<ON_PX_EVENT> px_event;
2560  // this point must be outside the closed region
2561  ON_2dPoint out = _pt + ON_2dVector(m_bboxA.Diagonal());
2562  ON_LineCurve linecurve(_pt, out);
2563  ON_SimpleArray<ON_X_EVENT> x_event;
2564  ON_Intersect(&linecurve, m_event->m_curveA, x_event);
2565  int count = x_event.Count();
2566  for (int i = 0; i < x_event.Count(); i++) {
2567  // Find tangent intersections.
2568  // What should we do if it's ccx_overlap?
2569  if (m_event->m_curveA->TangentAt(x_event[i].m_a[0]).IsParallelTo(linecurve.m_line.Direction())) {
2570  count++;
2571  }
2572  }
2573  return count % 2 ? true : false;
2574  }
2575 
2576  // return false if the point is on the boundary
2577  bool IsPointOut(const ON_2dPoint &_pt) const
2578  {
2579  return !IsPointIn(_pt);
2580  }
2581 
2582  // Test whether there are intersections between the box boundaries
2583  // and the boundary curve.
2584  // If there are not intersection, there can be 3 cases:
2585  // 1) The box is completely in.
2586  // 2) The box is completely out.
2587  // 3) The box completely contains the closed region bounded by that curve.
2588  bool IsBoxIntersected(const ON_2dPoint &_min, const ON_2dPoint &_max) const
2589  {
2590  ON_2dPoint corner[4] = {
2591  _min,
2592  ON_2dPoint(_max.x, _min.y),
2593  _max,
2594  ON_2dPoint(_min.x, _max.y)
2595  };
2596  ON_LineCurve linecurve[4] = {
2597  ON_LineCurve(corner[0], corner[1]),
2598  ON_LineCurve(corner[1], corner[2]),
2599  ON_LineCurve(corner[2], corner[3]),
2600  ON_LineCurve(corner[3], corner[0])
2601  };
2602  ON_SimpleArray<ON_X_EVENT> x_event[4];
2603  for (int i = 0; i < 4; i++) {
2604  if (ON_Intersect(&linecurve[i], m_event->m_curveA, x_event[i])) {
2605  return true;
2606  }
2607  }
2608  return false;
2609  }
2610 
2611  bool IsBoxCompletelyIn(const ON_2dPoint &_min, const ON_2dPoint &_max) const
2612  {
2613  // there should not be intersections between the box boundaries
2614  // and the boundary curve
2615  if (IsBoxIntersected(_min, _max)) {
2616  return false;
2617  }
2618  ON_2dPoint center = (_min + _max) * 0.5;
2619  return m_bboxA.Includes(ON_BoundingBox(_min, _max)) && IsPointIn(center);
2620  }
2621 
2622  bool IsBoxCompletelyOut(const ON_2dPoint &_min, const ON_2dPoint &_max) const
2623  {
2624  if (IsBoxIntersected(_min, _max)) {
2625  return false;
2626  }
2627  ON_2dPoint center = (_min + _max) * 0.5;
2628  return !IsPointIn(center) && !ON_BoundingBox(_min, _max).Includes(m_bboxA);
2629  }
2630 
2631  bool IsCurveCompletelyIn(const ON_Curve *_curve) const
2632  {
2633  ON_SimpleArray<ON_X_EVENT> x_event;
2634  return !ON_Intersect(m_event->m_curve3d, _curve, x_event) && IsPointIn(_curve->PointAtStart());
2635  }
2636 };
2637 
2638 
2639 HIDDEN bool
2641  ON_2dPoint surf1_pt,
2642  const ON_Surface *surf1,
2643  const ON_Surface *surf2,
2644  Subsurface *surf2_tree)
2645 {
2646  if (!surf1->Domain(0).Includes(surf1_pt.x) ||
2647  !surf1->Domain(1).Includes(surf1_pt.y))
2648  {
2649  return false;
2650  }
2651  bool surf1_pt_intersects_surf2, surfs_parallel_at_pt = false;
2652  ON_ClassArray<ON_PX_EVENT> px_event;
2653 
2654  surf1_pt_intersects_surf2 = ON_Intersect(surf1->PointAt(surf1_pt.x, surf1_pt.y),
2655  *surf2, px_event, 1.0e-5, 0, 0, surf2_tree);
2656 
2657  if (surf1_pt_intersects_surf2) {
2658  surfs_parallel_at_pt = surf1->NormalAt(surf1_pt.x, surf1_pt.y).IsParallelTo(
2659  surf2->NormalAt(px_event[0].m_b[0], px_event[0].m_b[1]));
2660  }
2661  return surf1_pt_intersects_surf2 && surfs_parallel_at_pt;
2662 }
2663 
2664 
2665 HIDDEN bool
2667  const ON_SimpleArray<Overlapevent> &overlap_events,
2668  const Subsurface *subA,
2669  double isect_tolA)
2670 {
2671  ON_2dPoint bbox_min(subA->m_u.Min(), subA->m_v.Min());
2672  ON_2dPoint bbox_max(subA->m_u.Max(), subA->m_v.Max());
2673 
2674  // shrink bbox slightly
2675  ON_2dPoint offset(isect_tolA, isect_tolA);
2676  bbox_min += offset;
2677  bbox_max -= offset;
2678 
2679  // the bbox is considered to be inside an overlap if it's inside
2680  // an Overlapevent's outer loop and outside any of its inner loops
2681  bool inside_overlap = false;
2682  for (int i = 0; i < overlap_events.Count() && !inside_overlap; ++i) {
2683  const Overlapevent *outerloop = &overlap_events[i];
2684  if (outerloop->m_type == Overlapevent::outer &&
2685  outerloop->IsBoxCompletelyIn(bbox_min, bbox_max))
2686  {
2687  inside_overlap = true;
2688  for (size_t j = 0; j < outerloop->m_inside.size(); ++j) {
2689  const Overlapevent *innerloop = outerloop->m_inside[j];
2690  if (innerloop->m_type == Overlapevent::inner &&
2691  !innerloop->IsBoxCompletelyOut(bbox_min, bbox_max))
2692  {
2693  inside_overlap = false;
2694  break;
2695  }
2696  }
2697  }
2698  }
2699  return inside_overlap;
2700 }
2701 
2702 
2703 HIDDEN bool
2705  const ON_SimpleArray<Overlapevent> &overlap_events,
2706  const ON_2dPoint &pt)
2707 {
2708  // the pt is considered to be inside an overlap if it's inside
2709  // an Overlapevent's outer loop and outside any of its inner loops
2710  bool inside_overlap = false;
2711  for (int i = 0; i < overlap_events.Count() && !inside_overlap; ++i) {
2712  const Overlapevent *outerloop = &overlap_events[i];
2713  if (outerloop->m_type == Overlapevent::outer &&
2714  outerloop->IsPointIn(pt))
2715  {
2716  inside_overlap = true;
2717  for (size_t j = 0; j < outerloop->m_inside.size(); ++j) {
2718  const Overlapevent *innerloop = outerloop->m_inside[j];
2719  if (innerloop->m_type == Overlapevent::inner &&
2720  innerloop->IsPointIn(pt) &&
2721  !innerloop->IsPointOnBoundary(pt))
2722  {
2723  inside_overlap = false;
2724  break;
2725  }
2726  }
2727  }
2728  }
2729  return inside_overlap;
2730 }
2731 
2732 
2733 enum {
2737 };
2738 
2739 
2740 HIDDEN int
2742  const ON_Surface *surf1,
2743  double iso_t,
2744  int iso_dir,
2745  double overlap_start,
2746  double overlap_end,
2747  const ON_Surface *surf2,
2748  Subsurface *surf2_tree)
2749 {
2750  double test_distance = 0.01;
2751 
2752  // TODO: more sample points
2753  double u1, v1, u2, v2;
2754  double midpt = (overlap_start + overlap_end) * 0.5;
2755  if (iso_dir == 0) {
2756  u1 = u2 = midpt;
2757  v1 = iso_t - test_distance;
2758  v2 = iso_t + test_distance;
2759  } else {
2760  u1 = iso_t - test_distance;
2761  u2 = iso_t + test_distance;
2762  v1 = v2 = midpt;
2763  }
2764 
2765  bool in1, in2;
2766  ON_ClassArray<ON_PX_EVENT> px_event1, px_event2;
2767 
2768  in1 = is_pt_in_surf_overlap(ON_2dPoint(u1, v1), surf1, surf2, surf2_tree);
2769  in2 = is_pt_in_surf_overlap(ON_2dPoint(u2, v2), surf1, surf2, surf2_tree);
2770 
2771  if (in1 && in2) {
2772  return INSIDE_OVERLAP;
2773  }
2774  if (!in1 && !in2) {
2775  return OUTSIDE_OVERLAP;
2776  }
2777  return ON_OVERLAP_BOUNDARY;
2778 }
2779 
2780 
2781 ON_ClassArray<ON_3dPointArray>
2783  const ON_SimpleArray<OverlapSegment *> &overlaps,
2784  const ON_Surface *surfA,
2785  const ON_Surface *surfB,
2786  Subsurface *treeA,
2787  Subsurface *treeB,
2788  double isect_tol,
2789  double isect_tolA,
2790  double isect_tolB)
2791 {
2792  int count_before = overlaps.Count();
2793 
2794  ON_ClassArray<ON_3dPointArray> params(count_before);
2795 
2796  // count must equal capacity for array copy to work as expected
2797  // when the result of the function is assigned
2798  params.SetCount(params.Capacity());
2799 
2800  for (int i = 0; i < count_before; i++) {
2801  // split the curves from the intersection points between them
2802  int count = params[i].Count();
2803  for (int j = i + 1; j < count_before; j++) {
2804  ON_SimpleArray<ON_X_EVENT> x_event;
2805  ON_Intersect(overlaps[i]->m_curve3d, overlaps[j]->m_curve3d, x_event, isect_tol);
2806  count += x_event.Count();
2807  for (int k = 0; k < x_event.Count(); k++) {
2808  ON_3dPoint param;
2809  ON_2dPoint uvA, uvB;
2810  ON_ClassArray<ON_PX_EVENT> psxA, psxB, pcxA, pcxB;
2811 
2812  ON_3dPoint start_i = x_event[k].m_A[0];
2813  ON_3dPoint start_j = x_event[k].m_B[0];
2814  ON_3dPoint end_i = x_event[k].m_A[1];
2815  ON_3dPoint end_j = x_event[k].m_B[1];
2816  double start_ti = x_event[k].m_a[0];
2817  double start_tj = x_event[k].m_b[0];
2818  double end_ti = x_event[k].m_a[1];
2819  double end_tj = x_event[k].m_b[1];
2820 
2821  // the curves may have opposite directions, so make
2822  // sure the ends are consistent
2823  if (x_event[k].m_type == ON_X_EVENT::ccx_overlap &&
2824  start_i.DistanceTo(start_j) > start_i.DistanceTo(end_j)) {
2825  std::swap(start_j, end_j);
2826  std::swap(start_tj, end_tj);
2827  }
2828  // ensure that the curve intersection point also
2829  // intersects both surfaces (i.e. the point is in the
2830  // overlap)
2831  if (ON_Intersect(start_i, *surfA, psxA, isect_tolA, 0, 0, treeA)
2832  && ON_Intersect(start_j, *surfB, psxB, isect_tolB, 0, 0, treeB)) {
2833  // pull the 3D curve back to the 2D space
2834  uvA = psxA[0].m_b;
2835  uvB = psxB[0].m_b;
2836 
2837  // Convert the uv points into the equivalent curve
2838  // parameters by intersecting them with the 2D
2839  // versions of the curves. The 2D intersection is
2840  // done in 3D; note that uvA/uvB are implicitly
2841  // converted to ON_3dPoint with z=0.0.
2842  if (ON_Intersect(uvA, *(overlaps[i]->m_curveA), pcxA, isect_tolA)
2843  && ON_Intersect(uvB, *(overlaps[i]->m_curveB), pcxB, isect_tolB)) {
2844  param.x = start_ti;
2845  param.y = pcxA[0].m_b[0];
2846  param.z = pcxB[0].m_b[0];
2847  params[i].Append(param);
2848  }
2849  pcxA.SetCount(0);
2850  pcxB.SetCount(0);
2851  if (ON_Intersect(uvA, *(overlaps[j]->m_curveA), pcxA, isect_tolA)
2852  && ON_Intersect(uvB, *(overlaps[j]->m_curveB), pcxB, isect_tolB)) {
2853  // the same routine for overlaps[j]
2854  param.x = start_tj;
2855  param.y = pcxA[0].m_b[0];
2856  param.z = pcxB[0].m_b[0];
2857  params[j].Append(param);
2858  }
2859  }
2860  if (x_event[k].m_type == ON_X_EVENT::ccx_overlap) {
2861  // for overlap, have to do the same with the intersection end points
2862  psxA.SetCount(0);
2863  psxB.SetCount(0);
2864  if (ON_Intersect(end_i, *surfA, psxA, isect_tol, 0, 0, treeA)
2865  && ON_Intersect(end_j, *surfB, psxB, isect_tol, 0, 0, treeB)) {
2866  // pull the 3D curve back to the 2D space
2867  uvA = psxA[0].m_b;
2868  uvB = psxB[0].m_b;
2869  pcxA.SetCount(0);
2870  pcxB.SetCount(0);
2871  if (ON_Intersect(uvA, *(overlaps[i]->m_curveA), pcxA, isect_tolA)
2872  && ON_Intersect(uvB, *(overlaps[i]->m_curveB), pcxB, isect_tolB)) {
2873  param.x = end_ti;
2874  param.y = pcxA[0].m_b[0];
2875  param.z = pcxB[0].m_b[0];
2876  params[i].Append(param);
2877  }
2878  pcxA.SetCount(0);
2879  pcxB.SetCount(0);
2880  if (ON_Intersect(uvA, *(overlaps[j]->m_curveA), pcxA, isect_tolA)
2881  && ON_Intersect(uvB, *(overlaps[j]->m_curveB), pcxB, isect_tolB)) {
2882  // the same routine for overlaps[j]
2883  param.x = end_tj;
2884  param.y = pcxA[0].m_b[0];
2885  param.z = pcxB[0].m_b[0];
2886  params[j].Append(param);
2887  }
2888  }
2889  }
2890  }
2891  }
2892  }
2893 
2894  return params;
2895 }
2896 
2897 
2898 HIDDEN bool
2900 {
2901  if (overlap &&
2902  overlap->m_curve3d &&
2903  overlap->m_curveA &&
2904  overlap->m_curveB) {
2905  return true;
2906  }
2907  return false;
2908 }
2909 
2910 
2911 HIDDEN void
2913  ON_SimpleArray<OverlapSegment *> &overlaps,
2914  const ON_Surface *surfA,
2915  const ON_Surface *surfB,
2916  Subsurface *treeA,
2917  Subsurface *treeB,
2918  double isect_tol,
2919  double isect_tolA,
2920  double isect_tolB)
2921 {
2922  // Not points, but a bundle of params:
2923  // params[i] - the separating params on overlaps[i]
2924  // x - curve3d, y - curveA, z - curveB
2925  ON_ClassArray<ON_3dPointArray> params = get_overlap_intersection_parameters(
2926  overlaps, surfA, surfB, treeA, treeB, isect_tol,
2927  isect_tolA, isect_tolB);
2928 
2929  // split overlap curves at intersection parameters into subcurves
2930  for (int i = 0; i < params.Count(); i++) {
2931  params[i].QuickSort(ON_CompareIncreasing);
2932  int start = 0;
2933  bool splitted = false;
2934  for (int j = 1; j < params[i].Count(); j++) {
2935  if (params[i][j].x - params[i][start].x < isect_tol) {
2936  continue;
2937  }
2938  ON_Curve *subcurveA = NULL;
2939  try {
2940  subcurveA = sub_curve(overlaps[i]->m_curveA,
2941  params[i][start].y, params[i][j].y);
2942  } catch (InvalidInterval &) {
2943  continue;
2944  }
2945  bool isvalid = false, isreversed = false;
2946  double test_distance = 0.01;
2947  // TODO: more sample points
2948  ON_2dPoint uv1, uv2;
2949  uv1 = uv2 = subcurveA->PointAt(subcurveA->Domain().Mid());
2950  ON_3dVector normal = ON_CrossProduct(subcurveA->TangentAt(subcurveA->Domain().Mid()), ON_3dVector::ZAxis);
2951  normal.Unitize();
2952  uv1 -= normal * test_distance; // left
2953  uv2 += normal * test_distance; // right
2954  bool in1 = is_pt_in_surf_overlap(uv1, surfA, surfB, treeB);
2955  bool in2 = is_pt_in_surf_overlap(uv2, surfA, surfB, treeB);
2956  if (in1 && !in2) {
2957  isvalid = true;
2958  } else if (!in1 && in2) {
2959  // the right side is overlapped
2960  isvalid = true;
2961  isreversed = true;
2962  }
2963  if (isvalid) {
2965  seg->m_curve3d = sub_curve(overlaps[i]->m_curve3d, params[i][start].x, params[i][j].x);
2966  seg->m_curveA = subcurveA;
2967  seg->m_curveB = sub_curve(overlaps[i]->m_curveB, params[i][start].z, params[i][j].z);
2968  if (isreversed) {
2969  seg->m_curve3d->Reverse();
2970  seg->m_curveA->Reverse();
2971  seg->m_curveB->Reverse();
2972  }
2973  seg->m_dir = overlaps[i]->m_dir;
2974  seg->m_fix = overlaps[i]->m_fix;
2975  overlaps.Append(seg);
2976  } else {
2977  delete subcurveA;
2978  }
2979  start = j;
2980  splitted = true;
2981  }
2982  if (splitted) {
2983  delete overlaps[i];
2984  overlaps[i] = NULL;
2985  }
2986  }
2987 
2988  for (int i = 0; i < overlaps.Count(); i++) {
2989  for (int j = i + 1; j < overlaps.Count(); j++) {
2990  if (!is_valid_overlap(overlaps[i]) || !is_valid_overlap(overlaps[j])) {
2991  continue;
2992  }
2993  // eliminate duplications
2994  ON_SimpleArray<ON_X_EVENT> x_event1, x_event2;
2995  if (ON_Intersect(overlaps[i]->m_curveA, overlaps[j]->m_curveA, x_event1, isect_tolA)
2996  && ON_Intersect(overlaps[i]->m_curveB, overlaps[j]->m_curveB, x_event2, isect_tolB)) {
2997  if (x_event1[0].m_type == ON_X_EVENT::ccx_overlap
2998  && x_event2[0].m_type == ON_X_EVENT::ccx_overlap) {
2999  if (ON_NearZero(x_event1[0].m_a[0] - overlaps[i]->m_curveA->Domain().Min())
3000  && ON_NearZero(x_event1[0].m_a[1] - overlaps[i]->m_curveA->Domain().Max())
3001  && ON_NearZero(x_event2[0].m_a[0] - overlaps[i]->m_curveB->Domain().Min())
3002  && ON_NearZero(x_event2[0].m_a[1] - overlaps[i]->m_curveB->Domain().Max())) {
3003  // overlaps[i] is completely inside overlaps[j]
3004  delete overlaps[i];
3005  overlaps[i] = NULL;
3006  } else if (ON_NearZero(x_event1[0].m_b[0] - overlaps[j]->m_curveA->Domain().Min())
3007  && ON_NearZero(x_event1[0].m_b[1] - overlaps[j]->m_curveA->Domain().Max())
3008  && ON_NearZero(x_event2[0].m_b[0] - overlaps[j]->m_curveB->Domain().Min())
3009  && ON_NearZero(x_event2[0].m_b[1] - overlaps[j]->m_curveB->Domain().Max())) {
3010  // overlaps[j] is completely inside overlaps[i]
3011  delete overlaps[j];
3012  overlaps[j] = NULL;
3013  }
3014  }
3015  }
3016  }
3017  }
3018 }
3019 
3020 
3021 HIDDEN std::vector<Subsurface *>
3022 subdivide_subsurface(Subsurface *sub)
3023 {
3024  std::vector<Subsurface *> parts;
3025  if (sub->Split() != 0) {
3026  parts.push_back(sub);
3027  } else {
3028  for (int i = 0; i < 4; ++i) {
3029  parts.push_back(sub->m_children[i]);
3030  }
3031  }
3032  return parts;
3033 }
3034 
3035 
3036 HIDDEN std::pair<Triangle, Triangle>
3037 get_subsurface_triangles(const Subsurface *sub, const ON_Surface *surf)
3038 {
3039  double min_u = sub->m_u.Min();
3040  double max_u = sub->m_u.Max();
3041  double min_v = sub->m_v.Min();
3042  double max_v = sub->m_v.Max();
3043 
3044  ON_2dPoint corners2d[4] = {
3045  ON_2dPoint(min_u, min_v),
3046  ON_2dPoint(min_u, max_v),
3047  ON_2dPoint(max_u, min_v),
3048  ON_2dPoint(max_u, max_v)};
3049 
3050  ON_3dPoint corners3d[4];
3051  for (int i = 0; i < 4; ++i) {
3052  corners3d[i] = surf->PointAt(corners2d[i].x, corners2d[i].y);
3053  }
3054 
3055  std::pair<Triangle, Triangle> triangles;
3056  triangles.first.a_2d = corners2d[0];
3057  triangles.first.b_2d = corners2d[1];
3058  triangles.first.c_2d = corners2d[2];
3059  triangles.first.a = corners3d[0];
3060  triangles.first.b = corners3d[1];
3061  triangles.first.c = corners3d[2];
3062  triangles.second.a_2d = corners2d[1];
3063  triangles.second.b_2d = corners2d[2];
3064  triangles.second.c_2d = corners2d[3];
3065  triangles.second.a = corners3d[1];
3066  triangles.second.b = corners3d[2];
3067  triangles.second.c = corners3d[3];
3068 
3069  return triangles;
3070 }
3071 
3072 
3073 HIDDEN ON_2dPoint
3074 barycentric_to_uv(const ON_3dPoint &bc, const Triangle &tri)
3075 {
3076  ON_3dVector vertx(tri.a_2d.x, tri.b_2d.x, tri.c_2d.x);
3077  ON_3dVector verty(tri.a_2d.y, tri.b_2d.y, tri.c_2d.y);
3078  ON_2dPoint uv(ON_DotProduct(bc, vertx), ON_DotProduct(bc, verty));
3079  return uv;
3080 }
3081 
3082 
3083 HIDDEN bool
3084 join_continuous_ppair_to_polyline(ON_SimpleArray<int> *polyline, const PointPair &ppair)
3085 {
3086  int polyline_first = (*polyline)[0];
3087  int polyline_last = (*polyline)[polyline->Count() - 1];
3088 
3089  if (polyline_first == ppair.indexA) {
3090  polyline->Insert(0, ppair.indexB);
3091  return true;
3092  }
3093  if (polyline_first == ppair.indexB) {
3094  polyline->Insert(0, ppair.indexA);
3095  return true;
3096  }
3097  if (polyline_last == ppair.indexA) {
3098  polyline->Append(ppair.indexB);
3099  return true;
3100  }
3101  if (polyline_last == ppair.indexB) {
3102  polyline->Append(ppair.indexA);
3103  return true;
3104  }
3105  return false;
3106 }
3107 
3108 
3109 // Algorithm Overview
3110 //
3111 // 1) Find overlap intersections (regions where the two surfaces are
3112 // flush to one another).
3113 // a) Intersect the isocurves of surfA with surfB and vice versa
3114 // (any intersection points found are saved for later).
3115 // b) We save any overlap intersection curves that appear to be
3116 // boundary curves (curves which contain part of the boundary
3117 // of the overlap region).
3118 // i) If the overlap intersection curve is at the edge of the
3119 // surface, it's automatically a boundary curve.
3120 // ii) If the surfaces overlap on only one side of the overlap
3121 // intersection curve, it's a boundary curve (or at least
3122 // contains a boundary curve).
3123 // c) Because our boundary test is rather simple, not all of the
3124 // saved overlap curves actually overlap the surfaces over their
3125 // whole domain, so we do additional processing to get the real
3126 // overlaps.
3127 // i) The overlap curves are intersected with each other.
3128 // ii) The curves are split into subcurves at the points where
3129 // they intersect.
3130 // iii) We retest the subcurves to see if they are truly part of
3131 // the overlap region boundary, discarding those that aren't.
3132 // iv) Any overlap curve which is completely contained by another
3133 // overlap curve is discarded, hopefully leaving a set of
3134 // genuine and unique overlap curves.
3135 // d) Overlap curves that share endpoints are merged together.
3136 // e) Any of the resulting overlap curves which don't form closed
3137 // loops are discarded.
3138 // f) The overlap loops are turned into overlap intersection
3139 // events, with the orientation of outer and inner loops
3140 // adjusted as needed to ensure the overlapping regions are
3141 // correctly defined.
3142 // 2) Find non-overlap intersections.
3143 // a) The two surfaces are repeatedly subdivided into subsurfaces
3144 // until the MAX_SSI_DEPTH is reached. The subsurface pairs
3145 // whose bounding boxes intersect are saved.
3146 // b) Each pair of potentially intersecting subsurfaces is
3147 // intersected.
3148 // i) Each subsurface is approximated by two triangles which
3149 // share an edge.
3150 // ii) The triangles are intersected.
3151 // iii) The Newton solver is used to get the center point of the
3152 // intersection, using the averaged center point of the
3153 // triangle intersections as the initial guess.
3154 // c) The center points of the subsurface intersections are
3155 // combined with the isocurve/surface intersection points found
3156 // in step 1) above.
3157 // d) For surface seams that are closed, additional points are
3158 // generated on the seams (3d points near the seams are copied,
3159 // and their uv values are adjusted to put them over the seam).
3160 // e) The list of intersection points is reduced to just the unique
3161 // intersection points that are not inside any overlap regions.
3162 // f) Polylines are created from the intersection points.
3163 // i) Each point starts as its own polyline, and serves as both
3164 // endpoints of the line.
3165 // ii) The closest pair of endpoints are joined together, then
3166 // the next closest, and so on, growing the length of the
3167 // polylines.
3168 // iv) Endpoints that are too far apart (distance > max_dist) are
3169 // not joined, and some intersection points may remain by
3170 // themselves.
3171 // g) Intersection curves are generated from the polylines, and are
3172 // used to make intersection events.
3173 // h) Any single intersection points found in step 2-f-iv) above
3174 // that aren't contained by one of the generated intersection
3175 // curves are used to make additional point intersection
3176 // events.
3177 int
3178 ON_Intersect(const ON_Surface *surfA,
3179  const ON_Surface *surfB,
3180  ON_ClassArray<ON_SSX_EVENT> &x,
3181  double isect_tol,
3182  double overlap_tol,
3183  double fitting_tol,
3184  const ON_Interval *surfaceA_udomain,
3185  const ON_Interval *surfaceA_vdomain,
3186  const ON_Interval *surfaceB_udomain,
3187  const ON_Interval *surfaceB_vdomain,
3188  Subsurface *treeA,
3189  Subsurface *treeB)
3190 {
3191  if (surfA == NULL || surfB == NULL) {
3192  return 0;
3193  }
3194 
3195  int original_count = x.Count();
3196 
3197  if (isect_tol <= 0.0) {
3198  isect_tol = SSI_DEFAULT_TOLERANCE;
3199  }
3200  if (overlap_tol < isect_tol) {
3201  overlap_tol = 2.0 * isect_tol;
3202  }
3203  V_MAX(fitting_tol, isect_tol);
3204 
3205  check_domain(surfaceA_udomain, surfA->Domain(0), "surfaceA_udomain");
3206  check_domain(surfaceA_vdomain, surfA->Domain(1), "surfaceA_vdomain");
3207  check_domain(surfaceB_udomain, surfB->Domain(0), "surfaceB_udomain");
3208  check_domain(surfaceB_vdomain, surfB->Domain(1), "surfaceB_vdomain");
3209 
3210  /* First step: Initialize the first two Subsurface.
3211  * It's just like getting the root of the surface tree.
3212  */
3213  Subsurface *rootA = treeA;
3214  Subsurface *rootB = treeB;
3215  try {
3216  if (rootA == NULL) {
3217  rootA = get_surface_root(surfA, surfaceA_udomain,
3218  surfaceA_vdomain);
3219  }
3220  if (rootB == NULL) {
3221  rootB = get_surface_root(surfB, surfaceB_udomain,
3222  surfaceB_vdomain);
3223  }
3224  } catch (std::bad_alloc &e) {
3225  bu_log("%s", e.what());
3226  if (treeA == NULL) {
3227  delete rootA;
3228  }
3229  return 0;
3230  }
3231 
3232  if (!rootA->Intersect(*rootB, isect_tol)) {
3233  if (treeA == NULL) {
3234  delete rootA;
3235  }
3236  if (treeB == NULL) {
3237  delete rootB;
3238  }
3239  return 0;
3240  }
3241 
3242  double isect_tolA =
3243  tolerance_2d_from_3d(isect_tol, rootA, surfA,
3244  surfaceA_udomain, surfaceA_vdomain);
3245 
3246  double fitting_tolA = tolerance_2d_from_3d(fitting_tol,
3247  rootA, surfA, surfaceA_udomain, surfaceA_vdomain);
3248 
3249  double isect_tolB =
3250  tolerance_2d_from_3d(isect_tol, rootB, surfB,
3251  surfaceB_udomain, surfaceB_vdomain);
3252 
3253  double fitting_tolB = tolerance_2d_from_3d(fitting_tol,
3254  rootB, surfB, surfaceB_udomain, surfaceB_vdomain);
3255 
3256  ON_3dPointArray curvept, tmp_curvept;
3257  ON_2dPointArray curve_uvA, curve_uvB, tmp_curve_uvA, tmp_curve_uvB;
3258 
3259  // Overlap detection:
3260  // According to the Theorem 3 in paper:
3261  // http://libgen.org/scimag1/10.1016/S0010-4485%252896%252900099-1.pdf
3262  // For two Bezier patches, if they overlap over a region, then the overlap
3263  // region must be bounded by parts of the boundaries of the two patches.
3264  // The Bezier patches for a NURBS surface are always bounded at the knots.
3265  // In other words, given two pairs of neighbor knots (in two dimensions),
3266  // we can get a Bezier patch between them.
3267  // (See ON_NurbsSurface::ConvertSpanToBezier()).
3268  // So we actually don't need to generate the Bezier patches explicitly,
3269  // and we can get the boundaries of them using IsoCurve() on knots.
3270  // Deal with boundaries with curve-surface intersections.
3271 
3272  ON_SimpleArray<OverlapSegment *> overlaps;
3273  for (int i = 0; i < 4; i++) {
3274  const ON_Surface *surf1 = i >= 2 ? surfB : surfA;
3275  const ON_Surface *surf2 = i >= 2 ? surfA : surfB;
3276  Subsurface *tree2 = i >= 2 ? treeA : treeB;
3277  ON_2dPointArray &ptarray1 = i >= 2 ? tmp_curve_uvB : tmp_curve_uvA;
3278  ON_2dPointArray &ptarray2 = i >= 2 ? tmp_curve_uvA : tmp_curve_uvB;
3279  int knot_dir = i % 2;
3280  int surf_dir = 1 - knot_dir;
3281  int knot_count = surf1->SpanCount(surf_dir) + 1;
3282  double *surf1_knots = new double[knot_count];
3283  surf1->GetSpanVector(surf_dir, surf1_knots);
3284  // knots that can be boundaries of Bezier patches
3285  ON_SimpleArray<double> surf1_bknots;
3286  surf1_bknots.Append(surf1_knots[0]);
3287  for (int j = 1; j < knot_count; j++) {
3288  if (surf1_knots[j] > *(surf1_bknots.Last())) {
3289  surf1_bknots.Append(surf1_knots[j]);
3290  }
3291  }
3292  if (surf1->IsClosed(surf_dir)) {
3293  surf1_bknots.Remove();
3294  }
3295 
3296  for (int j = 0; j < surf1_bknots.Count(); j++) {
3297  double surf1_knot = surf1_bknots[j];
3298  ON_Curve *surf1_boundary_iso = surf1->IsoCurve(knot_dir, surf1_knot);
3299  ON_SimpleArray<ON_X_EVENT> x_event;
3300  ON_CurveArray overlap2d;
3301  ON_Intersect(surf1_boundary_iso, surf2, x_event, isect_tol,
3302  overlap_tol, 0, 0, 0, &overlap2d);
3303 
3304  // stash surf1 points and surf2 parameters
3305  for (int k = 0; k < x_event.Count(); k++) {
3306  ON_X_EVENT &event = x_event[k];
3307 
3308  tmp_curvept.Append(event.m_A[0]);
3309  ptarray2.Append(ON_2dPoint(event.m_b[0], event.m_b[1]));
3310 
3311  if (event.m_type == ON_X_EVENT::csx_overlap) {
3312  tmp_curvept.Append(event.m_A[1]);
3313  ptarray2.Append(ON_2dPoint(event.m_b[2], event.m_b[3]));
3314  }
3315  }
3316 
3317  for (int k = 0; k < x_event.Count(); k++) {
3318  ON_X_EVENT &event = x_event[k];
3319 
3320  ON_2dPoint iso_pt1;
3321  iso_pt1.x = knot_dir ? surf1_knot : event.m_a[0];
3322  iso_pt1.y = knot_dir ? event.m_a[0] : surf1_knot;
3323  ptarray1.Append(iso_pt1);
3324 
3325  if (event.m_type == ON_X_EVENT::csx_overlap) {
3326  ON_2dPoint iso_pt2;
3327  iso_pt2.x = knot_dir ? surf1_knot : event.m_a[1];
3328  iso_pt2.y = knot_dir ? event.m_a[1] : surf1_knot;
3329  ptarray1.Append(iso_pt2);
3330  // get the overlap curve
3331  if (event.m_a[0] < event.m_a[1]) {
3332  bool curve_on_overlap_boundary = false;
3333  if (j == 0 || (j == surf1_bknots.Count() - 1 && !surf1->IsClosed(surf_dir))) {
3334  curve_on_overlap_boundary = true;
3335  } else {
3336  double overlap_start = event.m_a[0];
3337  double overlap_end = event.m_a[1];
3338  int location = isocurve_surface_overlap_location(
3339  surf1, surf1_knot, knot_dir, overlap_start, overlap_end, surf2, tree2);
3340  curve_on_overlap_boundary = (location == ON_OVERLAP_BOUNDARY);
3341  }
3342  if (curve_on_overlap_boundary) {
3343  // one side of it is shared and the other side is non-shared
3345  try {
3346  seg->m_curve3d = sub_curve(surf1_boundary_iso,
3347  event.m_a[0], event.m_a[1]);
3348  if (i < 2) {
3349  seg->m_curveA = new ON_LineCurve(iso_pt1, iso_pt2);
3350  seg->m_curveB = overlap2d[k];
3351  } else {
3352  seg->m_curveB = new ON_LineCurve(iso_pt1, iso_pt2);
3353  seg->m_curveA = overlap2d[k];
3354  }
3355  seg->m_dir = surf_dir;
3356  seg->m_fix = surf1_knot;
3357  overlaps.Append(seg);
3358  if (j == 0 && surf1->IsClosed(surf_dir)) {
3359  // Something like close_domain().
3360  // If the domain is closed, the iso-curve on the
3361  // first knot and the last knot is the same, so
3362  // we don't need to compute the intersections twice.
3363  seg = new OverlapSegment;
3364  iso_pt1.x = knot_dir ?
3365  surf1_knots[knot_count - 1] : event.m_a[0];
3366  iso_pt1.y = knot_dir ?
3367  event.m_a[0] : surf1_knots[knot_count - 1];
3368  iso_pt2.x = knot_dir ?
3369  surf1_knots[knot_count - 1] : event.m_a[1];
3370  iso_pt2.y = knot_dir ?
3371  event.m_a[1] : surf1_knots[knot_count - 1];
3372  seg->m_curve3d = (*overlaps.Last())->m_curve3d->Duplicate();
3373  if (i < 2) {
3374  seg->m_curveA = new ON_LineCurve(iso_pt1, iso_pt2);
3375  seg->m_curveB = overlap2d[k]->Duplicate();
3376  } else {
3377  seg->m_curveB = new ON_LineCurve(iso_pt1, iso_pt2);
3378  seg->m_curveA = overlap2d[k]->Duplicate();
3379  }
3380  seg->m_dir = surf_dir;
3381  seg->m_fix = surf1_knots[knot_count - 1];
3382  overlaps.Append(seg);
3383  }
3384  // We set overlap2d[k] to NULL in case the curve
3385  // is delete by the destructor of overlap2d. (See ~ON_CurveArray())
3386  overlap2d[k] = NULL;
3387  } catch (InvalidInterval &e) {
3388  bu_log("%s", e.what());
3389  delete seg;
3390  }
3391  }
3392  }
3393  }
3394  }
3395  delete surf1_boundary_iso;
3396  }
3397  delete []surf1_knots;
3398  }
3399 
3400  split_overlaps_at_intersections(overlaps, surfA, surfB, treeA, treeB,
3401  isect_tol, isect_tolA, isect_tolB);
3402 
3403  // find the neighbors for every overlap segment
3404  ON_SimpleArray<bool> start_linked(overlaps.Count()), end_linked(overlaps.Count());
3405  for (int i = 0; i < overlaps.Count(); i++) {
3406  start_linked[i] = end_linked[i] = false;
3407  }
3408  for (int i = 0; i < overlaps.Count(); i++) {
3409  if (!is_valid_overlap(overlaps[i])) {
3410  continue;
3411  }
3412 
3413  if (overlaps[i]->m_curve3d->IsClosed() && overlaps[i]->m_curveA->IsClosed() && overlaps[i]->m_curveB->IsClosed()) {
3414  start_linked[i] = end_linked[i] = true;
3415  }
3416 
3417  for (int j = i + 1; j < overlaps.Count(); j++) {
3418  if (!is_valid_overlap(overlaps[j])) {
3419  continue;
3420  }
3421 
3422  // check whether the start and end points are linked to
3423  // another curve
3424 #define OVERLAPS_LINKED(from, to) \
3425  overlaps[i]->m_curve3d->PointAt##from().DistanceTo(overlaps[j]->m_curve3d->PointAt##to()) < isect_tol && \
3426  overlaps[i]->m_curveA->PointAt##from().DistanceTo(overlaps[j]->m_curveA->PointAt##to()) < isect_tolA && \
3427  overlaps[i]->m_curveB->PointAt##from().DistanceTo(overlaps[j]->m_curveB->PointAt##to()) < isect_tolB
3428 
3429  if (OVERLAPS_LINKED(Start, End)) {
3430  start_linked[i] = end_linked[j] = true;
3431  } else if (OVERLAPS_LINKED(End, Start)) {
3432  start_linked[j] = end_linked[i] = true;
3433  } else if (OVERLAPS_LINKED(Start, Start)) {
3434  start_linked[i] = start_linked[j] = true;
3435  } else if (OVERLAPS_LINKED(End, End)) {
3436  end_linked[i] = end_linked[j] = true;
3437  }
3438  }
3439  }
3440 
3441  for (int i = 0; i < overlaps.Count(); i++) {
3442  if (!is_valid_overlap(overlaps[i])) {
3443  continue;
3444  }
3445  if (!start_linked[i] || !end_linked[i]) {
3446  delete overlaps[i];
3447  overlaps[i] = NULL;
3448  }
3449  }
3450 
3451  for (int i = 0; i < overlaps.Count(); i++) {
3452  if (!is_valid_overlap(overlaps[i])) {
3453  continue;
3454  }
3455 
3456  for (int j = 0; j <= overlaps.Count(); j++) {
3457  if (overlaps[i]->m_curve3d->IsClosed() && overlaps[i]->m_curveA->IsClosed() && overlaps[i]->m_curveB->IsClosed()) {
3458  // i-th curve is a closed loop, completely bounding an overlap region
3459  ON_SSX_EVENT event;
3460  event.m_curve3d = curve_fitting(overlaps[i]->m_curve3d, fitting_tol);
3461  event.m_curveA = curve_fitting(overlaps[i]->m_curveA, fitting_tolA);
3462  event.m_curveB = curve_fitting(overlaps[i]->m_curveB, fitting_tolB);
3463  event.m_type = ON_SSX_EVENT::ssx_overlap;
3464  // normalize the curves
3465  event.m_curve3d->SetDomain(ON_Interval(0.0, 1.0));
3466  event.m_curveA->SetDomain(ON_Interval(0.0, 1.0));
3467  event.m_curveB->SetDomain(ON_Interval(0.0, 1.0));
3468  x.Append(event);
3469  // Set the curves to NULL in case they are deleted by
3470  // ~ON_SSX_EVENT() or ~ON_CurveArray().
3471  event.m_curve3d = event.m_curveA = event.m_curveB = NULL;
3472  overlaps[i]->SetCurvesToNull();
3473  delete overlaps[i];
3474  overlaps[i] = NULL;
3475  break;
3476  }
3477 
3478  if (j == overlaps.Count() || j == i || !is_valid_overlap(overlaps[j])) {
3479  continue;
3480  }
3481 
3482  // merge the curves that link together
3483  if (OVERLAPS_LINKED(Start, End)) {
3484  overlaps[i]->m_curve3d = link_curves(overlaps[j]->m_curve3d, overlaps[i]->m_curve3d);
3485  overlaps[i]->m_curveA = link_curves(overlaps[j]->m_curveA, overlaps[i]->m_curveA);
3486  overlaps[i]->m_curveB = link_curves(overlaps[j]->m_curveB, overlaps[i]->m_curveB);
3487  } else if (OVERLAPS_LINKED(End, Start)) {
3488  overlaps[i]->m_curve3d = link_curves(overlaps[i]->m_curve3d, overlaps[j]->m_curve3d);
3489  overlaps[i]->m_curveA = link_curves(overlaps[i]->m_curveA, overlaps[j]->m_curveA);
3490  overlaps[i]->m_curveB = link_curves(overlaps[i]->m_curveB, overlaps[j]->m_curveB);
3491  } else if (OVERLAPS_LINKED(Start, Start)) {
3492  if (overlaps[i]->m_curve3d->Reverse() && overlaps[i]->m_curveA->Reverse() && overlaps[i]->m_curveB->Reverse()) {
3493  overlaps[i]->m_curve3d = link_curves(overlaps[i]->m_curve3d, overlaps[j]->m_curve3d);
3494  overlaps[i]->m_curveA = link_curves(overlaps[i]->m_curveA, overlaps[j]->m_curveA);
3495  overlaps[i]->m_curveB = link_curves(overlaps[i]->m_curveB, overlaps[j]->m_curveB);
3496  }
3497  } else if (OVERLAPS_LINKED(End, End)) {
3498  if (overlaps[j]->m_curve3d->Reverse() && overlaps[j]->m_curveA->Reverse() && overlaps[j]->m_curveB->Reverse()) {
3499  overlaps[i]->m_curve3d = link_curves(overlaps[i]->m_curve3d, overlaps[j]->m_curve3d);
3500  overlaps[i]->m_curveA = link_curves(overlaps[i]->m_curveA, overlaps[j]->m_curveA);
3501  overlaps[i]->m_curveB = link_curves(overlaps[i]->m_curveB, overlaps[j]->m_curveB);
3502  }
3503  }
3504  if (!is_valid_overlap(overlaps[j])) {
3505  delete overlaps[j];
3506  overlaps[j] = NULL;
3507  }
3508  }
3509  }
3510 
3511  // generate Overlapevents
3512  ON_SimpleArray<Overlapevent> overlap_events;
3513  for (int i = original_count; i < x.Count(); i++) {
3514  overlap_events.Append(Overlapevent(&x[i]));
3515  }
3516 
3517  for (int i = original_count; i < x.Count(); i++) {
3518  // The overlap region should be to the LEFT of that *m_curveA*.
3519  // (see opennurbs/opennurbs_x.h)
3520 
3521  // find a point with a single tangent
3522  double start_t = x[i].m_curveA->Domain().Mid();
3523  if (!x[i].m_curveA->IsContinuous(ON::G1_continuous, start_t)) {
3524  // if mid point doesn't work, try a couple other points
3525  // TODO: what happens if none of them work?
3526  start_t = x[i].m_curveA->Domain().NormalizedParameterAt(1.0 / 3.0);
3527  if (!x[i].m_curveA->IsContinuous(ON::G1_continuous, start_t)) {
3528  start_t = x[i].m_curveA->Domain().NormalizedParameterAt(2.0 / 3.0);
3529  }
3530  }
3531 
3532  ON_3dPoint start_ptA = x[i].m_curveA->PointAt(start_t);
3533  ON_3dVector normA = ON_CrossProduct(ON_3dVector::ZAxis,
3534  x[i].m_curveA->TangentAt(start_t));
3535  double line_len = 1.0 + x[i].m_curveA->BoundingBox().Diagonal().Length();
3536 
3537  // get a line that should be extending left, through the loop
3538  ON_3dPoint left_ptA = start_ptA + normA * line_len;
3539  ON_LineCurve inside_lineA(start_ptA, left_ptA);
3540 
3541  // get a line that should be extending right, away from the loop
3542  ON_3dPoint right_ptA = start_ptA - normA * line_len;
3543  ON_LineCurve outside_lineA(start_ptA, right_ptA);
3544 
3545  // outside line should intersect at start only
3546  // inside line should intersect at start and exit
3547  // otherwise we've got them backwards
3548  // TODO: seems like this should really be an even/odd test
3549  ON_SimpleArray<ON_X_EVENT> inside_events, outside_events;
3550  ON_Intersect(&inside_lineA, x[i].m_curveA, inside_events, isect_tol);
3551  ON_Intersect(&outside_lineA, x[i].m_curveA, outside_events, isect_tol);
3552  std::vector<double> line_t;
3553  if (inside_events.Count() != 2 && outside_events.Count() == 2) {
3554  // the direction of the curve should be opposite
3555  x[i].m_curve3d->Reverse();
3556  x[i].m_curveA->Reverse();
3557  x[i].m_curveB->Reverse();
3558  inside_lineA = outside_lineA;
3559  line_t.push_back(outside_events[0].m_a[0]);
3560  line_t.push_back(outside_events[1].m_a[0]);
3561  } else if (inside_events.Count() == 2 && outside_events.Count() != 2) {
3562  line_t.push_back(inside_events[0].m_a[0]);
3563  line_t.push_back(inside_events[1].m_a[0]);
3564  } else {
3565  bu_log("error: Cannot determine the correct direction of overlap event %d\n", i);
3566  continue;
3567  }
3568 
3569  // need to reverse inner loops to indicate overlap is outside
3570  // the closed region
3571  for (int j = original_count; j < x.Count(); j++) {
3572  // any curves that intersect the line crossing through the
3573  // loop may be inside the loop
3574  // FIXME: What about curves inside the loop that the line
3575  // happens to miss?
3576  if (i == j) {
3577  continue;
3578  }
3579  ON_SimpleArray<ON_X_EVENT> x_event;
3580  ON_Intersect(&inside_lineA, x[j].m_curveA, x_event, isect_tol);
3581  for (int k = 0; k < x_event.Count(); k++) {
3582  line_t.push_back(x_event[k].m_a[0]);
3583  }
3584  }
3585 
3586  // normalize all params in line_t
3587  for (size_t j = 0; j < line_t.size(); j++) {
3588  line_t[j] = inside_lineA.Domain().NormalizedParameterAt(line_t[j]);
3589  }
3590 
3591  std::sort(line_t.begin(), line_t.end());
3592  if (!ON_NearZero(line_t[0])) {
3593  bu_log("Error: param 0 is not in line_t!\n");
3594  continue;
3595  }
3596 
3597  // get a test point inside the region
3598  ON_3dPoint test_pt2d = inside_lineA.PointAt(line_t[1] * 0.5);
3599  ON_3dPoint test_pt = surfA->PointAt(test_pt2d.x, test_pt2d.y);
3600  ON_ClassArray<ON_PX_EVENT> px_event;
3601  if (!ON_Intersect(test_pt, *surfB, px_event, overlap_tol, 0, 0, treeB)) {
3602  // the test point is not overlapped
3603  x[i].m_curve3d->Reverse();
3604  x[i].m_curveA->Reverse();
3605  x[i].m_curveB->Reverse();
3606  overlap_events[i].m_type = Overlapevent::inner;
3607  } else {
3608  // the test point inside that region is an overlap point
3609  overlap_events[i].m_type = Overlapevent::outer;
3610  }
3611  }
3612 
3613  // find the inner events inside each outer event
3614  for (int i = 0; i < overlap_events.Count(); i++) {
3615  for (int j = 0; j < overlap_events.Count(); j++) {
3616  if (i == j || overlap_events[i].m_type != Overlapevent::outer || overlap_events[j].m_type != Overlapevent::inner) {
3617  continue;
3618  }
3619  if (overlap_events[i].IsCurveCompletelyIn(overlap_events[j].m_event->m_curve3d)) {
3620  overlap_events[i].m_inside.push_back(&(overlap_events[j]));
3621  }
3622  }
3623  }
3624 
3625  if (DEBUG_BREP_INTERSECT) {
3626  bu_log("%d overlap events.\n", overlap_events.Count());
3627  }
3628 
3629  if (surfA->IsPlanar() && surfB->IsPlanar() && overlap_events.Count()) {
3630  return x.Count() - original_count;
3631  }
3632 
3633  /* Second step: calculate the intersection of the bounding boxes.
3634  * Only the children of intersecting b-box pairs need to be considered.
3635  * The children will be generated only when they are needed, using the
3636  * method of splitting a NURBS surface.
3637  * So finally only a small subset of the surface tree is created.
3638  */
3639  typedef std::vector<std::pair<Subsurface *, Subsurface *> > NodePairs;
3640  NodePairs candidates, next_candidates;
3641  candidates.push_back(std::make_pair(rootA, rootB));
3642 
3643  for (int h = 0; h <= MAX_SSI_DEPTH && !candidates.empty(); h++) {
3644  next_candidates.clear();
3645  for (NodePairs::iterator i = candidates.begin(); i != candidates.end(); i++) {
3646  // If the Subsurface from either surfA or surfB (we're
3647  // checking surfA in this case) is completely inside an
3648  // overlap region, then we know that there is no
3649  // (non-overlap) intersection between this pair of
3650  // Subsurfaces, and thus no point in subdividing them
3651  // further.
3652  if (is_subsurfaceA_completely_inside_overlap(overlap_events, i->first, isect_tolA)) {
3653  continue;
3654  }
3655 
3656  // subdivide and add the parts whose bounding boxes
3657  // intersect to the new list of candidates
3658  std::vector<Subsurface *> partsA = subdivide_subsurface(i->first);
3659  std::vector<Subsurface *> partsB = subdivide_subsurface(i->second);
3660  for (size_t j = 0; j < partsA.size(); ++j) {
3661  for (size_t k = 0; k < partsB.size(); ++k) {
3662  if (partsB[k]->Intersect(*partsA[j], isect_tol)) {
3663  next_candidates.push_back(std::make_pair(partsA[j], partsB[k]));
3664  }
3665  }
3666  }
3667  }
3668  candidates = next_candidates;
3669  }
3670  if (DEBUG_BREP_INTERSECT) {
3671  bu_log("We get %lu intersection bounding boxes.\n", candidates.size());
3672  }
3673 
3674  /* Third step: get the intersection points using triangular approximation,
3675  * and then Newton iterations.
3676  */
3677  for (NodePairs::iterator i = candidates.begin(); i != candidates.end(); i++) {
3678  // We have arrived at the bottom of the trees.
3679  // Get an estimate of the intersection point lying on the intersection curve.
3680 
3681  // Get the corners of each surface sub-patch inside the bounding-box.
3682  ON_BoundingBox box_intersect;
3683  i->first->Intersect(*(i->second), isect_tol, &box_intersect);
3684 
3685  /* We approximate each surface sub-patch inside the bounding-box with two
3686  * triangles that share an edge.
3687  * The intersection of the surface sub-patches is approximated as the
3688  * intersection of triangles.
3689  */
3690  std::pair<Triangle, Triangle> trianglesA = get_subsurface_triangles(i->first, surfA);
3691  std::pair<Triangle, Triangle> trianglesB = get_subsurface_triangles(i->second, surfB);
3692 
3693  // calculate the intersection centers' uv coordinates
3694  int num_intersects = 0;
3695  ON_3dPoint average(0.0, 0.0, 0.0);
3696  ON_2dPoint avg_uvA(0.0, 0.0), avg_uvB(0.0, 0.0);
3697  for (int j = 0; j < 2; ++j) {
3698  const Triangle &triA = j ? trianglesA.second : trianglesA.first;
3699 
3700  for (int k = 0; k < 2; ++k) {
3701  const Triangle &triB = k ? trianglesB.second : trianglesB.first;
3702 
3703  ON_3dPoint intersection_center;
3704  if (triangle_intersection(triA, triB, intersection_center)) {
3705  ON_3dPoint bcA = triA.BarycentricCoordinate(intersection_center);
3706  ON_3dPoint bcB = triB.BarycentricCoordinate(intersection_center);
3707  ON_2dPoint uvA = barycentric_to_uv(bcA, triA);
3708  ON_2dPoint uvB = barycentric_to_uv(bcB, triB);
3709  average += intersection_center;
3710  avg_uvA += uvA;
3711  avg_uvB += uvB;
3712  ++num_intersects;
3713  }
3714  }
3715  }
3716 
3717  // The centroid of these intersection centers is the
3718  // intersection point we want.
3719  if (num_intersects) {
3720  average /= num_intersects;
3721  avg_uvA /= num_intersects;
3722  avg_uvB /= num_intersects;
3723  if (box_intersect.IsPointIn(average)) {
3724  // use Newton iterations to get an accurate intersection
3725  if (newton_ssi(avg_uvA.x, avg_uvA.y, avg_uvB.x, avg_uvB.y, surfA, surfB, isect_tol)) {
3726  average = surfA->PointAt(avg_uvA.x, avg_uvA.y);
3727  tmp_curvept.Append(average);
3728  tmp_curve_uvA.Append(avg_uvA);
3729  tmp_curve_uvB.Append(avg_uvB);
3730  }
3731  }
3732  }
3733  }
3734 
3735  add_points_to_closed_seams(surfA, surfB, tmp_curvept, tmp_curve_uvA, tmp_curve_uvB);
3736 
3737  // get just the unique, non-overlap intersection points
3738  for (int i = 0; i < tmp_curvept.Count(); i++) {
3739  bool unique_pt = true;
3740  for (int j = 0; j < curvept.Count(); j++) {
3741  if (tmp_curvept[i].DistanceTo(curvept[j]) < isect_tol &&
3742  tmp_curve_uvA[i].DistanceTo(curve_uvA[j]) < isect_tolA &&
3743  tmp_curve_uvB[i].DistanceTo(curve_uvB[j]) < isect_tolB)
3744  {
3745  unique_pt = false;
3746  break;
3747  }
3748  }
3749  if (unique_pt) {
3750  if (!is_uvA_completely_inside_overlap(overlap_events, tmp_curve_uvA[i])) {
3751  curvept.Append(tmp_curvept[i]);
3752  curve_uvA.Append(tmp_curve_uvA[i]);
3753  curve_uvB.Append(tmp_curve_uvB[i]);
3754  }
3755  }
3756  }
3757  if (DEBUG_BREP_INTERSECT) {
3758  bu_log("%d points on the intersection curves.\n", curvept.Count());
3759  }
3760 
3761  if (!curvept.Count()) {
3762  if (treeA == NULL) {
3763  delete rootA;
3764  }
3765  if (treeB == NULL) {
3766  delete rootB;
3767  }
3768 
3769  // should not return 0 as there might be overlap events
3770  return x.Count() - original_count;
3771  }
3772 
3773  // Fourth step: Fit the points in curvept into NURBS curves.
3774  // Here we use polyline approximation.
3775  // TODO: Find a better fitting algorithm unless this is good enough.
3776  // Time complexity: O(n^2*log(n)), really slow when n is large.
3777  // (n is the number of points generated above)
3778 
3779  // need to automatically generate a threshold
3780  // TODO: need more tests to find a better threshold
3781  double max_dist = std::min(rootA->m_surf->BoundingBox().Diagonal().Length(), rootB->m_surf->BoundingBox().Diagonal().Length()) * 0.1;
3782 
3783  double max_dist_uA = surfA->Domain(0).Length() * 0.05;
3784  double max_dist_vA = surfA->Domain(1).Length() * 0.05;
3785  double max_dist_uB = surfB->Domain(0).Length() * 0.05;
3786  double max_dist_vB = surfB->Domain(1).Length() * 0.05;
3787  if (DEBUG_BREP_INTERSECT) {
3788  bu_log("max_dist: %f\n", max_dist);
3789  bu_log("max_dist_uA: %f\n", max_dist_uA);
3790  bu_log("max_dist_vA: %f\n", max_dist_vA);
3791  bu_log("max_dist_uB: %f\n", max_dist_uB);
3792  bu_log("max_dist_vB: %f\n", max_dist_vB);
3793  }
3794 
3795  // identify the neighbors of each point
3796  std::vector<PointPair> ptpairs;
3797  for (int i = 0; i < curvept.Count(); i++) {
3798  for (int j = i + 1; j < curvept.Count(); j++) {
3799  PointPair ppair;
3800  ppair.distance3d = curvept[i].DistanceTo(curvept[j]);
3801  ppair.dist_uA = fabs(curve_uvA[i].x - curve_uvA[j].x);
3802  ppair.dist_vA = fabs(curve_uvA[i].y - curve_uvA[j].y);
3803  ppair.dist_uB = fabs(curve_uvB[i].x - curve_uvB[j].x);
3804  ppair.dist_vB = fabs(curve_uvB[i].y - curve_uvB[j].y);
3805  ppair.tol = isect_tol;
3806  if (ppair.dist_uA < max_dist_uA && ppair.dist_vA < max_dist_vA && ppair.dist_uB < max_dist_uB && ppair.dist_vB < max_dist_vB && ppair.distance3d < max_dist) {
3807  ppair.indexA = i;
3808  ppair.indexB = j;
3809  ptpairs.push_back(ppair);
3810  }
3811  }
3812  }
3813  std::sort(ptpairs.begin(), ptpairs.end());
3814 
3815  std::vector<ON_SimpleArray<int>*> polylines(curvept.Count());
3816 
3817  // polyline_of_terminal[i] = j means curvept[i] is a startpoint/endpoint of polylines[j]
3818  int *polyline_of_terminal = (int *)bu_malloc(curvept.Count() * sizeof(int), "int");
3819 
3820  // index of startpoints and endpoints of polylines[i]
3821  int *startpt = (int *)bu_malloc(curvept.Count() * sizeof(int), "int");
3822  int *endpt = (int *)bu_malloc(curvept.Count() * sizeof(int), "int");
3823 
3824  // initialize each polyline with only one point
3825  for (int i = 0; i < curvept.Count(); i++) {
3826  ON_SimpleArray<int> *single = new ON_SimpleArray<int>();
3827  single->Append(i);
3828  polylines[i] = single;
3829  polyline_of_terminal[i] = i;
3830  startpt[i] = i;
3831  endpt[i] = i;
3832  }
3833 
3834  // merge polylines with distance less than max_dist
3835  for (size_t i = 0; i < ptpairs.size(); i++) {
3836  int terminal1 = ptpairs[i].indexA;
3837  int terminal2 = ptpairs[i].indexB;
3838  int poly1 = polyline_of_terminal[terminal1];
3839  int poly2 = polyline_of_terminal[terminal2];
3840  if (poly1 == -1 || poly2 == -1) {
3841  continue;
3842  }
3843  // to merge poly2 into poly1, first start by making start/end points belong to poly1
3844  polyline_of_terminal[startpt[poly1]] = polyline_of_terminal[endpt[poly1]] = poly1;
3845  polyline_of_terminal[startpt[poly2]] = polyline_of_terminal[endpt[poly2]] = poly1;
3846 
3847  ON_SimpleArray<int> *line1 = polylines[poly1];
3848  ON_SimpleArray<int> *line2 = polylines[poly2];
3849  if (line1 != NULL && line2 != NULL && line1 != line2) {
3850  // Join the two polylines so that terminal1 and terminal2
3851  // are adjacent in the new unionline, reflecting the
3852  // spatial adjacency of the corresponding curve points.
3853  ON_SimpleArray<int> *unionline = new ON_SimpleArray<int>();
3854  bool start1 = (*line1)[0] == terminal1;
3855  bool start2 = (*line2)[0] == terminal2;
3856  if (start1) {
3857  if (start2) {
3858  // line1: [terminal1, ..., end1]
3859  // line2: [terminal2, ..., end2]
3860  // unionline: [end1, ..., terminal1, terminal2, ..., end2]
3861  line1->Reverse();
3862  unionline->Append(line1->Count(), line1->Array());
3863  unionline->Append(line2->Count(), line2->Array());
3864  startpt[poly1] = endpt[poly1];
3865  endpt[poly1] = endpt[poly2];
3866  } else {
3867  // line1: [terminal1, ..., end1]
3868  // line2: [start2, ..., terminal2]
3869  // unionline: [start2, ..., terminal2, terminal1, ..., end1]
3870  unionline->Append(line2->Count(), line2->Array());
3871  unionline->Append(line1->Count(), line1->Array());
3872  startpt[poly1] = startpt[poly2];
3873  endpt[poly1] = endpt[poly1];
3874  }
3875  } else {
3876  if (start2) {
3877  // line1: [start1, ..., terminal1]
3878  // line2: [terminal2, ..., end2]
3879  // unionline: [start1, ..., terminal1, terminal2, ..., end2]
3880  unionline->Append(line1->Count(), line1->Array());
3881  unionline->Append(line2->Count(), line2->Array());
3882  startpt[poly1] = startpt[poly1];
3883  endpt[poly1] = endpt[poly2];
3884  } else {
3885  // line1: [start1, ..., terminal1]
3886  // line2: [start2, ..., terminal2]
3887  // unionline: [start1, ..., terminal1, terminal2, ..., start2]
3888  unionline->Append(line1->Count(), line1->Array());
3889  line2->Reverse();
3890  unionline->Append(line2->Count(), line2->Array());
3891  startpt[poly1] = startpt[poly1];
3892  endpt[poly1] = startpt[poly2];
3893  }
3894  }
3895  polylines[poly1] = unionline;
3896  polylines[poly2] = NULL;
3897 
3898  // If lineN has more than one point, then joining to
3899  // terminalN has made it an interior point and it must be
3900  // invalidated as a terminal.
3901  if (line1->Count() >= 2) {
3902  polyline_of_terminal[terminal1] = -1;
3903  }
3904  if (line2->Count() >= 2) {
3905  polyline_of_terminal[terminal2] = -1;
3906  }
3907  delete line1;
3908  delete line2;
3909  }
3910  }
3911 
3912  // In some cases, the intersection curves will intersect with each
3913  // other. But with our merging mechanism, a point can only belong
3914  // to one curve. If curves need to share an intersection point, we
3915  // need some "seaming" segments to handle it.
3916  size_t num_curves = polylines.size();
3917  for (size_t i = 0; i < num_curves; i++) {
3918  if (!polylines[i]) {
3919  continue;
3920  }
3921  for (size_t j = i + 1; j < num_curves; j++) {
3922  if (!polylines[j]) {
3923  continue;
3924  }
3925  ON_SimpleArray<int> &polyline1 = *polylines[i];
3926  ON_SimpleArray<int> &polyline2 = *polylines[j];
3927 
3928  // find the closest pair of adjacent points between the two polylines
3929  int point_count1 = polyline1.Count();
3930  int point_count2 = polyline2.Count();
3931  PointPair pair;
3932  pair.distance3d = DBL_MAX;
3933  for (int k = 0; k < point_count1; k++) {
3934  for (int m = 0; m < point_count2; m++) {
3935  PointPair newpair;
3936  int start = polyline1[k], end = polyline2[m];
3937  newpair.distance3d = curvept[start].DistanceTo(curvept[end]);
3938  newpair.dist_uA = fabs(curve_uvA[start].x - curve_uvA[end].x);
3939  newpair.dist_vA = fabs(curve_uvA[start].y - curve_uvA[end].y);
3940  newpair.dist_uB = fabs(curve_uvB[start].x - curve_uvB[end].x);
3941  newpair.dist_vB = fabs(curve_uvB[start].y - curve_uvB[end].y);
3942  newpair.tol = isect_tol;
3943  if (newpair.dist_uA < max_dist_uA && newpair.dist_vA < max_dist_vA && newpair.dist_uB < max_dist_uB && newpair.dist_vB < max_dist_vB) {
3944  if (newpair < pair) {
3945  newpair.indexA = start;
3946  newpair.indexB = end;
3947  pair = newpair;
3948  }
3949  }
3950  }
3951  }
3952  if (pair.distance3d < max_dist) {
3953  // Generate a seaming curve if there are currently no
3954  // intersections.
3955  // TODO: These curve-curve intersections are
3956  // expensive. Is this really necessary?
3957  ON_3dPointArray uvA1, uvA2, uvB1, uvB2;
3958  for (int k = 0; k < polyline1.Count(); k++) {
3959  uvA1.Append(curve_uvA[polyline1[k]]);
3960  uvB1.Append(curve_uvB[polyline1[k]]);
3961  }
3962  for (int k = 0; k < polyline2.Count(); k++) {
3963  uvA2.Append(curve_uvA[polyline2[k]]);
3964  uvB2.Append(curve_uvB[polyline2[k]]);
3965  }
3966  ON_PolylineCurve curveA1(uvA1), curveA2(uvA2), curveB1(uvB1), curveB2(uvB2);
3967  ON_SimpleArray<ON_X_EVENT> x_event1, x_event2;
3968  if (ON_Intersect(&curveA1, &curveA2, x_event1, isect_tol) &&
3969  ON_Intersect(&curveB1, &curveB2, x_event2, isect_tol)) {
3970  continue;
3971  }
3972 
3973  // If the seaming curve is continuous to polylines[i]
3974  // or polylines[j], we can just merge the curve with
3975  // the polyline, otherwise we generate a new segment.
3976  if (!join_continuous_ppair_to_polyline(polylines[i], pair) &&
3977  !join_continuous_ppair_to_polyline(polylines[j], pair))
3978  {
3979  ON_SimpleArray<int> *seam = new ON_SimpleArray<int>;
3980  seam->Append(pair.indexA);
3981  seam->Append(pair.indexB);
3982  polylines.push_back(seam);
3983  }
3984  }
3985  }
3986  }
3987 
3988  // generate ON_Curves from the polylines
3989  ON_SimpleArray<ON_Curve *> intersect3d, intersect_uvA, intersect_uvB;
3990  ON_SimpleArray<int> single_pts;
3991  for (size_t i = 0; i < polylines.size(); i++) {
3992  if (polylines[i] == NULL) {
3993  continue;
3994  }
3995  bool is_seam = (i >= num_curves);
3996 
3997  ON_SimpleArray<int> &polyline = *polylines[i];
3998  int startpoint = polyline[0];
3999  int endpoint = polyline[polyline.Count() - 1];
4000 
4001  if (polyline.Count() == 1) {
4002  single_pts.Append(startpoint);
4003  delete polylines[i];
4004  continue;
4005  }
4006 
4007  // curve in 3D space
4008  ON_3dPointArray ptarray;
4009  for (int j = 0; j < polyline.Count(); j++) {
4010  ptarray.Append(curvept[polyline[j]]);
4011  }
4012  // close curve if it forms a loop
4013  if (!is_seam && curvept[startpoint].DistanceTo(curvept[endpoint]) < max_dist) {
4014  ptarray.Append(curvept[startpoint]);
4015  }
4016  ON_PolylineCurve *curve = new ON_PolylineCurve(ptarray);
4017  intersect3d.Append(curve);
4018 
4019  // curve in UV space (surfA)
4020  ptarray.Empty();
4021  for (int j = 0; j < polyline.Count(); j++) {
4022  ON_2dPoint &pt2d = curve_uvA[polyline[j]];
4023  ptarray.Append(ON_3dPoint(pt2d.x, pt2d.y, 0.0));
4024  }
4025  // close curve if it forms a loop (happens rarely compared to 3D)
4026  if (!is_seam &&
4027  fabs(curve_uvA[startpoint].x - curve_uvA[endpoint].x) < max_dist_uA &&
4028  fabs(curve_uvA[startpoint].y - curve_uvA[endpoint].y) < max_dist_vA)
4029  {
4030  ON_2dPoint &pt2d = curve_uvA[startpoint];
4031  ptarray.Append(ON_3dPoint(pt2d.x, pt2d.y, 0.0));
4032  }
4033  curve = new ON_PolylineCurve(ptarray);
4034  curve->ChangeDimension(2);
4035  intersect_uvA.Append(curve_fitting(curve, fitting_tolA));
4036 
4037  // curve in UV space (surfB)
4038  ptarray.Empty();
4039  for (int j = 0; j < polyline.Count(); j++) {
4040  ON_2dPoint &pt2d = curve_uvB[polyline[j]];
4041  ptarray.Append(ON_3dPoint(pt2d.x, pt2d.y, 0.0));
4042  }
4043  // close curve if it forms a loop (happens rarely compared to 3D)
4044  if (!is_seam &&
4045  fabs(curve_uvB[startpoint].x - curve_uvB[endpoint].x) < max_dist_uB &&
4046  fabs(curve_uvB[startpoint].y - curve_uvB[endpoint].y) < max_dist_vB)
4047  {
4048  ON_2dPoint &pt2d = curve_uvB[startpoint];
4049  ptarray.Append(ON_3dPoint(pt2d.x, pt2d.y, 0.0));
4050  }
4051  curve = new ON_PolylineCurve(ptarray);
4052  curve->ChangeDimension(2);
4053  intersect_uvB.Append(curve_fitting(curve, fitting_tolB));
4054 
4055  delete polylines[i];
4056  }
4057 
4058  if (DEBUG_BREP_INTERSECT) {
4059  bu_log("%d curve segments and %d single points.\n", intersect3d.Count(), single_pts.Count());
4060  }
4061  bu_free(polyline_of_terminal, "int");
4062  bu_free(startpt, "int");
4063  bu_free(endpt, "int");
4064 
4065  // generate ON_SSX_EVENTs
4066  for (int i = 0; i < intersect3d.Count(); i++) {
4067  ON_SSX_EVENT event;
4068  event.m_curve3d = intersect3d[i];
4069  event.m_curveA = intersect_uvA[i];
4070  event.m_curveB = intersect_uvB[i];
4071  // Normalize the curves, so that their domains are the same,
4072  // which is required by ON_SSX_EVENT::IsValid().
4073  event.m_curve3d->SetDomain(ON_Interval(0.0, 1.0));
4074  event.m_curveA->SetDomain(ON_Interval(0.0, 1.0));
4075  event.m_curveB->SetDomain(ON_Interval(0.0, 1.0));
4076  // If the surfA and surfB normals of all points are
4077  // parallel, the intersection is considered tangent.
4078  bool is_tangent = true;
4079  int count = std::min(event.m_curveA->SpanCount(), event.m_curveB->SpanCount());
4080  for (int j = 0; j <= count; ++j) {
4081  ON_3dVector normalA, normalB;
4082  ON_3dPoint pointA = event.m_curveA->PointAt((double)j / count);
4083  ON_3dPoint pointB = event.m_curveB->PointAt((double)j / count);
4084  if (!(surfA->EvNormal(pointA.x, pointA.y, normalA) &&
4085  surfB->EvNormal(pointB.x, pointB.y, normalB) &&
4086  normalA.IsParallelTo(normalB)))
4087  {
4088  is_tangent = false;
4089  break;
4090  }
4091  }
4092  if (is_tangent) {
4093  // For ssx_tangent events, the 3d curve direction may
4094  // not agree with SurfaceNormalA x SurfaceNormalB
4095  // (See opennurbs/opennurbs_x.h).
4096  ON_3dVector direction = event.m_curve3d->TangentAt(0);
4097  ON_3dVector SurfaceNormalA = surfA->NormalAt(
4098  event.m_curveA->PointAtStart().x,
4099  event.m_curveA->PointAtStart().y);
4100  ON_3dVector SurfaceNormalB = surfB->NormalAt(
4101  event.m_curveB->PointAtStart().x,
4102  event.m_curveB->PointAtStart().y);
4103  if (ON_DotProduct(direction, ON_CrossProduct(SurfaceNormalB, SurfaceNormalA)) < 0) {
4104  if (!(event.m_curve3d->Reverse() &&
4105  event.m_curveA->Reverse() &&
4106  event.m_curveB->Reverse()))
4107  {
4108  bu_log("warning: reverse failed. The direction of %d might be wrong.\n",
4109  x.Count() - original_count);
4110  }
4111  }
4112  event.m_type = ON_SSX_EVENT::ssx_tangent;
4113  } else {
4114  event.m_type = ON_SSX_EVENT::ssx_transverse;
4115  }
4116  // ssx_overlap is handled above
4117  x.Append(event);
4118  // set the curves to NULL in case they are deleted by ~ON_SSX_EVENT()
4119  event.m_curve3d = event.m_curveA = event.m_curveB = NULL;
4120  }
4121 
4122  for (int i = 0; i < single_pts.Count(); i++) {
4123  bool unique_pt = true;
4124  for (int j = 0; j < intersect3d.Count(); ++j) {
4125  ON_ClassArray<ON_PX_EVENT> px_event;
4126  if (ON_Intersect(curvept[single_pts[i]], *intersect3d[j], px_event)) {
4127  unique_pt = false;
4128  break;
4129  }
4130  }
4131  if (unique_pt) {
4132  ON_SSX_EVENT event;
4133  event.m_point3d = curvept[single_pts[i]];
4134  event.m_pointA = curve_uvA[single_pts[i]];
4135  event.m_pointB = curve_uvB[single_pts[i]];
4136  // If the surfA and surfB normals are parallel, the
4137  // intersection is considered tangent.
4138  ON_3dVector normalA, normalB;
4139  if (surfA->EvNormal(event.m_pointA.x, event.m_pointA.y, normalA) &&
4140  surfB->EvNormal(event.m_pointB.x, event.m_pointB.y, normalB) &&
4141  normalA.IsParallelTo(normalB))
4142  {
4143  event.m_type = ON_SSX_EVENT::ssx_tangent_point;
4144  } else {
4145  event.m_type = ON_SSX_EVENT::ssx_transverse_point;
4146  }
4147  x.Append(event);
4148  }
4149  }
4150 
4151  if (treeA == NULL) {
4152  delete rootA;
4153  }
4154  if (treeB == NULL) {
4155  delete rootB;
4156  }
4157  return x.Count() - original_count;
4158 }
4159 
4160 
4161 // Local Variables:
4162 // tab-width: 8
4163 // mode: C++
4164 // c-basic-offset: 4
4165 // indent-tabs-mode: t
4166 // c-file-style: "stroustrup"
4167 // End:
4168 // ex: shiftwidth=4 tabstop=8
void SetAOverlapRange(const ON_Interval &interval)
Definition: intersect.cpp:139
ON_3dPoint c
Definition: intersect.cpp:2159
bool ON_Intersect(const ON_3dPoint &pointA, const ON_3dPoint &pointB, ON_ClassArray< ON_PX_EVENT > &x, double tol)
Definition: intersect.cpp:647
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define CCI_OVERLAP_TEST_POINTS
Definition: intersect.cpp:993
bool IsPointIn(const ON_2dPoint &_pt) const
Definition: intersect.cpp:2551
#define OVERLAPS_LINKED(from, to)
ON_2dPoint a_2d
Definition: intersect.cpp:2160
HIDDEN ON_2dPoint barycentric_to_uv(const ON_3dPoint &bc, const Triangle &tri)
Definition: intersect.cpp:3074
XEventProxy(ON_X_EVENT::TYPE type)
Definition: intersect.cpp:102
Overlapevent(ON_SSX_EVENT *_e)
Definition: intersect.cpp:2534
void SetBPoints(const ON_3dPoint &start, const ON_3dPoint &end)
Definition: intersect.cpp:131
Definition: clone.c:90
void SetBSurfaceParameters(double u1, double v1, double u2, double v2)
Definition: intersect.cpp:197
HIDDEN bool newton_ssi(double &uA, double &vA, double &uB, double &vB, const ON_Surface *surfA, const ON_Surface *surfB, double isect_tol)
Definition: intersect.cpp:2390
HIDDEN void newton_cci(double &t_a, double &t_b, const ON_Curve *curveA, const ON_Curve *curveB, double isect_tol)
Definition: intersect.cpp:996
Definition: nmg_tri.c:75
HIDDEN int isocurve_surface_overlap_location(const ON_Surface *surf1, double iso_t, int iso_dir, double overlap_start, double overlap_end, const ON_Surface *surf2, Subsurface *surf2_tree)
Definition: intersect.cpp:2741
ON_3dPoint BarycentricCoordinate(const ON_3dPoint &pt) const
Definition: intersect.cpp:2167
OverlapSegment(void)
Definition: intersect.cpp:2487
Definition: raytrace.h:368
void SetCurvesToNull()
Definition: intersect.cpp:2497
bool IsBoxCompletelyIn(const ON_2dPoint &_min, const ON_2dPoint &_max) const
Definition: intersect.cpp:2611
bool IsBoxCompletelyOut(const ON_2dPoint &_min, const ON_2dPoint &_max) const
Definition: intersect.cpp:2622
bool IsPointOut(const ON_2dPoint &_pt) const
Definition: intersect.cpp:2577
ON_2dPoint c_2d
Definition: intersect.cpp:2160
HIDDEN ON_Curve * link_curves(ON_Curve *&c1, ON_Curve *&c2)
Definition: intersect.cpp:2456
HIDDEN bool build_curve_root(const ON_Curve *curve, const ON_Interval *domain, Subcurve &root)
Definition: intersect.cpp:350
ON_3dPoint a
Definition: intersect.cpp:2159
Header file for the BRL-CAD common definitions.
ON_X_EVENT Event(void)
Definition: intersect.cpp:214
HIDDEN bool newton_pci(double &t, const ON_3dPoint &pointA, const ON_Curve &curveB, double tol)
Definition: intersect.cpp:686
#define PCI_MAX_ITERATIONS
Definition: intersect.cpp:67
#define HIDDEN
Definition: common.h:86
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
#define CCI_DEFAULT_TOLERANCE
Definition: intersect.cpp:59
#define SSI_DEFAULT_TOLERANCE
Definition: intersect.cpp:61
double dist_vA
Definition: intersect.cpp:2333
void GetLineSegments(ON_Line line[3]) const
Definition: intersect.cpp:2181
bool operator<(const PointPair &pp) const
Definition: intersect.cpp:2335
HIDDEN void split_overlaps_at_intersections(ON_SimpleArray< OverlapSegment * > &overlaps, const ON_Surface *surfA, const ON_Surface *surfB, Subsurface *treeA, Subsurface *treeB, double isect_tol, double isect_tolA, double isect_tolB)
Definition: intersect.cpp:2912
HIDDEN bool is_pt_in_surf_overlap(ON_2dPoint surf1_pt, const ON_Surface *surf1, const ON_Surface *surf2, Subsurface *surf2_tree)
Definition: intersect.cpp:2640
void SetBCurveParameter(double t)
Definition: intersect.cpp:190
ON_3dPoint b
Definition: intersect.cpp:2159
HIDDEN std::pair< Triangle, Triangle > get_subsurface_triangles(const Subsurface *sub, const ON_Surface *surf)
Definition: intersect.cpp:3037
bool IsCurveCompletelyIn(const ON_Curve *_curve) const
Definition: intersect.cpp:2631
#define MAX_PSI_DEPTH
Definition: intersect.cpp:47
#define CSI_DEFAULT_TOLERANCE
Definition: intersect.cpp:60
HIDDEN bool triangle_intersection(const struct Triangle &triA, const struct Triangle &triB, ON_3dPoint &center)
Definition: intersect.cpp:2191
HIDDEN ON_Curve * curve_fitting(ON_Curve *in, double fitting_tol=ON_ZERO_TOLERANCE, bool delete_curve=true)
Definition: intersect.cpp:477
enum Overlapevent::TYPE m_type
double distance3d
Definition: intersect.cpp:2333
ON_Curve * m_curveA
Definition: intersect.cpp:2478
double dist_uA
Definition: intersect.cpp:2333
HIDDEN bool build_surface_root(const ON_Surface *surf, const ON_Interval *u_domain, const ON_Interval *v_domain, Subsurface &root)
Definition: intersect.cpp:399
void SetACurveParameter(double t)
Definition: intersect.cpp:161
HIDDEN void newton_csi(double &t, double &u, double &v, const ON_Curve *curveA, const ON_Surface *surfB, double isect_tol, Subsurface *tree)
Definition: intersect.cpp:1518
goto out
Definition: nmg_mod.c:3846
HIDDEN Subsurface * get_surface_root(const ON_Surface *surf, const ON_Interval *u_domain, const ON_Interval *v_domain)
Definition: intersect.cpp:459
HIDDEN bool is_uvA_completely_inside_overlap(const ON_SimpleArray< Overlapevent > &overlap_events, const ON_2dPoint &pt)
Definition: intersect.cpp:2704
bool IsBoxIntersected(const ON_2dPoint &_min, const ON_2dPoint &_max) const
Definition: intersect.cpp:2588
void CreateFromPoints(ON_3dPoint &pA, ON_3dPoint &pB, ON_3dPoint &pC)
Definition: intersect.cpp:2161
ON_Surface * sub_surface(const ON_Surface *in, int dir, double a, double b)
Definition: intersect.cpp:286
#define MAX_CCI_DEPTH
Definition: intersect.cpp:48
void SetAPoints(const ON_3dPoint &start, const ON_3dPoint &end)
Definition: intersect.cpp:116
#define PSI_MAX_ITERATIONS
Definition: intersect.cpp:68
ON_Curve * m_curveB
Definition: intersect.cpp:2478
bool IsPointOnBoundary(const ON_2dPoint &_pt) const
Definition: intersect.cpp:2540
#define ZERO(val)
Definition: units.c:38
int overlap(struct application *ap, struct partition *pp, struct region *reg1, struct region *reg2, struct partition *hp)
Definition: gqa.c:843
ON_Curve * sub_curve(const ON_Curve *in, double a, double b)
Definition: intersect.cpp:220
HIDDEN Subcurve * get_curve_root(const ON_Curve *curve, const ON_Interval *domain)
Definition: intersect.cpp:382
double tol
Definition: intersect.cpp:2334
HIDDEN bool join_continuous_ppair_to_polyline(ON_SimpleArray< int > *polyline, const PointPair &ppair)
Definition: intersect.cpp:3084
#define MAX_CSI_DEPTH
Definition: intersect.cpp:49
void SetBOverlapRange(const ON_Interval &interval)
Definition: intersect.cpp:168
#define CSI_MAX_ITERATIONS
Definition: intersect.cpp:70
HIDDEN std::vector< Subsurface * > subdivide_subsurface(Subsurface *sub)
Definition: intersect.cpp:3022
#define MAX_SSI_DEPTH
Definition: intersect.cpp:50
ON_ClassArray< ON_3dPointArray > get_overlap_intersection_parameters(const ON_SimpleArray< OverlapSegment * > &overlaps, const ON_Surface *surfA, const ON_Surface *surfB, Subsurface *treeA, Subsurface *treeB, double isect_tol, double isect_tolA, double isect_tolB)
Definition: intersect.cpp:2782
ON_2dPoint Get2DParam(double t)
Definition: intersect.cpp:2492
HIDDEN bool newton_psi(double &u, double &v, const ON_3dPoint &pointA, const ON_Surface &surfB, double tol)
Definition: intersect.cpp:835
HIDDEN void add_points_to_closed_seams(const ON_Surface *surfA, const ON_Surface *surfB, ON_3dPointArray &curvept, ON_2dPointArray &curve_uvA, ON_2dPointArray &curve_uvB)
Definition: intersect.cpp:2346
double dist_uB
Definition: intersect.cpp:2333
#define DEBUG_BREP_INTERSECT
Definition: intersect.cpp:41
void SetAPoint(const ON_3dPoint &pt)
Definition: intersect.cpp:109
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
ON_2dPoint b_2d
Definition: intersect.cpp:2160
#define PCI_DEFAULT_TOLERANCE
Definition: intersect.cpp:57
#define CSI_OVERLAP_TEST_POINTS
Definition: intersect.cpp:1515
#define MAX_PCI_DEPTH
Definition: intersect.cpp:46
std::vector< Overlapevent * > m_inside
Definition: intersect.cpp:2530
#define CCI_MAX_ITERATIONS
Definition: intersect.cpp:69
void SetBPoint(const ON_3dPoint &pt)
Definition: intersect.cpp:124
ON_BoundingBox m_bboxA
Definition: intersect.cpp:2523
HIDDEN const ON_Interval check_domain(const ON_Interval *in, const ON_Interval &domain, const char *name)
Definition: intersect.cpp:325
bool is_closed(const ON_Surface *surf)
HIDDEN const point_t delta
Definition: sh_prj.c:618
#define SSI_MAX_ITERATIONS
Definition: intersect.cpp:71
ON_SSX_EVENT * m_event
Definition: intersect.cpp:2521
ON_Curve * m_curve3d
Definition: intersect.cpp:2478
double tolerance_2d_from_3d(double tol_3d, Subcurve *root, const ON_Curve *curve, const ON_Interval *curve_domain)
Definition: intersect.cpp:1078
HIDDEN bool is_subsurfaceA_completely_inside_overlap(const ON_SimpleArray< Overlapevent > &overlap_events, const Subsurface *subA, double isect_tolA)
Definition: intersect.cpp:2666
double dist_vB
Definition: intersect.cpp:2333
HIDDEN bool is_valid_overlap(const OverlapSegment *overlap)
Definition: intersect.cpp:2899
bool ON_NearZero(double val, double epsilon)
Return truthfully whether a value is within a specified epsilon distance from zero.
void SetBSurfaceParameter(double u, double v)
Definition: intersect.cpp:207