BRL-CAD
PullbackCurve.cpp
Go to the documentation of this file.
1 /* P U L L B A C K C U R V E . C P P
2  * BRL-CAD
3  *
4  * Copyright (c) 2009-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 step/PullbackCurve.cpp
21  *
22  * Pull curve back into UV space from 3D space
23  *
24  */
25 
26 #include "common.h"
27 
28 #include "vmath.h"
29 #include "dvec.h"
30 
31 #include <assert.h>
32 #include <vector>
33 #include <list>
34 #include <limits>
35 #include <set>
36 #include <map>
37 #include <string>
38 
39 #include "brep.h"
40 #include "bu/parallel.h"
41 
42 /* interface header */
43 #include "PullbackCurve.h"
44 
45 #define RANGE_HI 0.55
46 #define RANGE_LO 0.45
47 #define UNIVERSAL_SAMPLE_COUNT 1001
48 
49 #define _DEBUG_TESTING_
50 #ifdef _DEBUG_TESTING_
51 bool _debug_print_ = false;
55 #endif
56 
57 /* FIXME: duplicated with opennurbs_ext.cpp */
58 class BSpline
59 {
60 public:
61  int p; // degree
62  int m; // num_knots-1
63  int n; // num_samples-1 (aka number of control points)
64  std::vector<double> params;
65  std::vector<double> knots;
66  ON_2dPointArray controls;
67 };
68 
69 bool
70 isFlat(const ON_2dPoint& p1, const ON_2dPoint& m, const ON_2dPoint& p2, double flatness)
71 {
72  ON_Line line = ON_Line(ON_3dPoint(p1), ON_3dPoint(p2));
73  return line.DistanceTo(ON_3dPoint(m)) <= flatness;
74 }
75 
76 void
77 utah_ray_planes(const ON_Ray &r, ON_3dVector &p1, double &p1d, ON_3dVector &p2, double &p2d)
78 {
79  ON_3dPoint ro(r.m_origin);
80  ON_3dVector rdir(r.m_dir);
81  double rdx, rdy, rdz;
82  double rdxmag, rdymag, rdzmag;
83 
84  rdx = rdir.x;
85  rdy = rdir.y;
86  rdz = rdir.z;
87 
88  rdxmag = fabs(rdx);
89  rdymag = fabs(rdy);
90  rdzmag = fabs(rdz);
91 
92  if (rdxmag > rdymag && rdxmag > rdzmag)
93  p1 = ON_3dVector(rdy, -rdx, 0);
94  else
95  p1 = ON_3dVector(0, rdz, -rdy);
96  p1.Unitize();
97 
98  p2 = ON_CrossProduct(p1, rdir);
99 
100  p1d = -(p1 * ro);
101  p2d = -(p2 * ro);
102 }
103 
104 
105 enum seam_direction
106 seam_direction(ON_2dPoint uv1, ON_2dPoint uv2)
107 {
108  if (NEAR_EQUAL(uv1.x, 0.0, PBC_TOL) && NEAR_EQUAL(uv2.x, 0.0, PBC_TOL)) {
109  return WEST_SEAM;
110  } else if (NEAR_EQUAL(uv1.x, 1.0, PBC_TOL) && NEAR_EQUAL(uv2.x, 1.0, PBC_TOL)) {
111  return EAST_SEAM;
112  } else if (NEAR_EQUAL(uv1.y, 0.0, PBC_TOL) && NEAR_EQUAL(uv2.y, 0.0, PBC_TOL)) {
113  return SOUTH_SEAM;
114  } else if (NEAR_EQUAL(uv1.y, 1.0, PBC_TOL) && NEAR_EQUAL(uv2.y, 1.0, PBC_TOL)) {
115  return NORTH_SEAM;
116  } else {
117  return UNKNOWN_SEAM_DIRECTION;
118  }
119 }
120 
121 
122 ON_BOOL32
124  const ON_Surface *surf,
125  const ON_Interval &u_interval,
126  const ON_Interval &v_interval,
127  ON_Interval domSplits[2][2]
128  )
129 {
130  ON_3dPoint min, max;
131  ON_Interval dom[2];
132  double length[2];
133 
134  // initialize intervals
135  for (int i = 0; i < 2; i++) {
136  for (int j = 0; j < 2; j++) {
137  domSplits[i][j] = ON_Interval::EmptyInterval;
138  }
139  }
140 
141  min[0] = u_interval.m_t[0];
142  min[1] = v_interval.m_t[0];
143  max[0] = u_interval.m_t[1];
144  max[1] = v_interval.m_t[1];
145 
146  for (int i=0; i<2; i++) {
147  if (surf->IsClosed(i)) {
148  dom[i] = surf->Domain(i);
149  length[i] = dom[i].Length();
150  if ((max[i] - min[i]) >= length[i]) {
151  domSplits[i][0] = dom[i];
152  } else {
153  double minbound = min[i];
154  double maxbound = max[i];
155  while (minbound < dom[i].m_t[0]) {
156  minbound += length[i];
157  maxbound += length[i];
158  }
159  while (minbound > dom[i].m_t[1]) {
160  minbound -= length[i];
161  maxbound -= length[i];
162  }
163  if (maxbound > dom[i].m_t[1]) {
164  domSplits[i][0].m_t[0] = minbound;
165  domSplits[i][0].m_t[1] = dom[i].m_t[1];
166  domSplits[i][1].m_t[0] = dom[i].m_t[0];
167  domSplits[i][1].m_t[1] = maxbound - length[i];
168  } else {
169  domSplits[i][0].m_t[0] = minbound;
170  domSplits[i][0].m_t[1] = maxbound;
171  }
172  }
173  } else {
174  domSplits[i][0].m_t[0] = min[i];
175  domSplits[i][0].m_t[1] = max[i];
176  }
177  }
178 
179  return true;
180 }
181 
182 
183 /*
184  * Wrapped OpenNURBS 'EvNormal()' because it fails when at surface singularity
185  * but not on edge of domain. If fails and at singularity this wrapper will
186  * reevaluate at domain edge.
187  */
188 ON_BOOL32
189 surface_EvNormal( // returns false if unable to evaluate
190  const ON_Surface *surf,
191  double s, double t, // evaluation parameters (s,t)
192  ON_3dPoint& point, // returns value of surface
193  ON_3dVector& normal, // unit normal
194  int side, // optional - determines which side to evaluate from
195  // 0 = default
196  // 1 from NE quadrant
197  // 2 from NW quadrant
198  // 3 from SW quadrant
199  // 4 from SE quadrant
200  int* hint // optional - evaluation hint (int[2]) used to speed
201  // repeated evaluations
202  )
203 {
204  ON_BOOL32 rc = false;
205 
206  if (!(rc=surf->EvNormal(s, t, point, normal, side, hint))) {
207  side = IsAtSingularity(surf, s, t, PBC_SEAM_TOL);// 0 = south, 1 = east, 2 = north, 3 = west
208  if (side >= 0) {
209  ON_Interval u = surf->Domain(0);
210  ON_Interval v = surf->Domain(1);
211  if (side == 0) {
212  rc=surf->EvNormal(u.m_t[0], v.m_t[0], point, normal, side, hint);
213  } else if (side == 1) {
214  rc=surf->EvNormal(u.m_t[1], v.m_t[1], point, normal, side, hint);
215  } else if (side == 2) {
216  rc=surf->EvNormal(u.m_t[1], v.m_t[1], point, normal, side, hint);
217  } else if (side == 3) {
218  rc=surf->EvNormal(u.m_t[0], v.m_t[0], point, normal, side, hint);
219  }
220  } else {
221  /*
222  * brute force and try to solve from each side of the surface domain
223  */
224  ON_Interval u = surf->Domain(0);
225  ON_Interval v = surf->Domain(1);
226  for(int iside=1; iside <= 4; iside++) {
227  rc=surf->EvNormal(s, t, point, normal, iside, hint);
228  if (rc)
229  break;
230  }
231  }
232  }
233  return rc;
234 }
235 
236 
237 static bool locals_initialized[MAX_PSW] = {false};
238 static ON_RevSurface *rev_surface[MAX_PSW] = {NULL};
239 static ON_NurbsSurface *nurbs_surface[MAX_PSW] = {NULL};
240 static ON_Extrusion *extr_surface[MAX_PSW] = {NULL};
241 static ON_PlaneSurface *plane_surface[MAX_PSW] = {NULL};
242 static ON_SumSurface *sum_surface[MAX_PSW] = {NULL};
243 static ON_SurfaceProxy *proxy_surface[MAX_PSW] = {NULL};
244 
245 ON_BOOL32
247  const ON_Surface *surf,
248  const ON_Interval &u_interval,
249  const ON_Interval &v_interval,
250  ON_BoundingBox& bbox,
251  ON_BOOL32 bGrowBox
252  )
253 {
254  int p = bu_parallel_id();
255 
256  if (!locals_initialized[p]) {
258  rev_surface[p] = ON_RevSurface::New();
259  nurbs_surface[p] = ON_NurbsSurface::New();
260  extr_surface[p] = new ON_Extrusion();
261  plane_surface[p] = new ON_PlaneSurface();
262  sum_surface[p] = ON_SumSurface::New();
263  proxy_surface[p] = new ON_SurfaceProxy();
264  locals_initialized[p] = true;
266  }
267 
268  ON_Interval domSplits[2][2] = { { ON_Interval::EmptyInterval, ON_Interval::EmptyInterval }, { ON_Interval::EmptyInterval, ON_Interval::EmptyInterval }};
269  if (!GetDomainSplits(surf,u_interval,v_interval,domSplits)) {
270  return false;
271  }
272 
273  bool growcurrent = bGrowBox;
274  for (int i=0; i<2; i++) {
275  if (domSplits[0][i] == ON_Interval::EmptyInterval)
276  continue;
277 
278  for (int j=0; j<2; j++) {
279  if (domSplits[1][j] != ON_Interval::EmptyInterval) {
280  if (dynamic_cast<ON_RevSurface * >(const_cast<ON_Surface *>(surf)) != NULL) {
281  *rev_surface[p] = *dynamic_cast<ON_RevSurface * >(const_cast<ON_Surface *>(surf));
282  if (rev_surface[p]->Trim(0, domSplits[0][i]) && rev_surface[p]->Trim(1, domSplits[1][j])) {
283  if (!rev_surface[p]->GetBoundingBox(bbox, growcurrent)) {
284  return false;
285  }
286  growcurrent = true;
287  }
288  } else if (dynamic_cast<ON_NurbsSurface * >(const_cast<ON_Surface *>(surf)) != NULL) {
289  *nurbs_surface[p] = *dynamic_cast<ON_NurbsSurface * >(const_cast<ON_Surface *>(surf));
290  if (nurbs_surface[p]->Trim(0, domSplits[0][i]) && nurbs_surface[p]->Trim(1, domSplits[1][j])) {
291  if (!nurbs_surface[p]->GetBoundingBox(bbox, growcurrent)) {
292  return false;
293  }
294  }
295  growcurrent = true;
296  } else if (dynamic_cast<ON_Extrusion * >(const_cast<ON_Surface *>(surf)) != NULL) {
297  *extr_surface[p] = *dynamic_cast<ON_Extrusion * >(const_cast<ON_Surface *>(surf));
298  if (extr_surface[p]->Trim(0, domSplits[0][i]) && extr_surface[p]->Trim(1, domSplits[1][j])) {
299  if (!extr_surface[p]->GetBoundingBox(bbox, growcurrent)) {
300  return false;
301  }
302  }
303  growcurrent = true;
304  } else if (dynamic_cast<ON_PlaneSurface * >(const_cast<ON_Surface *>(surf)) != NULL) {
305  *(plane_surface[p]) = *dynamic_cast<ON_PlaneSurface * >(const_cast<ON_Surface *>(surf));
306  if (plane_surface[p]->Trim(0, domSplits[0][i]) && plane_surface[p]->Trim(1, domSplits[1][j])) {
307  if (!plane_surface[p]->GetBoundingBox(bbox, growcurrent)) {
308  return false;
309  }
310  }
311  growcurrent = true;
312  } else if (dynamic_cast<ON_SumSurface * >(const_cast<ON_Surface *>(surf)) != NULL) {
313  *sum_surface[p] = *dynamic_cast<ON_SumSurface * >(const_cast<ON_Surface *>(surf));
314  if (sum_surface[p]->Trim(0, domSplits[0][i]) && sum_surface[p]->Trim(1, domSplits[1][j])) {
315  if (!sum_surface[p]->GetBoundingBox(bbox, growcurrent)) {
316  return false;
317  }
318  }
319  growcurrent = true;
320  } else if (dynamic_cast<ON_SurfaceProxy * >(const_cast<ON_Surface *>(surf)) != NULL) {
321  *proxy_surface[p] = *dynamic_cast<ON_SurfaceProxy * >(const_cast<ON_Surface *>(surf));
322  if (proxy_surface[p]->Trim(0, domSplits[0][i]) && proxy_surface[p]->Trim(1, domSplits[1][j])) {
323  if (!proxy_surface[p]->GetBoundingBox(bbox, growcurrent)) {
324  return false;
325  }
326  }
327  growcurrent = true;
328  } else {
329  std::cerr << "Unknown Surface Type" << std::endl;
330  }
331  }
332  }
333  }
334 
335  return true;
336 }
337 
338 
339 ON_BOOL32
341  const ON_BrepFace &face,
342  ON_BoundingBox& bbox,
343  ON_BOOL32 bGrowBox
344  )
345 {
346  const ON_Surface *surf = face.SurfaceOf();
347 
348  // may be a smaller trimmed subset of surface so worth getting
349  // face boundary
350  bool growcurrent = bGrowBox;
351  ON_3dPoint min, max;
352  for (int li = 0; li < face.LoopCount(); li++) {
353  for (int ti = 0; ti < face.Loop(li)->TrimCount(); ti++) {
354  ON_BrepTrim *trim = face.Loop(li)->Trim(ti);
355  trim->GetBoundingBox(min, max, growcurrent);
356  growcurrent = true;
357  }
358  }
359 
360  ON_Interval u_interval(min.x, max.x);
361  ON_Interval v_interval(min.y, max.y);
362  if (!surface_GetBoundingBox(surf, u_interval,v_interval,bbox,growcurrent)) {
363  return false;
364  }
365 
366  return true;
367 }
368 
369 
370 ON_BOOL32
372  const ON_3dPoint& p,
373  ON_BoundingBox &bbox,
374  double &min_distance,
375  double &max_distance
376  )
377 {
378  min_distance = bbox.MinimumDistanceTo(p);
379  max_distance = bbox.MaximumDistanceTo(p);
380  return true;
381 }
382 
383 
384 ON_BOOL32
386  const ON_Surface *surf,
387  const ON_3dPoint& p,
388  const ON_Interval &u_interval,
389  const ON_Interval &v_interval,
390  double &min_distance,
391  double &max_distance
392  )
393 {
394  ON_BoundingBox bbox;
395 
396  if (surface_GetBoundingBox(surf,u_interval,v_interval,bbox, false)) {
397  min_distance = bbox.MinimumDistanceTo(p);
398 
399  max_distance = bbox.MaximumDistanceTo(p);
400  return true;
401  }
402  return false;
403 }
404 
405 
406 double
407 surface_GetOptimalNormalUSplit(const ON_Surface *surf, const ON_Interval &u_interval, const ON_Interval &v_interval,double tol)
408 {
409  ON_3dVector normal[4];
410  double u = u_interval.Mid();
411 
412  if ((normal[0] = surf->NormalAt(u_interval.m_t[0],v_interval.m_t[0])) &&
413  (normal[2] = surf->NormalAt(u_interval.m_t[0],v_interval.m_t[1]))) {
414  double step = u_interval.Length()/2.0;
415  double stepdir = 1.0;
416  u = u_interval.m_t[0] + stepdir * step;
417 
418  while (step > tol) {
419  if ((normal[1] = surf->NormalAt(u,v_interval.m_t[0])) &&
420  (normal[3] = surf->NormalAt(u,v_interval.m_t[1]))) {
421  double udot_1 = normal[0] * normal[1];
422  double udot_2 = normal[2] * normal[3];
423  if ((udot_1 < 0.0) || (udot_2 < 0.0)) {
424  stepdir = -1.0;
425  } else {
426  stepdir = 1.0;
427  }
428  step = step / 2.0;
429  u = u + stepdir * step;
430  }
431  }
432  }
433  return u;
434 }
435 
436 
437 double
438 surface_GetOptimalNormalVSplit(const ON_Surface *surf, const ON_Interval &u_interval, const ON_Interval &v_interval,double tol)
439 {
440  ON_3dVector normal[4];
441  double v = v_interval.Mid();
442 
443  if ((normal[0] = surf->NormalAt(u_interval.m_t[0],v_interval.m_t[0])) &&
444  (normal[1] = surf->NormalAt(u_interval.m_t[1],v_interval.m_t[0]))) {
445  double step = v_interval.Length()/2.0;
446  double stepdir = 1.0;
447  v = v_interval.m_t[0] + stepdir * step;
448 
449  while (step > tol) {
450  if ((normal[2] = surf->NormalAt(u_interval.m_t[0],v)) &&
451  (normal[3] = surf->NormalAt(u_interval.m_t[1],v))) {
452  double vdot_1 = normal[0] * normal[2];
453  double vdot_2 = normal[1] * normal[3];
454  if ((vdot_1 < 0.0) || (vdot_2 < 0.0)) {
455  stepdir = -1.0;
456  } else {
457  stepdir = 1.0;
458  }
459  step = step / 2.0;
460  v = v + stepdir * step;
461  }
462  }
463  }
464  return v;
465 }
466 
467 
468 //forward for cyclic
469 double surface_GetClosestPoint3dFirstOrderByRange(const ON_Surface *surf,const ON_3dPoint& p,const ON_Interval& u_range,
470  const ON_Interval& v_range,double current_closest_dist,ON_2dPoint& p2d,ON_3dPoint& p3d,double same_point_tol, double within_distance_tol,int level);
471 
473  const ON_3dPoint& p, const ON_Interval &u_interval, double u, const ON_Interval &v_interval, double v,
474  double current_closest_dist, ON_2dPoint& p2d, ON_3dPoint& p3d,
475  double same_point_tol, double within_distance_tol, int level)
476 {
477  double min_distance = 0;
478  double max_distance = 0;
479  ON_Interval new_u_interval = u_interval;
480  ON_Interval new_v_interval = v_interval;
481 
482  for (int iu = 0; iu < 2; iu++) {
483  new_u_interval.m_t[iu] = u_interval.m_t[iu];
484  new_u_interval.m_t[1 - iu] = u;
485  for (int iv = 0; iv < 2; iv++) {
486  new_v_interval.m_t[iv] = v_interval.m_t[iv];
487  new_v_interval.m_t[1 - iv] = v;
488  if (surface_GetIntervalMinMaxDistance(surf, p, new_u_interval,new_v_interval, min_distance, max_distance)) {
489  double distance = DBL_MAX;
490  if ((min_distance < current_closest_dist) && NEAR_ZERO(min_distance,within_distance_tol)) {
491  /////////////////////////////////////////
492  // Could check normals and CV angles here
493  /////////////////////////////////////////
494  ON_3dVector normal[4];
495  if ((normal[0] = surf->NormalAt(new_u_interval.m_t[0],new_v_interval.m_t[0])) &&
496  (normal[1] = surf->NormalAt(new_u_interval.m_t[1],new_v_interval.m_t[0])) &&
497  (normal[2] = surf->NormalAt(new_u_interval.m_t[0],new_v_interval.m_t[1])) &&
498  (normal[3] = surf->NormalAt(new_u_interval.m_t[1],new_v_interval.m_t[1]))) {
499 
500  double udot_1 = normal[0] * normal[1];
501  double udot_2 = normal[2] * normal[3];
502  double vdot_1 = normal[0] * normal[2];
503  double vdot_2 = normal[1] * normal[3];
504 
505  if ((udot_1 < 0.0) || (udot_2 < 0.0) || (vdot_1 < 0.0) || (vdot_2 < 0.0)) {
506  double u_split,v_split;
507  if ((udot_1 < 0.0) || (udot_2 < 0.0)) {
508  //get optimal U split
509  u_split = surface_GetOptimalNormalUSplit(surf,new_u_interval,new_v_interval,same_point_tol);
510  } else {
511  u_split = new_u_interval.Mid();
512  }
513  if ((vdot_1 < 0.0) || (vdot_2 < 0.0)) {
514  //get optimal V split
515  v_split = surface_GetOptimalNormalVSplit(surf,new_u_interval,new_v_interval,same_point_tol);
516  } else {
517  v_split = new_v_interval.Mid();
518  }
519  distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,new_u_interval,u_split,new_v_interval,v_split,current_closest_dist,p2d,p3d,same_point_tol,within_distance_tol,level++);
520  if (distance < current_closest_dist) {
521  current_closest_dist = distance;
522  if (current_closest_dist < same_point_tol)
523  return current_closest_dist;
524  }
525  } else {
526  distance = surface_GetClosestPoint3dFirstOrderByRange(surf,p,new_u_interval,new_v_interval,current_closest_dist,p2d,p3d,same_point_tol,within_distance_tol,level++);
527  if (distance < current_closest_dist) {
528  current_closest_dist = distance;
529  if (current_closest_dist < same_point_tol)
530  return current_closest_dist;
531  }
532  }
533  }
534  }
535  }
536  }
537  }
538  return current_closest_dist;
539 }
540 
541 
542 double
544  const ON_Surface *surf,
545  const ON_3dPoint& p,
546  const ON_Interval& u_range,
547  const ON_Interval& v_range,
548  double current_closest_dist,
549  ON_2dPoint& p2d,
550  ON_3dPoint& p3d,
551  double same_point_tol,
552  double within_distance_tol,
553  int level
554  )
555 {
556  ON_3dPoint p0;
557  ON_2dPoint p2d0;
558  ON_3dVector ds, dt, dss, dst, dtt;
559  ON_2dPoint working_p2d;
560  ON_3dPoint working_p3d;
561  ON_3dPoint P;
562  ON_3dVector Ds, Dt, Dss, Dst, Dtt;
563  bool notdone = true;
564  double previous_distance = DBL_MAX;
565  double distance;
566  int errcnt = 0;
567 
568  p2d0.x = u_range.Mid();
569  p2d0.y = v_range.Mid();
570 
571  while (notdone && (surf->Ev2Der(p2d0.x, p2d0.y, p0, ds, dt, dss, dst, dtt))) {
572  if ((distance = p0.DistanceTo(p)) >= previous_distance) {
573  if (++errcnt <= 10) {
574  p2d0 = (p2d0 + working_p2d) / 2.0;
575  continue;
576  } else {
577  ///////////////////////////
578  // Don't Subdivide just not getting any closer
579  ///////////////////////////
580  /*
581  double distance =
582  surface_GetClosestPoint3dFirstOrderSubdivision(surf, p,
583  u_range, u_range.Mid(), v_range, v_range.Mid(),
584  current_closest_dist, p2d, p3d, tol, level++);
585  if (distance < current_closest_dist) {
586  current_closest_dist = distance;
587  if (current_closest_dist < tol)
588  return current_closest_dist;
589  }
590  */
591  break;
592  }
593  } else {
594  if (distance < current_closest_dist) {
595  current_closest_dist = distance;
596  p3d = p0;
597  p2d = p2d0;
598  if (current_closest_dist < same_point_tol)
599  return current_closest_dist;
600  }
601  previous_distance = distance;
602  working_p3d = p0;
603  working_p2d = p2d0;
604  errcnt = 0;
605  }
606  ON_3dVector N = ON_CrossProduct(ds, dt);
607  N.Unitize();
608  ON_Plane plane(p0, N);
609  ON_3dPoint q = plane.ClosestPointTo(p);
610  ON_2dVector pullback;
611  ON_3dVector vector = q - p0;
612  double vlength = vector.Length();
613 
614  if (vlength > 0.0) {
615  if (ON_Pullback3dVector(vector, 0.0, ds, dt, dss, dst, dtt,
616  pullback)) {
617  p2d0 = p2d0 + pullback;
618  if (!u_range.Includes(p2d0.x, false)) {
619  int i = (u_range.m_t[0] <= u_range.m_t[1]) ? 0 : 1;
620  p2d0.x =
621  (p2d0.x < u_range.m_t[i]) ?
622  u_range.m_t[i] : u_range.m_t[1 - i];
623  }
624  if (!v_range.Includes(p2d0.y, false)) {
625  int i = (v_range.m_t[0] <= v_range.m_t[1]) ? 0 : 1;
626  p2d0.y =
627  (p2d0.y < v_range.m_t[i]) ?
628  v_range.m_t[i] : v_range.m_t[1 - i];
629  }
630  } else {
631  ///////////////////////////
632  // Subdivide and go again
633  ///////////////////////////
634  notdone = false;
635  distance =
637  u_range, u_range.Mid(), v_range, v_range.Mid(),
638  current_closest_dist, p2d, p3d, same_point_tol, within_distance_tol, level++);
639  if (distance < current_closest_dist) {
640  current_closest_dist = distance;
641  if (current_closest_dist < same_point_tol)
642  return current_closest_dist;
643  }
644  break;
645  }
646  } else {
647  // can't get any closer
648  notdone = false;
649  break;
650  }
651  }
652  if (previous_distance < current_closest_dist) {
653  current_closest_dist = previous_distance;
654  p3d = working_p3d;
655  p2d = working_p2d;
656  }
657 
658  return current_closest_dist;
659 }
660 
661 
663  const ON_Surface *surf,
664  const ON_3dPoint& p,
665  ON_2dPoint& p2d,
666  ON_3dPoint& p3d,
667  double &current_distance,
668  int quadrant, // optional - determines which quadrant to evaluate from
669  // 0 = default
670  // 1 from NE quadrant
671  // 2 from NW quadrant
672  // 3 from SW quadrant
673  // 4 from SE quadrant
674  double same_point_tol,
675  double within_distance_tol
676  )
677 {
678  ON_3dPoint p0;
679  ON_2dPoint p2d0;
680  ON_3dVector ds, dt, dss, dst, dtt;
681  ON_3dVector T, K;
682  bool rc = false;
683 
684  static const ON_Surface *prev_surface = NULL;
685  static int prev_u_spancnt = 0;
686  static int u_spancnt = 0;
687  static int v_spancnt = 0;
688  static double *uspan = NULL;
689  static double *vspan = NULL;
690  static ON_BoundingBox **bbox = NULL;
691  static double umid = 0.0;
692  static int umid_index = 0;
693  static double vmid = 0.0;
694  static int vmid_index = 0;
695 
696  current_distance = DBL_MAX;
697 
698  int prec = std::cerr.precision();
699  std::cerr.precision(15);
700 
701 
702  if (prev_surface != surf) {
703  if (uspan)
704  delete [] uspan;
705  if (vspan)
706  delete [] vspan;
707  if (bbox) {
708  for( int i = 0 ; i < (prev_u_spancnt + 2) ; i++ )
709  delete [] bbox[i] ;
710  delete [] bbox;
711  }
712  u_spancnt = prev_u_spancnt = surf->SpanCount(0);
713  v_spancnt = surf->SpanCount(1);
714  // adding 2 here because going to divide at midpoint
715  uspan = new double[u_spancnt + 2];
716  vspan = new double[v_spancnt + 2];
717  bbox = new ON_BoundingBox *[(u_spancnt + 2)];
718  for( int i = 0 ; i < (u_spancnt + 2) ; i++ )
719  bbox[i] = new ON_BoundingBox [v_spancnt + 2];
720 
721  if (surf->GetSpanVector(0, uspan) && surf->GetSpanVector(1, vspan)) {
722  prev_surface = surf;
723  umid = surf->Domain(0).Mid();
724  umid_index = u_spancnt/2;
725  for (int u_span_index = 0; u_span_index < u_spancnt + 1;u_span_index++) {
726  if (NEAR_EQUAL(uspan[u_span_index],umid,same_point_tol)) {
727  umid_index = u_span_index;
728  break;
729  } else if (uspan[u_span_index] > umid) {
730  for (u_span_index = u_spancnt + 1; u_span_index > 0;u_span_index--) {
731  if (uspan[u_span_index-1] < umid) {
732  uspan[u_span_index] = umid;
733  umid_index = u_span_index;
734  u_spancnt++;
735  u_span_index = u_spancnt+1;
736  break;
737  } else {
738  uspan[u_span_index] = uspan[u_span_index-1];
739  }
740  }
741  }
742  }
743  vmid = surf->Domain(1).Mid();
744  vmid_index = v_spancnt/2;
745  for (int v_span_index = 0; v_span_index < v_spancnt + 1;v_span_index++) {
746  if (NEAR_EQUAL(vspan[v_span_index],vmid,same_point_tol)) {
747  vmid_index = v_span_index;
748  break;
749  } else if (vspan[v_span_index] > vmid) {
750  for (v_span_index = v_spancnt + 1; v_span_index > 0;v_span_index--) {
751  if (vspan[v_span_index-1] < vmid) {
752  vspan[v_span_index] = vmid;
753  vmid_index = v_span_index;
754  v_spancnt++;
755  v_span_index = v_spancnt+1;
756  break;
757  } else {
758  vspan[v_span_index] = vspan[v_span_index-1];
759  }
760  }
761  }
762  }
763  for (int u_span_index = 1; u_span_index < u_spancnt + 1;
764  u_span_index++) {
765  for (int v_span_index = 1; v_span_index < v_spancnt + 1;
766  v_span_index++) {
767  ON_Interval u_interval(uspan[u_span_index - 1],
768  uspan[u_span_index]);
769  ON_Interval v_interval(vspan[v_span_index - 1],
770  vspan[v_span_index]);
771 
772  if (!surface_GetBoundingBox(surf,u_interval,v_interval,bbox[u_span_index-1][v_span_index-1], false)) {
773  std::cerr << "Error computing bounding box for surface interval" << std::endl;
774  }
775  }
776  }
777  } else {
778  prev_surface = NULL;
779  }
780  }
781  if (prev_surface == surf) {
782  if (quadrant == 0) {
783  for (int u_span_index = 1; u_span_index < u_spancnt + 1;
784  u_span_index++) {
785  for (int v_span_index = 1; v_span_index < v_spancnt + 1;
786  v_span_index++) {
787  ON_Interval u_interval(uspan[u_span_index - 1],
788  uspan[u_span_index]);
789  ON_Interval v_interval(vspan[v_span_index - 1],
790  vspan[v_span_index]);
791  double min_distance,max_distance;
792 
793  int level = 1;
794  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
795  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
796  /////////////////////////////////////////
797  // Could check normals and CV angles here
798  /////////////////////////////////////////
799  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
800  if (distance < current_distance) {
801  current_distance = distance;
802  if (current_distance < same_point_tol) {
803  rc = true;
804  goto cleanup;
805  }
806  }
807  }
808  }
809  }
810  }
811  if (current_distance < within_distance_tol) {
812  rc = true;
813  goto cleanup;
814  }
815  } else if (quadrant == 1) {
816  if (surf->IsClosed(0)) { //NE,SE,NW.SW
817  // 1 from NE quadrant
818  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
819  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
820  ON_Interval u_interval(uspan[u_span_index - 1],
821  uspan[u_span_index]);
822  ON_Interval v_interval(vspan[v_span_index - 1],
823  vspan[v_span_index]);
824  double min_distance,max_distance;
825 
826  int level = 1;
827  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
828  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
829  /////////////////////////////////////////
830  // Could check normals and CV angles here
831  /////////////////////////////////////////
832  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
833  if (distance < current_distance) {
834  current_distance = distance;
835  if (current_distance < same_point_tol) {
836  rc = true;
837  goto cleanup;
838  }
839  }
840  }
841  }
842  }
843  }
844  // 4 from SE quadrant
845  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
846  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
847  ON_Interval u_interval(uspan[u_span_index - 1],
848  uspan[u_span_index]);
849  ON_Interval v_interval(vspan[v_span_index - 1],
850  vspan[v_span_index]);
851  double min_distance,max_distance;
852 
853  int level = 1;
854  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
855  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
856  /////////////////////////////////////////
857  // Could check normals and CV angles here
858  /////////////////////////////////////////
859  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
860  if (distance < current_distance) {
861  current_distance = distance;
862  if (current_distance < same_point_tol) {
863  rc = true;
864  goto cleanup;
865  }
866  }
867  }
868  }
869  }
870  }
871  // 2 from NW quadrant
872  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
873  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
874  ON_Interval u_interval(uspan[u_span_index - 1],
875  uspan[u_span_index]);
876  ON_Interval v_interval(vspan[v_span_index - 1],
877  vspan[v_span_index]);
878  double min_distance,max_distance;
879 
880  int level = 1;
881  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
882  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
883  /////////////////////////////////////////
884  // Could check normals and CV angles here
885  /////////////////////////////////////////
886  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
887  if (distance < current_distance) {
888  current_distance = distance;
889  if (current_distance < same_point_tol) {
890  rc = true;
891  goto cleanup;
892  }
893  }
894  }
895  }
896  }
897  }
898  // 3 from SW quadrant
899  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
900  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
901  ON_Interval u_interval(uspan[u_span_index - 1],
902  uspan[u_span_index]);
903  ON_Interval v_interval(vspan[v_span_index - 1],
904  vspan[v_span_index]);
905  double min_distance,max_distance;
906 
907  int level = 1;
908  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
909  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
910  /////////////////////////////////////////
911  // Could check normals and CV angles here
912  /////////////////////////////////////////
913  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
914  if (distance < current_distance) {
915  current_distance = distance;
916  if (current_distance < same_point_tol) {
917  rc = true;
918  goto cleanup;
919  }
920  }
921  }
922  }
923  }
924  }
925  } else { //NE,NW,SW,SE
926  // 1 from NE quadrant
927  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
928  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
929  ON_Interval u_interval(uspan[u_span_index - 1],
930  uspan[u_span_index]);
931  ON_Interval v_interval(vspan[v_span_index - 1],
932  vspan[v_span_index]);
933  double min_distance,max_distance;
934 
935  int level = 1;
936  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
937  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
938  /////////////////////////////////////////
939  // Could check normals and CV angles here
940  /////////////////////////////////////////
941  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
942  if (distance < current_distance) {
943  current_distance = distance;
944  if (current_distance < same_point_tol) {
945  rc = true;
946  goto cleanup;
947  }
948  }
949  }
950  }
951  }
952  }
953  // 2 from NW quadrant
954  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
955  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
956  ON_Interval u_interval(uspan[u_span_index - 1],
957  uspan[u_span_index]);
958  ON_Interval v_interval(vspan[v_span_index - 1],
959  vspan[v_span_index]);
960  double min_distance,max_distance;
961 
962  int level = 1;
963  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
964  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
965  /////////////////////////////////////////
966  // Could check normals and CV angles here
967  /////////////////////////////////////////
968  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
969  if (distance < current_distance) {
970  current_distance = distance;
971  if (current_distance < same_point_tol) {
972  rc = true;
973  goto cleanup;
974  }
975  }
976  }
977  }
978  }
979  }
980  // 3 from SW quadrant
981  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
982  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
983  ON_Interval u_interval(uspan[u_span_index - 1],
984  uspan[u_span_index]);
985  ON_Interval v_interval(vspan[v_span_index - 1],
986  vspan[v_span_index]);
987  double min_distance,max_distance;
988 
989  int level = 1;
990  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
991  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
992  /////////////////////////////////////////
993  // Could check normals and CV angles here
994  /////////////////////////////////////////
995  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
996  if (distance < current_distance) {
997  current_distance = distance;
998  if (current_distance < same_point_tol) {
999  rc = true;
1000  goto cleanup;
1001  }
1002  }
1003  }
1004  }
1005  }
1006  }
1007  // 4 from SE quadrant
1008  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1009  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1010  ON_Interval u_interval(uspan[u_span_index - 1],
1011  uspan[u_span_index]);
1012  ON_Interval v_interval(vspan[v_span_index - 1],
1013  vspan[v_span_index]);
1014  double min_distance,max_distance;
1015 
1016  int level = 1;
1017  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1018  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1019  /////////////////////////////////////////
1020  // Could check normals and CV angles here
1021  /////////////////////////////////////////
1022  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1023  if (distance < current_distance) {
1024  current_distance = distance;
1025  if (current_distance < same_point_tol) {
1026  rc = true;
1027  goto cleanup;
1028  }
1029  }
1030  }
1031  }
1032  }
1033  }
1034  }
1035  if (current_distance < within_distance_tol) {
1036  rc = true;
1037  goto cleanup;
1038  }
1039  } else if (quadrant == 2) {
1040  if (surf->IsClosed(0)) { // NW,SW,NE,SE
1041  // 2 from NW quadrant
1042  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1043  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1044  ON_Interval u_interval(uspan[u_span_index - 1],
1045  uspan[u_span_index]);
1046  ON_Interval v_interval(vspan[v_span_index - 1],
1047  vspan[v_span_index]);
1048  double min_distance,max_distance;
1049 
1050  int level = 1;
1051  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1052  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1053  /////////////////////////////////////////
1054  // Could check normals and CV angles here
1055  /////////////////////////////////////////
1056  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1057  if (distance < current_distance) {
1058  current_distance = distance;
1059  if (current_distance < same_point_tol) {
1060  rc = true;
1061  goto cleanup;
1062  }
1063  }
1064  }
1065  }
1066  }
1067  }
1068  // 3 from SW quadrant
1069  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1070  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1071  ON_Interval u_interval(uspan[u_span_index - 1],
1072  uspan[u_span_index]);
1073  ON_Interval v_interval(vspan[v_span_index - 1],
1074  vspan[v_span_index]);
1075  double min_distance,max_distance;
1076 
1077  int level = 1;
1078  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1079  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1080  /////////////////////////////////////////
1081  // Could check normals and CV angles here
1082  /////////////////////////////////////////
1083  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1084  if (distance < current_distance) {
1085  current_distance = distance;
1086  if (current_distance < same_point_tol) {
1087  rc = true;
1088  goto cleanup;
1089  }
1090  }
1091  }
1092  }
1093  }
1094  }
1095  // 1 from NE quadrant
1096  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1097  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1098  ON_Interval u_interval(uspan[u_span_index - 1],
1099  uspan[u_span_index]);
1100  ON_Interval v_interval(vspan[v_span_index - 1],
1101  vspan[v_span_index]);
1102  double min_distance,max_distance;
1103 
1104  int level = 1;
1105  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1106  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1107  /////////////////////////////////////////
1108  // Could check normals and CV angles here
1109  /////////////////////////////////////////
1110  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1111  if (distance < current_distance) {
1112  current_distance = distance;
1113  if (current_distance < same_point_tol) {
1114  rc = true;
1115  goto cleanup;
1116  }
1117  }
1118  }
1119  }
1120  }
1121  }
1122  // 4 from SE quadrant
1123  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1124  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1125  ON_Interval u_interval(uspan[u_span_index - 1],
1126  uspan[u_span_index]);
1127  ON_Interval v_interval(vspan[v_span_index - 1],
1128  vspan[v_span_index]);
1129  double min_distance,max_distance;
1130 
1131  int level = 1;
1132  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1133  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1134  /////////////////////////////////////////
1135  // Could check normals and CV angles here
1136  /////////////////////////////////////////
1137  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1138  if (distance < current_distance) {
1139  current_distance = distance;
1140  if (current_distance < same_point_tol) {
1141  rc = true;
1142  goto cleanup;
1143  }
1144  }
1145  }
1146  }
1147  }
1148  }
1149  } else { // NW,NE,SW,SE
1150  // 2 from NW quadrant
1151  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1152  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1153  ON_Interval u_interval(uspan[u_span_index - 1],
1154  uspan[u_span_index]);
1155  ON_Interval v_interval(vspan[v_span_index - 1],
1156  vspan[v_span_index]);
1157  double min_distance,max_distance;
1158 
1159  int level = 1;
1160  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1161  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1162  /////////////////////////////////////////
1163  // Could check normals and CV angles here
1164  /////////////////////////////////////////
1165  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1166  if (distance < current_distance) {
1167  current_distance = distance;
1168  if (current_distance < same_point_tol) {
1169  rc = true;
1170  goto cleanup;
1171  }
1172  }
1173  }
1174  }
1175  }
1176  }
1177  // 1 from NE quadrant
1178  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1179  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1180  ON_Interval u_interval(uspan[u_span_index - 1],
1181  uspan[u_span_index]);
1182  ON_Interval v_interval(vspan[v_span_index - 1],
1183  vspan[v_span_index]);
1184  double min_distance,max_distance;
1185 
1186  int level = 1;
1187  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1188  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1189  /////////////////////////////////////////
1190  // Could check normals and CV angles here
1191  /////////////////////////////////////////
1192  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1193  if (distance < current_distance) {
1194  current_distance = distance;
1195  if (current_distance < same_point_tol) {
1196  rc = true;
1197  goto cleanup;
1198  }
1199  }
1200  }
1201  }
1202  }
1203  }
1204  // 3 from SW quadrant
1205  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1206  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1207  ON_Interval u_interval(uspan[u_span_index - 1],
1208  uspan[u_span_index]);
1209  ON_Interval v_interval(vspan[v_span_index - 1],
1210  vspan[v_span_index]);
1211  double min_distance,max_distance;
1212 
1213  int level = 1;
1214  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1215  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1216  /////////////////////////////////////////
1217  // Could check normals and CV angles here
1218  /////////////////////////////////////////
1219  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1220  if (distance < current_distance) {
1221  current_distance = distance;
1222  if (current_distance < same_point_tol) {
1223  rc = true;
1224  goto cleanup;
1225  }
1226  }
1227  }
1228  }
1229  }
1230  }
1231  // 4 from SE quadrant
1232  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1233  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1234  ON_Interval u_interval(uspan[u_span_index - 1],
1235  uspan[u_span_index]);
1236  ON_Interval v_interval(vspan[v_span_index - 1],
1237  vspan[v_span_index]);
1238  double min_distance,max_distance;
1239 
1240  int level = 1;
1241  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1242  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1243  /////////////////////////////////////////
1244  // Could check normals and CV angles here
1245  /////////////////////////////////////////
1246  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1247  if (distance < current_distance) {
1248  current_distance = distance;
1249  if (current_distance < same_point_tol) {
1250  rc = true;
1251  goto cleanup;
1252  }
1253  }
1254  }
1255  }
1256  }
1257  }
1258  }
1259  if (current_distance < within_distance_tol) {
1260  rc = true;
1261  goto cleanup;
1262  }
1263  } else if (quadrant == 3) {
1264  if (surf->IsClosed(0)) { // SW,NW,SE,NE
1265  // 3 from SW quadrant
1266  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1267  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1268  ON_Interval u_interval(uspan[u_span_index - 1],
1269  uspan[u_span_index]);
1270  ON_Interval v_interval(vspan[v_span_index - 1],
1271  vspan[v_span_index]);
1272  double min_distance,max_distance;
1273 
1274  int level = 1;
1275  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1276  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1277  /////////////////////////////////////////
1278  // Could check normals and CV angles here
1279  /////////////////////////////////////////
1280  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1281  if (distance < current_distance) {
1282  current_distance = distance;
1283  if (current_distance < same_point_tol) {
1284  rc = true;
1285  goto cleanup;
1286  }
1287  }
1288  }
1289  }
1290  }
1291  }
1292  // 2 from NW quadrant
1293  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1294  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1295  ON_Interval u_interval(uspan[u_span_index - 1],
1296  uspan[u_span_index]);
1297  ON_Interval v_interval(vspan[v_span_index - 1],
1298  vspan[v_span_index]);
1299  double min_distance,max_distance;
1300 
1301  int level = 1;
1302  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1303  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1304  /////////////////////////////////////////
1305  // Could check normals and CV angles here
1306  /////////////////////////////////////////
1307  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1308  if (distance < current_distance) {
1309  current_distance = distance;
1310  if (current_distance < same_point_tol) {
1311  rc = true;
1312  goto cleanup;
1313  }
1314  }
1315  }
1316  }
1317  }
1318  }
1319  // 4 from SE quadrant
1320  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1321  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1322  ON_Interval u_interval(uspan[u_span_index - 1],
1323  uspan[u_span_index]);
1324  ON_Interval v_interval(vspan[v_span_index - 1],
1325  vspan[v_span_index]);
1326  double min_distance,max_distance;
1327 
1328  int level = 1;
1329  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1330  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1331  /////////////////////////////////////////
1332  // Could check normals and CV angles here
1333  /////////////////////////////////////////
1334  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1335  if (distance < current_distance) {
1336  current_distance = distance;
1337  if (current_distance < same_point_tol) {
1338  rc = true;
1339  goto cleanup;
1340  }
1341  }
1342  }
1343  }
1344  }
1345  }
1346  // 1 from NE quadrant
1347  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1348  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1349  ON_Interval u_interval(uspan[u_span_index - 1],
1350  uspan[u_span_index]);
1351  ON_Interval v_interval(vspan[v_span_index - 1],
1352  vspan[v_span_index]);
1353  double min_distance,max_distance;
1354 
1355  int level = 1;
1356  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1357  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1358  /////////////////////////////////////////
1359  // Could check normals and CV angles here
1360  /////////////////////////////////////////
1361  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1362  if (distance < current_distance) {
1363  current_distance = distance;
1364  if (current_distance < same_point_tol) {
1365  rc = true;
1366  goto cleanup;
1367  }
1368  }
1369  }
1370  }
1371  }
1372  }
1373  } else { // SW,SE,NW,NE
1374  // 3 from SW quadrant
1375  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1376  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1377  ON_Interval u_interval(uspan[u_span_index - 1],
1378  uspan[u_span_index]);
1379  ON_Interval v_interval(vspan[v_span_index - 1],
1380  vspan[v_span_index]);
1381  double min_distance,max_distance;
1382 
1383  int level = 1;
1384  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1385  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1386  /////////////////////////////////////////
1387  // Could check normals and CV angles here
1388  /////////////////////////////////////////
1389  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1390  if (distance < current_distance) {
1391  current_distance = distance;
1392  if (current_distance < same_point_tol) {
1393  rc = true;
1394  goto cleanup;
1395  }
1396  }
1397  }
1398  }
1399  }
1400  }
1401  // 4 from SE quadrant
1402  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1403  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1404  ON_Interval u_interval(uspan[u_span_index - 1],
1405  uspan[u_span_index]);
1406  ON_Interval v_interval(vspan[v_span_index - 1],
1407  vspan[v_span_index]);
1408  double min_distance,max_distance;
1409 
1410  int level = 1;
1411  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1412  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1413  /////////////////////////////////////////
1414  // Could check normals and CV angles here
1415  /////////////////////////////////////////
1416  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1417  if (distance < current_distance) {
1418  current_distance = distance;
1419  if (current_distance < same_point_tol) {
1420  rc = true;
1421  goto cleanup;
1422  }
1423  }
1424  }
1425  }
1426  }
1427  }
1428  // 2 from NW quadrant
1429  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1430  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1431  ON_Interval u_interval(uspan[u_span_index - 1],
1432  uspan[u_span_index]);
1433  ON_Interval v_interval(vspan[v_span_index - 1],
1434  vspan[v_span_index]);
1435  double min_distance,max_distance;
1436 
1437  int level = 1;
1438  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1439  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1440  /////////////////////////////////////////
1441  // Could check normals and CV angles here
1442  /////////////////////////////////////////
1443  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1444  if (distance < current_distance) {
1445  current_distance = distance;
1446  if (current_distance < same_point_tol) {
1447  rc = true;
1448  goto cleanup;
1449  }
1450  }
1451  }
1452  }
1453  }
1454  }
1455  // 1 from NE quadrant
1456  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1457  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1458  ON_Interval u_interval(uspan[u_span_index - 1],
1459  uspan[u_span_index]);
1460  ON_Interval v_interval(vspan[v_span_index - 1],
1461  vspan[v_span_index]);
1462  double min_distance,max_distance;
1463 
1464  int level = 1;
1465  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1466  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1467  /////////////////////////////////////////
1468  // Could check normals and CV angles here
1469  /////////////////////////////////////////
1470  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1471  if (distance < current_distance) {
1472  current_distance = distance;
1473  if (current_distance < same_point_tol) {
1474  rc = true;
1475  goto cleanup;
1476  }
1477  }
1478  }
1479  }
1480  }
1481  }
1482  }
1483  if (current_distance < within_distance_tol) {
1484  rc = true;
1485  goto cleanup;
1486  }
1487  } else if (quadrant == 4) {
1488  if (surf->IsClosed(0)) { // SE,NE,SW,NW
1489  // 4 from SE quadrant
1490  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1491  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1492  ON_Interval u_interval(uspan[u_span_index - 1],
1493  uspan[u_span_index]);
1494  ON_Interval v_interval(vspan[v_span_index - 1],
1495  vspan[v_span_index]);
1496  double min_distance,max_distance;
1497 
1498  int level = 1;
1499  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1500  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1501  /////////////////////////////////////////
1502  // Could check normals and CV angles here
1503  /////////////////////////////////////////
1504  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1505  if (distance < current_distance) {
1506  current_distance = distance;
1507  if (current_distance < same_point_tol) {
1508  rc = true;
1509  goto cleanup;
1510  }
1511  }
1512  }
1513  }
1514  }
1515  }
1516  // 1 from NE quadrant
1517  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1518  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1519  ON_Interval u_interval(uspan[u_span_index - 1],
1520  uspan[u_span_index]);
1521  ON_Interval v_interval(vspan[v_span_index - 1],
1522  vspan[v_span_index]);
1523  double min_distance,max_distance;
1524 
1525  int level = 1;
1526  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1527  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1528  /////////////////////////////////////////
1529  // Could check normals and CV angles here
1530  /////////////////////////////////////////
1531  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1532  if (distance < current_distance) {
1533  current_distance = distance;
1534  if (current_distance < same_point_tol) {
1535  rc = true;
1536  goto cleanup;
1537  }
1538  }
1539  }
1540  }
1541  }
1542  }
1543  // 3 from SW quadrant
1544  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1545  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1546  ON_Interval u_interval(uspan[u_span_index - 1],
1547  uspan[u_span_index]);
1548  ON_Interval v_interval(vspan[v_span_index - 1],
1549  vspan[v_span_index]);
1550  double min_distance,max_distance;
1551 
1552  int level = 1;
1553  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1554  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1555  /////////////////////////////////////////
1556  // Could check normals and CV angles here
1557  /////////////////////////////////////////
1558  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1559  if (distance < current_distance) {
1560  current_distance = distance;
1561  if (current_distance < same_point_tol) {
1562  rc = true;
1563  goto cleanup;
1564  }
1565  }
1566  }
1567  }
1568  }
1569  }
1570  // 2 from NW quadrant
1571  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1572  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1573  ON_Interval u_interval(uspan[u_span_index - 1],
1574  uspan[u_span_index]);
1575  ON_Interval v_interval(vspan[v_span_index - 1],
1576  vspan[v_span_index]);
1577  double min_distance,max_distance;
1578 
1579  int level = 1;
1580  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1581  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1582  /////////////////////////////////////////
1583  // Could check normals and CV angles here
1584  /////////////////////////////////////////
1585  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1586  if (distance < current_distance) {
1587  current_distance = distance;
1588  if (current_distance < same_point_tol) {
1589  rc = true;
1590  goto cleanup;
1591  }
1592  }
1593  }
1594  }
1595  }
1596  }
1597  } else { // SE,SW,NE,NW
1598  // 4 from SE quadrant
1599  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1600  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1601  ON_Interval u_interval(uspan[u_span_index - 1],
1602  uspan[u_span_index]);
1603  ON_Interval v_interval(vspan[v_span_index - 1],
1604  vspan[v_span_index]);
1605  double min_distance,max_distance;
1606 
1607  int level = 1;
1608  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1609  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1610  /////////////////////////////////////////
1611  // Could check normals and CV angles here
1612  /////////////////////////////////////////
1613  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1614  if (distance < current_distance) {
1615  current_distance = distance;
1616  if (current_distance < same_point_tol) {
1617  rc = true;
1618  goto cleanup;
1619  }
1620  }
1621  }
1622  }
1623  }
1624  }
1625  // 3 from SW quadrant
1626  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1627  for (int v_span_index = 1; v_span_index < vmid_index+1; v_span_index++) {
1628  ON_Interval u_interval(uspan[u_span_index - 1],
1629  uspan[u_span_index]);
1630  ON_Interval v_interval(vspan[v_span_index - 1],
1631  vspan[v_span_index]);
1632  double min_distance,max_distance;
1633 
1634  int level = 1;
1635  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1636  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1637  /////////////////////////////////////////
1638  // Could check normals and CV angles here
1639  /////////////////////////////////////////
1640  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1641  if (distance < current_distance) {
1642  current_distance = distance;
1643  if (current_distance < same_point_tol) {
1644  rc = true;
1645  goto cleanup;
1646  }
1647  }
1648  }
1649  }
1650  }
1651  }
1652  // 1 from NE quadrant
1653  for (int u_span_index = umid_index+1; u_span_index < u_spancnt + 1; u_span_index++) {
1654  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1655  ON_Interval u_interval(uspan[u_span_index - 1],
1656  uspan[u_span_index]);
1657  ON_Interval v_interval(vspan[v_span_index - 1],
1658  vspan[v_span_index]);
1659  double min_distance,max_distance;
1660 
1661  int level = 1;
1662  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1663  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1664  /////////////////////////////////////////
1665  // Could check normals and CV angles here
1666  /////////////////////////////////////////
1667  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1668  if (distance < current_distance) {
1669  current_distance = distance;
1670  if (current_distance < same_point_tol) {
1671  rc = true;
1672  goto cleanup;
1673  }
1674  }
1675  }
1676  }
1677  }
1678  }
1679  // 2 from NW quadrant
1680  for (int u_span_index = 1; u_span_index < umid_index+1; u_span_index++) {
1681  for (int v_span_index = vmid_index+1; v_span_index < v_spancnt + 1; v_span_index++) {
1682  ON_Interval u_interval(uspan[u_span_index - 1],
1683  uspan[u_span_index]);
1684  ON_Interval v_interval(vspan[v_span_index - 1],
1685  vspan[v_span_index]);
1686  double min_distance,max_distance;
1687 
1688  int level = 1;
1689  if (surface_GetIntervalMinMaxDistance(p,bbox[u_span_index-1][v_span_index-1],min_distance,max_distance)) {
1690  if ((min_distance < current_distance) && NEAR_ZERO(min_distance,within_distance_tol)) {
1691  /////////////////////////////////////////
1692  // Could check normals and CV angles here
1693  /////////////////////////////////////////
1694  double distance = surface_GetClosestPoint3dFirstOrderSubdivision(surf,p,u_interval,u_interval.Mid(),v_interval,v_interval.Mid(),current_distance,p2d,p3d,same_point_tol,within_distance_tol,level++);
1695  if (distance < current_distance) {
1696  current_distance = distance;
1697  if (current_distance < same_point_tol) {
1698  rc = true;
1699  goto cleanup;
1700  }
1701  }
1702  }
1703  }
1704  }
1705  }
1706  }
1707  if (current_distance < within_distance_tol) {
1708  rc = true;
1709  goto cleanup;
1710  }
1711  }
1712  }
1713 cleanup:
1714  std::cerr.precision(prec);
1715  return rc;
1716 }
1717 
1718 
1720  const ON_BrepTrim& trim,
1721  const ON_3dPoint& p,
1722  ON_2dPoint& p2d,
1723  double& t,
1724  double& distance,
1725  const ON_Interval* interval,
1726  double same_point_tol,
1727  double within_distance_tol
1728  )
1729 {
1730  bool rc = false;
1731  const ON_Surface *surf = trim.SurfaceOf();
1732 
1733  double t0 = interval->Mid();
1734  ON_3dPoint p3d;
1735  ON_3dPoint p0;
1736  ON_3dVector ds,dt,dss,dst,dtt;
1737  ON_3dVector T,K;
1738  int prec = std::cerr.precision();
1739  ON_BoundingBox tight_bbox;
1740  std::vector<ON_BoundingBox> bbox;
1741  std::cerr.precision(15);
1742 
1743  ON_Curve *c = trim.Brep()->m_C2[trim.m_c2i];
1744  ON_NurbsCurve N;
1745  if ( 0 == c->GetNurbForm(N) )
1746  return false;
1747  if ( N.m_order < 2 || N.m_cv_count < N.m_order )
1748  return false;
1749 
1750  p2d = trim.PointAt(t);
1751  int quadrant = 0; // optional - determines which quadrant to evaluate from
1752  // 0 = default
1753  // 1 from NE quadrant
1754  // 2 from NW quadrant
1755  // 3 from SW quadrant
1756  // 4 from SE quadrant
1757  ON_Interval u_interval = surf->Domain(0);
1758  ON_Interval v_interval = surf->Domain(1);
1759  if (p2d.y > v_interval.Mid()) {
1760  // North quadrants -> 1 or 2;
1761  if (p2d.x > u_interval.Mid()) {
1762  quadrant = 1; // NE
1763  } else {
1764  quadrant = 2; //NW
1765  }
1766  } else {
1767  // South quadrants -> 3 or 4;
1768  if (p2d.x > u_interval.Mid()) {
1769  quadrant = 4; // SE
1770  } else {
1771  quadrant = 3; //SW
1772  }
1773  }
1774  if (surface_GetClosestPoint3dFirstOrder(surf,p,p2d,p3d,distance,quadrant,same_point_tol,within_distance_tol)) {
1775  ON_BezierCurve B;
1776  bool bGrowBox = false;
1777  ON_3dVector d1,d2;
1778  double max_dist_to_closest_pt = DBL_MAX;
1779  ON_Interval *span_interval = new ON_Interval[N.m_cv_count - N.m_order + 1];
1780  double *min_distance = new double[N.m_cv_count - N.m_order + 1];
1781  double *max_distance = new double[N.m_cv_count - N.m_order + 1];
1782  bool *skip = new bool[N.m_cv_count - N.m_order + 1];
1783  bbox.resize(N.m_cv_count - N.m_order + 1);
1784  for ( int span_index = 0; span_index <= N.m_cv_count - N.m_order; span_index++ )
1785  {
1786  skip[span_index] = true;
1787  if ( !(N.m_knot[span_index + N.m_order-2] < N.m_knot[span_index + N.m_order-1]) )
1788  continue;
1789 
1790  // check for span out of interval
1791  int i = (interval->m_t[0] <= interval->m_t[1]) ? 0 : 1;
1792  if ( N.m_knot[span_index + N.m_order-2] > interval->m_t[1-i] )
1793  continue;
1794  if ( N.m_knot[span_index + N.m_order-1] < interval->m_t[i] )
1795  continue;
1796 
1797  if ( !N.ConvertSpanToBezier( span_index, B ) )
1798  continue;
1799  ON_Interval bi = B.Domain();
1800  if ( !B.GetTightBoundingBox(tight_bbox,bGrowBox,NULL) )
1801  continue;
1802  bbox[span_index] = tight_bbox;
1803  d1 = tight_bbox.m_min - p2d;
1804  d2 = tight_bbox.m_max - p2d;
1805  min_distance[span_index] = tight_bbox.MinimumDistanceTo(p2d);
1806 
1807  if (min_distance[span_index] > max_dist_to_closest_pt) {
1808  max_distance[span_index] = DBL_MAX;
1809  continue;
1810  }
1811  skip[span_index] = false;
1812  span_interval[span_index].m_t[0] = ((N.m_knot[span_index + N.m_order-2]) < interval->m_t[i]) ? interval->m_t[i] : N.m_knot[span_index + N.m_order-2];
1813  span_interval[span_index].m_t[1] = ((N.m_knot[span_index + N.m_order-1]) > interval->m_t[1 -i]) ? interval->m_t[1 -i] : (N.m_knot[span_index + N.m_order-1]);
1814  ON_3dPoint d1sq(d1.x*d1.x,d1.y*d1.y,0.0),d2sq(d2.x*d2.x,d2.y*d2.y,0.0);
1815  double distancesq;
1816  if (d1sq.x < d2sq.x) {
1817  if (d1sq.y < d2sq.y) {
1818  if ((d1sq.x + d2sq.y) < (d2sq.x + d1sq.y)) {
1819  distancesq = d1sq.x + d2sq.y;
1820  } else {
1821  distancesq = d2sq.x + d1sq.y;
1822  }
1823  } else {
1824  if ((d1sq.x + d1sq.y) < (d2sq.x + d2sq.y)) {
1825  distancesq = d1sq.x + d1sq.y;
1826  } else {
1827  distancesq = d2sq.x + d2sq.y;
1828  }
1829  }
1830  } else {
1831  if (d1sq.y < d2sq.y) {
1832  if ((d1sq.x + d1sq.y) < (d2sq.x + d2sq.y)) {
1833  distancesq = d1sq.x + d1sq.y;
1834  } else {
1835  distancesq = d2sq.x + d2sq.y;
1836  }
1837  } else {
1838  if ((d1sq.x + d2sq.y) < (d2sq.x + d1sq.y)) {
1839  distancesq = d1sq.x + d2sq.y;
1840  } else {
1841  distancesq = d2sq.x + d1sq.y;
1842  }
1843  }
1844  }
1845  max_distance[span_index] = sqrt(distancesq);
1846  if (max_distance[span_index] < max_dist_to_closest_pt) {
1847  max_dist_to_closest_pt = max_distance[span_index];
1848  }
1849  if (max_distance[span_index] < min_distance[span_index]) {
1850  // should only be here for near equal fuzz
1851  min_distance[span_index] = max_distance[span_index];
1852  }
1853  }
1854  for ( int span_index = 0; span_index <= N.m_cv_count - N.m_order; span_index++ )
1855  {
1856 
1857  if ( skip[span_index] )
1858  continue;
1859 
1860  if (min_distance[span_index] > max_dist_to_closest_pt) {
1861  skip[span_index] = true;
1862  continue;
1863  }
1864 
1865  }
1866 
1867  ON_3dPoint q;
1868  ON_3dPoint point;
1869  double closest_distance = DBL_MAX;
1870  double closestT = DBL_MAX;
1871  for ( int span_index = 0; span_index <= N.m_cv_count - N.m_order; span_index++ )
1872  {
1873  if (skip[span_index]) {
1874  continue;
1875  }
1876  t0 = span_interval[span_index].Mid();
1877  bool closestfound = false;
1878  bool notdone = true;
1879  double span_distance = DBL_MAX;
1880  double previous_distance = DBL_MAX;
1881  ON_3dVector firstDervative, secondDervative;
1882  while (notdone
1883  && trim.Ev2Der(t0, point, firstDervative, secondDervative)
1884  && ON_EvCurvature(firstDervative, secondDervative, T, K)) {
1885  ON_Line line(point, point + 100.0 * T);
1886  q = line.ClosestPointTo(p2d);
1887  double delta_t = (firstDervative * (q - point))
1888  / (firstDervative * firstDervative);
1889  double new_t0 = t0 + delta_t;
1890  if (!span_interval[span_index].Includes(new_t0, false)) {
1891  // limit to interval
1892  int i = (span_interval[span_index].m_t[0] <= span_interval[span_index].m_t[1]) ? 0 : 1;
1893  new_t0 =
1894  (new_t0 < span_interval[span_index].m_t[i]) ?
1895  span_interval[span_index].m_t[i] : span_interval[span_index].m_t[1 - i];
1896  }
1897  delta_t = new_t0 - t0;
1898  t0 = new_t0;
1899  point = trim.PointAt(t0);
1900  span_distance = point.DistanceTo(p2d);
1901  if (span_distance < previous_distance) {
1902  closestfound = true;
1903  closestT = t0;
1904  previous_distance = span_distance;
1905  if (fabs(delta_t) < same_point_tol) {
1906  notdone = false;
1907  }
1908  } else {
1909  notdone = false;
1910  }
1911  }
1912  if (closestfound && (span_distance < closest_distance)) {
1913  closest_distance = span_distance;
1914  rc = true;
1915  t = closestT;
1916  }
1917  }
1918  delete [] span_interval;
1919  delete [] min_distance;
1920  delete [] max_distance;
1921  delete [] skip;
1922 
1923  }
1924  std::cerr.precision(prec);
1925 
1926  return rc;
1927 }
1928 
1929 
1930 bool
1931 toUV(PBCData& data, ON_2dPoint& out_pt, double t, double knudge = 0.0, double within_distance_tol = BREP_EDGE_MISS_TOLERANCE)
1932 {
1933  ON_3dPoint pointOnCurve = data.curve->PointAt(t);
1934  ON_3dPoint knudgedPointOnCurve = data.curve->PointAt(t + knudge);
1935 
1936  ON_2dPoint uv;
1937  if (data.surftree->getSurfacePoint((const ON_3dPoint&)pointOnCurve, uv, (const ON_3dPoint&)knudgedPointOnCurve, within_distance_tol) > 0) {
1938  out_pt.Set(uv.x, uv.y);
1939  return true;
1940  } else {
1941  return false;
1942  }
1943 }
1944 
1945 
1946 bool
1947 toUV(brlcad::SurfaceTree *surftree, const ON_Curve *curve, ON_2dPoint& out_pt, double t, double knudge = 0.0)
1948 {
1949  if (!surftree)
1950  return false;
1951 
1952  const ON_Surface *surf = surftree->getSurface();
1953  if (!surf)
1954  return false;
1955 
1956  ON_3dPoint pointOnCurve = curve->PointAt(t);
1957  ON_3dPoint knudgedPointOnCurve = curve->PointAt(t + knudge);
1958  ON_3dVector dt;
1959  curve->Ev1Der(t, pointOnCurve, dt);
1960  ON_3dVector tangent = curve->TangentAt(t);
1961  //data.surf->GetClosestPoint(pointOnCurve, &a, &b, 0.0001);
1962  ON_Ray r(pointOnCurve, tangent);
1963  plane_ray pr;
1964  brep_get_plane_ray(r, pr);
1965  ON_3dVector p1;
1966  double p1d;
1967  ON_3dVector p2;
1968  double p2d;
1969 
1970  utah_ray_planes(r, p1, p1d, p2, p2d);
1971 
1972  VMOVE(pr.n1, p1);
1973  pr.d1 = p1d;
1974  VMOVE(pr.n2, p2);
1975  pr.d2 = p2d;
1976 
1977  try {
1978  pt2d_t uv;
1979  ON_2dPoint uv2d = surftree->getClosestPointEstimate(knudgedPointOnCurve);
1980  move(uv, uv2d);
1981 
1982  ON_3dVector dir = surf->NormalAt(uv[0], uv[1]);
1983  dir.Reverse();
1984  ON_Ray ray(pointOnCurve, dir);
1985  brep_get_plane_ray(ray, pr);
1986  //know use this as guess to iterate to closer solution
1987  pt2d_t Rcurr;
1988  pt2d_t new_uv;
1989  ON_3dPoint pt;
1990  ON_3dVector su, sv;
1991 #ifdef SHOW_UNUSED
1992  bool found = false;
1993 #endif
1994  fastf_t Dlast = MAX_FASTF;
1995  for (int i = 0; i < 10; i++) {
1996  brep_r(surf, pr, uv, pt, su, sv, Rcurr);
1997  fastf_t d = v2mag(Rcurr);
1999  TRACE1("R:" << ON_PRINT2(Rcurr));
2000 #ifdef SHOW_UNUSED
2001  found = true;
2002  break;
2003 #endif
2004  } else if (d > Dlast) {
2005 #ifdef SHOW_UNUSED
2006  found = false; //break;
2007 #endif
2008  break;
2009  //return brep_edge_check(found, sbv, face, surf, ray, hits);
2010  }
2011  brep_newton_iterate(pr, Rcurr, su, sv, uv, new_uv);
2012  move(uv, new_uv);
2013  Dlast = d;
2014  }
2015 
2016 ///////////////////////////////////////
2017  out_pt.Set(uv[0], uv[1]);
2018  return true;
2019  } catch (...) {
2020  return false;
2021  }
2022 }
2023 
2024 
2025 double
2026 randomPointFromRange(PBCData& data, ON_2dPoint& out, double lo, double hi)
2027 {
2028  assert(lo < hi);
2029  double random_pos = drand48() * (RANGE_HI - RANGE_LO) + RANGE_LO;
2030  double newt = random_pos * (hi - lo) + lo;
2031  assert(toUV(data, out, newt));
2032  return newt;
2033 }
2034 
2035 
2036 bool
2037 sample(PBCData& data,
2038  double t1,
2039  double t2,
2040  const ON_2dPoint& p1,
2041  const ON_2dPoint& p2)
2042 {
2043  ON_2dPoint m;
2044  double t = randomPointFromRange(data, m, t1, t2);
2045  if (!data.segments.empty()) {
2046  ON_2dPointArray * samples = data.segments.back();
2047  if (isFlat(p1, m, p2, data.flatness)) {
2048  samples->Append(p2);
2049  } else {
2050  sample(data, t1, t, p1, m);
2051  sample(data, t, t2, m, p2);
2052  }
2053  return true;
2054  }
2055  return false;
2056 }
2057 
2058 
2059 void
2061 {
2062  int num_knots = bspline.m + 1;
2063  bspline.knots.resize(num_knots);
2064  for (int i = 0; i <= bspline.p; i++) {
2065  bspline.knots[i] = 0.0;
2066  }
2067  for (int i = bspline.m - bspline.p; i <= bspline.m; i++) {
2068  bspline.knots[i] = 1.0;
2069  }
2070  for (int i = 1; i <= bspline.n - bspline.p; i++) {
2071  bspline.knots[bspline.p + i] = (double)i / (bspline.n - bspline.p + 1.0);
2072  }
2073 }
2074 
2075 
2076 int
2077 getKnotInterval(BSpline& bspline, double u)
2078 {
2079  int k = 0;
2080  while (u >= bspline.knots[k]) k++;
2081  k = (k == 0) ? k : k - 1;
2082  return k;
2083 }
2084 
2085 ON_NurbsCurve*
2087 {
2088  int num_samples = Q.Count();
2089  int num_segments = Q.Count() - 1;
2090  int qsize = num_samples + 4;
2091  std::vector < ON_2dVector > qarray(qsize);
2092 
2093  for (int i = 1; i < Q.Count(); i++) {
2094  qarray[i + 1] = Q[i] - Q[i - 1];
2095  }
2096  qarray[1] = 2.0 * qarray[2] - qarray[3];
2097  qarray[0] = 2.0 * qarray[1] - qarray[2];
2098 
2099  qarray[num_samples + 1] = 2 * qarray[num_samples] - qarray[num_samples - 1];
2100  qarray[num_samples + 2] = 2 * qarray[num_samples + 1] - qarray[num_samples];
2101  qarray[num_samples + 3] = 2 * qarray[num_samples + 2] - qarray[num_samples + 1];
2102 
2103  std::vector < ON_2dVector > T(num_samples);
2104  std::vector<double> A(num_samples);
2105  for (int k = 0; k < num_samples; k++) {
2106  ON_3dVector a = ON_CrossProduct(qarray[k], qarray[k + 1]);
2107  ON_3dVector b = ON_CrossProduct(qarray[k + 2], qarray[k + 3]);
2108  double alength = a.Length();
2109  if (NEAR_ZERO(alength, PBC_TOL)) {
2110  A[k] = 1.0;
2111  } else {
2112  A[k] = (a.Length()) / (a.Length() + b.Length());
2113  }
2114  T[k] = (1.0 - A[k]) * qarray[k + 1] + A[k] * qarray[k + 2];
2115  T[k].Unitize();
2116  }
2117  std::vector < ON_2dPointArray > P(num_samples - 1);
2118  ON_2dPointArray control_points;
2119  control_points.Append(Q[0]);
2120  for (int i = 1; i < num_samples; i++) {
2121  ON_2dPoint P0 = Q[i - 1];
2122  ON_2dPoint P3 = Q[i];
2123  ON_2dVector T0 = T[i - 1];
2124  ON_2dVector T3 = T[i];
2125 
2126  double a, b, c;
2127 
2128  ON_2dVector vT0T3 = T0 + T3;
2129  ON_2dVector dP0P3 = P3 - P0;
2130  a = 16.0 - vT0T3.Length() * vT0T3.Length();
2131  b = 12.0 * (dP0P3 * vT0T3);
2132  c = -36.0 * dP0P3.Length() * dP0P3.Length();
2133 
2134  double alpha = (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
2135 
2136  ON_2dPoint P1 = P0 + (1.0 / 3.0) * alpha * T0;
2137  control_points.Append(P1);
2138  ON_2dPoint P2 = P3 - (1.0 / 3.0) * alpha * T3;
2139  control_points.Append(P2);
2140  P[i - 1].Append(P0);
2141  P[i - 1].Append(P1);
2142  P[i - 1].Append(P2);
2143  P[i - 1].Append(P3);
2144  }
2145  control_points.Append(Q[num_samples - 1]);
2146 
2147  std::vector<double> u(num_segments + 1);
2148  u[0] = 0.0;
2149  for (int k = 0; k < num_segments; k++) {
2150  u[k + 1] = u[k] + 3.0 * (P[k][1] - P[k][0]).Length();
2151  }
2152  int degree = 3;
2153  int n = control_points.Count();
2154  int p = degree;
2155  int m = n + p - 1;
2156  int dimension = 2;
2157  ON_NurbsCurve* c = ON_NurbsCurve::New(dimension, false, degree + 1, n);
2158  c->ReserveKnotCapacity(m);
2159  for (int i = 0; i < degree; i++) {
2160  c->SetKnot(i, 0.0);
2161  }
2162  for (int i = 1; i < num_segments; i++) {
2163  double knot_value = u[i] / u[num_segments];
2164  c->SetKnot(degree + 2 * (i - 1), knot_value);
2165  c->SetKnot(degree + 2 * (i - 1) + 1, knot_value);
2166  }
2167  for (int i = m - p; i < m; i++) {
2168  c->SetKnot(i, 1.0);
2169  }
2170 
2171  // insert the control points
2172  for (int i = 0; i < n; i++) {
2173  ON_3dPoint pnt = control_points[i];
2174  c->SetCV(i, pnt);
2175  }
2176  return c;
2177 }
2178 
2179 ON_NurbsCurve*
2180 interpolateLocalCubicCurve(const ON_3dPointArray &Q)
2181 {
2182  int num_samples = Q.Count();
2183  int num_segments = Q.Count() - 1;
2184  int qsize = num_samples + 3;
2185  std::vector<ON_3dVector> qarray(qsize + 1);
2186  ON_3dVector *q = &qarray[1];
2187 
2188  for (int i = 1; i < Q.Count(); i++) {
2189  q[i] = Q[i] - Q[i - 1];
2190  }
2191  q[0] = 2.0 * q[1] - q[2];
2192  q[-1] = 2.0 * q[0] - q[1];
2193 
2194  q[num_samples] = 2 * q[num_samples - 1] - q[num_samples - 2];
2195  q[num_samples + 1] = 2 * q[num_samples] - q[num_samples - 1];
2196  q[num_samples + 2] = 2 * q[num_samples + 1] - q[num_samples];
2197 
2198  std::vector<ON_3dVector> T(num_samples);
2199  std::vector<double> A(num_samples);
2200  for (int k = 0; k < num_samples; k++) {
2201  ON_3dVector avec = ON_CrossProduct(q[k - 1], q[k]);
2202  ON_3dVector bvec = ON_CrossProduct(q[k + 1], q[k + 2]);
2203  double alength = avec.Length();
2204  if (NEAR_ZERO(alength, PBC_TOL)) {
2205  A[k] = 1.0;
2206  } else {
2207  A[k] = (avec.Length()) / (avec.Length() + bvec.Length());
2208  }
2209  T[k] = (1.0 - A[k]) * q[k] + A[k] * q[k + 1];
2210  T[k].Unitize();
2211  }
2212  std::vector<ON_3dPointArray> P(num_samples - 1);
2213  ON_3dPointArray control_points;
2214  control_points.Append(Q[0]);
2215  for (int i = 1; i < num_samples; i++) {
2216  ON_3dPoint P0 = Q[i - 1];
2217  ON_3dPoint P3 = Q[i];
2218  ON_3dVector T0 = T[i - 1];
2219  ON_3dVector T3 = T[i];
2220 
2221  double a, b, c;
2222 
2223  ON_3dVector vT0T3 = T0 + T3;
2224  ON_3dVector dP0P3 = P3 - P0;
2225  a = 16.0 - vT0T3.Length() * vT0T3.Length();
2226  b = 12.0 * (dP0P3 * vT0T3);
2227  c = -36.0 * dP0P3.Length() * dP0P3.Length();
2228 
2229  double alpha = (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
2230 
2231  ON_3dPoint P1 = P0 + (1.0 / 3.0) * alpha * T0;
2232  control_points.Append(P1);
2233  ON_3dPoint P2 = P3 - (1.0 / 3.0) * alpha * T3;
2234  control_points.Append(P2);
2235  P[i - 1].Append(P0);
2236  P[i - 1].Append(P1);
2237  P[i - 1].Append(P2);
2238  P[i - 1].Append(P3);
2239  }
2240  control_points.Append(Q[num_samples - 1]);
2241 
2242  std::vector<double> u(num_segments + 1);
2243  u[0] = 0.0;
2244  for (int k = 0; k < num_segments; k++) {
2245  u[k + 1] = u[k] + 3.0 * (P[k][1] - P[k][0]).Length();
2246  }
2247  int degree = 3;
2248  int n = control_points.Count();
2249  int p = degree;
2250  int m = n + p - 1;
2251  int dimension = 3;
2252  ON_NurbsCurve* c = ON_NurbsCurve::New(dimension, false, degree + 1, n);
2253  c->ReserveKnotCapacity(m);
2254  for (int i = 0; i < degree; i++) {
2255  c->SetKnot(i, 0.0);
2256  }
2257  for (int i = 1; i < num_segments; i++) {
2258  double knot_value = u[i] / u[num_segments];
2259  c->SetKnot(degree + 2 * (i - 1), knot_value);
2260  c->SetKnot(degree + 2 * (i - 1) + 1, knot_value);
2261  }
2262  for (int i = m - p; i < m; i++) {
2263  c->SetKnot(i, 1.0);
2264  }
2265 
2266  // insert the control points
2267  for (int i = 0; i < n; i++) {
2268  ON_3dPoint pnt = control_points[i];
2269  c->SetCV(i, pnt);
2270  }
2271  return c;
2272 }
2273 
2274 
2275 ON_NurbsCurve*
2276 newNURBSCurve(BSpline& spline, int dimension = 3)
2277 {
2278  // we now have everything to complete our spline
2279  ON_NurbsCurve* c = ON_NurbsCurve::New(dimension,
2280  false,
2281  spline.p + 1,
2282  spline.n + 1);
2283  c->ReserveKnotCapacity(spline.knots.size() - 2);
2284  for (unsigned int i = 1; i < spline.knots.size() - 1; i++) {
2285  c->m_knot[i - 1] = spline.knots[i];
2286  }
2287 
2288  for (int i = 0; i < spline.controls.Count(); i++) {
2289  c->SetCV(i, ON_3dPoint(spline.controls[i]));
2290  }
2291 
2292  return c;
2293 }
2294 
2295 
2296 ON_Curve*
2297 interpolateCurve(ON_2dPointArray &samples)
2298 {
2299  ON_NurbsCurve* nurbs;
2300 
2301  if (samples.Count() == 2)
2302  // build a line
2303  return new ON_LineCurve(samples[0], samples[1]);
2304 
2305  // local vs. global interpolation for large point sampled curves
2306  nurbs = interpolateLocalCubicCurve(samples);
2307 
2308  return nurbs;
2309 }
2310 
2311 
2312 int
2313 IsAtSeam(const ON_Surface *surf, int dir, double u, double v, double tol)
2314 {
2315  int rc = 0;
2316  if (!surf->IsClosed(dir))
2317  return rc;
2318 
2319  double p = (dir) ? v : u;
2320  if (NEAR_EQUAL(p, surf->Domain(dir)[0], tol) || NEAR_EQUAL(p, surf->Domain(dir)[1], tol))
2321  rc += (dir + 1);
2322 
2323  return rc;
2324 }
2325 
2326 /*
2327  * Similar to openNURBS's surf->IsAtSeam() function but uses tolerance to do a near check versus
2328  * the floating point equality used by openNURBS.
2329  * rc = 0 Not on seam, 1 on East/West seam(umin/umax), 2 on North/South seam(vmin/vmax), 3 seam on both U/V boundaries
2330  */
2331 int
2332 IsAtSeam(const ON_Surface *surf, double u, double v, double tol)
2333 {
2334  int rc = 0;
2335  int i;
2336  for (i = 0; i < 2; i++) {
2337  rc += IsAtSeam(surf, i, u, v, tol);
2338  }
2339 
2340  return rc;
2341 }
2342 
2343 int
2344 IsAtSeam(const ON_Surface *surf, int dir, const ON_2dPoint &pt, double tol)
2345 {
2346  int rc = 0;
2347  ON_2dPoint unwrapped_pt = UnwrapUVPoint(surf,pt,tol);
2348  rc = IsAtSeam(surf,dir,unwrapped_pt.x,unwrapped_pt.y,tol);
2349 
2350  return rc;
2351 }
2352 
2353 /*
2354  * Similar to IsAtSeam(surf,u,v,tol) function but takes a ON_2dPoint
2355  * and unwraps any closed seam extents before passing on IsAtSeam(surf,u,v,tol)
2356  */
2357 int
2358 IsAtSeam(const ON_Surface *surf, const ON_2dPoint &pt, double tol)
2359 {
2360  int rc = 0;
2361  ON_2dPoint unwrapped_pt = UnwrapUVPoint(surf,pt,tol);
2362  rc = IsAtSeam(surf,unwrapped_pt.x,unwrapped_pt.y,tol);
2363 
2364  return rc;
2365 }
2366 
2367 
2368 int
2369 IsAtSingularity(const ON_Surface *surf, double u, double v, double tol)
2370 {
2371  // 0 = south, 1 = east, 2 = north, 3 = west
2372  //std::cerr << "IsAtSingularity = u, v - " << u << ", " << v << std::endl;
2373  //std::cerr << "surf->Domain(0) - " << surf->Domain(0)[0] << ", " << surf->Domain(0)[1] << std::endl;
2374  //std::cerr << "surf->Domain(1) - " << surf->Domain(1)[0] << ", " << surf->Domain(1)[1] << std::endl;
2375  if (NEAR_EQUAL(u, surf->Domain(0)[0], tol)) {
2376  if (surf->IsSingular(3))
2377  return 3;
2378  } else if (NEAR_EQUAL(u, surf->Domain(0)[1], tol)) {
2379  if (surf->IsSingular(1))
2380  return 1;
2381  }
2382  if (NEAR_EQUAL(v, surf->Domain(1)[0], tol)) {
2383  if (surf->IsSingular(0))
2384  return 0;
2385  } else if (NEAR_EQUAL(v, surf->Domain(1)[1], tol)) {
2386  if (surf->IsSingular(2))
2387  return 2;
2388  }
2389  return -1;
2390 }
2391 
2392 int
2393 IsAtSingularity(const ON_Surface *surf, const ON_2dPoint &pt, double tol)
2394 {
2395  int rc = 0;
2396  ON_2dPoint unwrapped_pt = UnwrapUVPoint(surf,pt,tol);
2397  rc = IsAtSingularity(surf,unwrapped_pt.x,unwrapped_pt.y,tol);
2398 
2399  return rc;
2400 }
2401 
2402 ON_2dPointArray *
2404  double t,
2405  double s,
2406  double same_point_tol,
2407  double within_distance_tol)
2408 {
2409  if (!data)
2410  return NULL;
2411 
2412  const ON_Curve* curve = data->curve;
2413  const ON_Surface* surf = data->surf;
2414  ON_2dPointArray *samples = new ON_2dPointArray();
2415  int numKnots = curve->SpanCount();
2416  double *knots = new double[numKnots + 1];
2417  curve->GetSpanVector(knots);
2418 
2419  int istart = 0;
2420  while (t >= knots[istart])
2421  istart++;
2422 
2423  if (istart > 0) {
2424  istart--;
2425  knots[istart] = t;
2426  }
2427 
2428  int istop = numKnots;
2429  while (s <= knots[istop])
2430  istop--;
2431 
2432  if (istop < numKnots) {
2433  istop++;
2434  knots[istop] = s;
2435  }
2436 
2437  int samplesperknotinterval;
2438  int degree = curve->Degree();
2439 
2440  if (degree > 1) {
2441  samplesperknotinterval = 3 * degree;
2442  } else {
2443  samplesperknotinterval = 18 * degree;
2444  }
2445  ON_2dPoint pt;
2446  ON_3dPoint p = ON_3dPoint::UnsetPoint;
2447  ON_3dPoint p3d = ON_3dPoint::UnsetPoint;
2448  double distance;
2449  for (int i = istart; i <= istop; i++) {
2450  if (i <= numKnots / 2) {
2451  if (i > 0) {
2452  double delta = (knots[i] - knots[i - 1]) / (double) samplesperknotinterval;
2453  for (int j = 1; j < samplesperknotinterval; j++) {
2454  p = curve->PointAt(knots[i - 1] + j * delta);
2455  p3d = ON_3dPoint::UnsetPoint;
2456  if (surface_GetClosestPoint3dFirstOrder(surf,p,pt,p3d,distance,0,same_point_tol,within_distance_tol)) {
2457  samples->Append(pt);
2458  }
2459  }
2460  }
2461  p = curve->PointAt(knots[i]);
2462  p3d = ON_3dPoint::UnsetPoint;
2463  if (surface_GetClosestPoint3dFirstOrder(surf,p,pt,p3d,distance,0,same_point_tol,within_distance_tol)) {
2464  samples->Append(pt);
2465  }
2466  } else {
2467  if (i > 0) {
2468  double delta = (knots[i] - knots[i - 1]) / (double) samplesperknotinterval;
2469  for (int j = 1; j < samplesperknotinterval; j++) {
2470  p = curve->PointAt(knots[i - 1] + j * delta);
2471  p3d = ON_3dPoint::UnsetPoint;
2472  if (surface_GetClosestPoint3dFirstOrder(surf,p,pt,p3d,distance,0,same_point_tol,within_distance_tol)) {
2473  samples->Append(pt);
2474  }
2475  }
2476  p = curve->PointAt(knots[i]);
2477  p3d = ON_3dPoint::UnsetPoint;
2478  if (surface_GetClosestPoint3dFirstOrder(surf,p,pt,p3d,distance,0,same_point_tol,within_distance_tol)) {
2479  samples->Append(pt);
2480  }
2481  }
2482  }
2483  }
2484  delete[] knots;
2485  return samples;
2486 }
2487 
2488 
2489 /*
2490  * Unwrap 2D UV point values to within actual surface UV. Points often wrap around the closed seam.
2491  */
2492 ON_2dPoint
2493 UnwrapUVPoint(const ON_Surface *surf,const ON_2dPoint &pt, double tol)
2494 {
2495  ON_2dPoint p = pt;
2496  for (int i=0; i<2; i++) {
2497  if (!surf->IsClosed(i))
2498  continue;
2499  while (p[i] < surf->Domain(i).m_t[0] - tol) {
2500  double length = surf->Domain(i).Length();
2501  if (i<=0) {
2502  p.x = p.x + length;
2503  } else {
2504  p.y = p.y + length;
2505  }
2506  }
2507  while (p[i] >= surf->Domain(i).m_t[1] + tol) {
2508  double length = surf->Domain(i).Length();
2509  if (i<=0) {
2510  p.x = p.x - length;
2511  } else {
2512  p.y = p.y - length;
2513  }
2514  }
2515  }
2516 
2517  return p;
2518 }
2519 
2520 
2521 double
2522 DistToNearestClosedSeam(const ON_Surface *surf,const ON_2dPoint &pt)
2523 {
2524  double dist = -1.0;
2525  ON_2dPoint unwrapped_pt = UnwrapUVPoint(surf,pt);
2526  for (int i=0; i<2; i++) {
2527  if (!surf->IsClosed(i))
2528  continue;
2529  dist = fabs(unwrapped_pt[i] - surf->Domain(i)[0]);
2530  V_MIN(dist,fabs(surf->Domain(i)[1]-unwrapped_pt[i]));
2531  }
2532  return dist;
2533 }
2534 
2535 
2536 void
2537 GetClosestExtendedPoint(const ON_Surface *surf,ON_2dPoint &pt,ON_2dPoint &prev_pt, double UNUSED(tol)) {
2538  if (surf->IsClosed(0)) {
2539  double length = surf->Domain(0).Length();
2540  double delta=pt.x-prev_pt.x;
2541  while (fabs(delta) > length/2.0) {
2542  if (delta > length/2.0) {
2543  pt.x = pt.x - length;
2544  delta=pt.x-prev_pt.x;
2545  } else {
2546  pt.x = pt.x + length;
2547  delta=pt.x - prev_pt.x;
2548  }
2549  }
2550  }
2551 
2552  if (surf->IsClosed(1)) {
2553  double length = surf->Domain(1).Length();
2554  double delta=pt.y-prev_pt.y;
2555  while (fabs(delta) > length/2.0) {
2556  if (delta > length/2.0) {
2557  pt.y = pt.y - length;
2558  delta=pt.y - prev_pt.y;
2559  } else {
2560  pt.y = pt.y + length;
2561  delta=pt.y -prev_pt.y;
2562  }
2563  }
2564  }
2565 }
2566 
2567 
2568 /*
2569  * Simple check to determine if two consecutive points pulled back from 3d curve sampling
2570  * to 2d UV parameter space crosses the seam of the closed UV. The assumption here is that
2571  * the sampling of the 3d curve is a small fraction of the UV domain.
2572  *
2573  * // dir - 0 = not crossing, 1 = south/east bound, 2 = north/west bound
2574  */
2575 bool
2576 ConsecutivePointsCrossClosedSeam(const ON_Surface *surf,const ON_2dPoint &pt,const ON_2dPoint &prev_pt, int &udir, int &vdir, double tol)
2577 {
2578  bool rc = false;
2579 
2580  /*
2581  * if one of the points is at a seam then not crossing
2582  */
2583  int dir =0;
2584  ON_2dPoint unwrapped_pt = UnwrapUVPoint(surf,pt);
2585  ON_2dPoint unwrapped_prev_pt = UnwrapUVPoint(surf,prev_pt);
2586 
2587  if (!IsAtSeam(surf,dir,pt,tol) && !IsAtSeam(surf,dir,prev_pt,tol)) {
2588  udir = vdir = 0;
2589  if (surf->IsClosed(0)) {
2590  double delta=unwrapped_pt.x-unwrapped_prev_pt.x;
2591  if (fabs(delta) > surf->Domain(0).Length()/2.0) {
2592  if (delta < 0.0) {
2593  udir = 1; // east bound
2594  } else {
2595  udir= 2; // west bound
2596  }
2597  rc = true;
2598  }
2599  }
2600  }
2601  dir = 1;
2602  if (!IsAtSeam(surf,dir,pt,tol) && !IsAtSeam(surf,dir,prev_pt,tol)) {
2603  if (surf->IsClosed(1)) {
2604  double delta=unwrapped_pt.y-unwrapped_prev_pt.y;
2605  if (fabs(delta) > surf->Domain(1).Length()/2.0) {
2606  if (delta < 0.0) {
2607  vdir = 2; // north bound
2608  } else {
2609  vdir= 1; // south bound
2610  }
2611  rc = true;
2612  }
2613  }
2614  }
2615 
2616  return rc;
2617 }
2618 
2619 /*
2620  * If within UV tolerance to a seam force force to actually seam value so surface
2621  * seam function can be used.
2622  */
2623 void
2624 ForceToClosestSeam(const ON_Surface *surf, ON_2dPoint &pt, double tol)
2625 {
2626  int seam;
2627  ON_2dPoint unwrapped_pt = UnwrapUVPoint(surf,pt,tol);
2628  ON_2dVector wrap = ON_2dVector::ZeroVector;
2629  ON_Interval dom[2] = { ON_Interval::EmptyInterval, ON_Interval::EmptyInterval };
2630  double length[2] = { ON_UNSET_VALUE, ON_UNSET_VALUE};
2631 
2632  for (int i=0; i<2; i++) {
2633  dom[i] = surf->Domain(i);
2634  length[i] = dom[i].Length();
2635  if (!surf->IsClosed(i))
2636  continue;
2637  if (pt[i] > (dom[i].m_t[1] + tol)) {
2638  ON_2dPoint p = pt;
2639  while (p[i] > (dom[i].m_t[1] + tol)) {
2640  p[i] -= length[i];
2641  wrap[i] += 1.0;
2642  }
2643  } else if (pt[i] < (dom[i].m_t[0] - tol)) {
2644  ON_2dPoint p = pt;
2645  while (p[i] < (dom[i].m_t[0] - tol)) {
2646  p[i] += length[i];
2647  wrap[i] -= 1.0;
2648  }
2649  }
2650  wrap[i] = wrap[i] * length[i];
2651  }
2652 
2653  if ((seam=IsAtSeam(surf, unwrapped_pt, tol)) > 0) {
2654  if (seam == 1) { // east/west seam
2655  if (fabs(unwrapped_pt.x - dom[0].m_t[0]) < length[0]/2.0) {
2656  unwrapped_pt.x = dom[0].m_t[0]; // on east swap to west seam
2657  } else {
2658  unwrapped_pt.x = dom[0].m_t[1]; // on west swap to east seam
2659  }
2660  } else if (seam == 2) { // north/south seam
2661  if (fabs(unwrapped_pt.y - dom[1].m_t[0]) < length[1]/2.0) {
2662  unwrapped_pt.y = dom[1].m_t[0]; // on north swap to south seam
2663  } else {
2664  unwrapped_pt.y = dom[1].m_t[1]; // on south swap to north seam
2665  }
2666  } else { //on both seams
2667  if (fabs(unwrapped_pt.x - dom[0].m_t[0]) < length[0]/2.0) {
2668  unwrapped_pt.x = dom[0].m_t[0]; // on east swap to west seam
2669  } else {
2670  unwrapped_pt.x = dom[0].m_t[1]; // on west swap to east seam
2671  }
2672  if (fabs(pt.y - dom[1].m_t[0]) < length[1]/2.0) {
2673  unwrapped_pt.y = dom[1].m_t[0]; // on north swap to south seam
2674  } else {
2675  unwrapped_pt.y = dom[1].m_t[1]; // on south swap to north seam
2676  }
2677  }
2678  }
2679  pt = unwrapped_pt + wrap;
2680 }
2681 
2682 
2683 /*
2684  * If point lies on a seam(s) swap to opposite side of UV.
2685  * hint = 1 swap E/W, 2 swap N/S or default 3 swap both
2686  */
2687 void
2688 SwapUVSeamPoint(const ON_Surface *surf, ON_2dPoint &p, int hint)
2689 {
2690  int seam;
2691  ON_Interval dom[2];
2692  dom[0] = surf->Domain(0);
2693  dom[1] = surf->Domain(1);
2694 
2695  if ((seam=surf->IsAtSeam(p.x,p.y)) > 0) {
2696  if (seam == 1) { // east/west seam
2697  if (fabs(p.x - dom[0].m_t[0]) > dom[0].Length()/2.0) {
2698  p.x = dom[0].m_t[0]; // on east swap to west seam
2699  } else {
2700  p.x = dom[0].m_t[1]; // on west swap to east seam
2701  }
2702  } else if (seam == 2) { // north/south seam
2703  if (fabs(p.y - dom[1].m_t[0]) > dom[1].Length()/2.0) {
2704  p.y = dom[1].m_t[0]; // on north swap to south seam
2705  } else {
2706  p.y = dom[1].m_t[1]; // on south swap to north seam
2707  }
2708  } else { //on both seams check hint 1=east/west only, 2=north/south, 3 = both
2709  if (hint == 1) {
2710  if (fabs(p.x - dom[0].m_t[0]) > dom[0].Length()/2.0) {
2711  p.x = dom[0].m_t[0]; // on east swap to west seam
2712  } else {
2713  p.x = dom[0].m_t[1]; // on west swap to east seam
2714  }
2715  } else if (hint == 2) {
2716  if (fabs(p.y - dom[1].m_t[0]) > dom[1].Length()/2.0) {
2717  p.y = dom[1].m_t[0]; // on north swap to south seam
2718  } else {
2719  p.y = dom[1].m_t[1]; // on south swap to north seam
2720  }
2721  } else {
2722  if (fabs(p.x - dom[0].m_t[0]) > dom[0].Length()/2.0) {
2723  p.x = dom[0].m_t[0]; // on east swap to west seam
2724  } else {
2725  p.x = dom[0].m_t[1]; // on west swap to east seam
2726  }
2727  if (fabs(p.y - dom[1].m_t[0]) > dom[1].Length()/2.0) {
2728  p.y = dom[1].m_t[0]; // on north swap to south seam
2729  } else {
2730  p.y = dom[1].m_t[1]; // on south swap to north seam
2731  }
2732  }
2733  }
2734  }
2735 }
2736 
2737 
2738 /*
2739  * Find where Pullback of 3d curve crosses closed seam of surface UV
2740  */
2741 bool
2742 Find3DCurveSeamCrossing(PBCData &data,double t0,double t1, double UNUSED(offset),double &seam_t,ON_2dPoint &from,ON_2dPoint &to,double tol,double same_point_tol,double within_distance_tol)
2743 {
2744  bool rc = true;
2745  const ON_Surface *surf = data.surf;
2746 
2747  // quick bail out is surface not closed
2748  if (surf->IsClosed(0) || surf->IsClosed(1)) {
2749  ON_2dPoint p0_2d = ON_2dPoint::UnsetPoint;
2750  ON_2dPoint p1_2d = ON_2dPoint::UnsetPoint;
2751  ON_3dPoint p0_3d = data.curve->PointAt(t0);
2752  ON_3dPoint p1_3d = data.curve->PointAt(t1);
2753  ON_Interval dom[2];
2754  double p0_distance;
2755  double p1_distance;
2756 
2757  dom[0] = surf->Domain(0);
2758  dom[1] = surf->Domain(1);
2759 
2760  ON_3dPoint check_pt_3d = ON_3dPoint::UnsetPoint;
2761  int udir=0;
2762  int vdir=0;
2763  if (surface_GetClosestPoint3dFirstOrder(surf,p0_3d,p0_2d,check_pt_3d,p0_distance,0,same_point_tol,within_distance_tol) &&
2764  surface_GetClosestPoint3dFirstOrder(surf,p1_3d,p1_2d,check_pt_3d,p1_distance,0,same_point_tol,within_distance_tol) ) {
2765  if (ConsecutivePointsCrossClosedSeam(surf,p0_2d,p1_2d,udir,vdir,tol)) {
2766  ON_2dPoint p_2d;
2767  //lets check to see if p0 || p1 are already on a seam
2768  int seam0=0;
2769  if ((seam0 = IsAtSeam(surf,p0_2d, tol)) > 0) {
2770  ForceToClosestSeam(surf, p0_2d, tol);
2771  }
2772  int seam1 = 0;
2773  if ((seam1 = IsAtSeam(surf,p1_2d, tol)) > 0) {
2774  ForceToClosestSeam(surf, p1_2d, tol);
2775  }
2776  if (seam0 > 0 ) {
2777  if (seam1 > 0) { // both p0 & p1 on seam shouldn't happen report error and return false
2778  rc = false;
2779  } else { // just p0 on seam
2780  from = to = p0_2d;
2781  seam_t = t0;
2782  SwapUVSeamPoint(surf, to);
2783  }
2784  } else if (seam1 > 0) { // only p1 on seam
2785  from = to = p1_2d;
2786  seam_t = t1;
2787  SwapUVSeamPoint(surf, from);
2788  } else { // crosses the seam somewhere in between the two points
2789  bool seam_not_found = true;
2790  while(seam_not_found) {
2791  double d0 = DistToNearestClosedSeam(surf,p0_2d);
2792  double d1 = DistToNearestClosedSeam(surf,p1_2d);
2793  if ((d0 > 0.0) && (d1 > 0.0)) {
2794  double t = t0 + (t1 - t0)*(d0/(d0+d1));
2795  int seam;
2796  ON_3dPoint p_3d = data.curve->PointAt(t);
2797  double distance;
2798 
2799  if (surface_GetClosestPoint3dFirstOrder(surf,p_3d,p_2d,check_pt_3d,distance,0,same_point_tol,within_distance_tol)) {
2800  if ((seam=IsAtSeam(surf,p_2d, tol)) > 0) {
2801  ForceToClosestSeam(surf, p_2d, tol);
2802  from = to = p_2d;
2803  seam_t = t;
2804  if (p0_2d.DistanceTo(p_2d) < p1_2d.DistanceTo(p_2d)) {
2805  SwapUVSeamPoint(surf, to);
2806  } else {
2807  SwapUVSeamPoint(surf, from);
2808  }
2809  seam_not_found=false;
2810  rc = true;
2811  } else {
2812  if (ConsecutivePointsCrossClosedSeam(surf,p0_2d,p_2d,udir,vdir,tol)) {
2813  p1_2d = p_2d;
2814  t1 = t;
2815  } else if (ConsecutivePointsCrossClosedSeam(surf,p_2d,p1_2d,udir,vdir,tol)) {
2816  p0_2d = p_2d;
2817  t0=t;
2818  } else {
2819  seam_not_found=false;
2820  rc = false;
2821  }
2822  }
2823  } else {
2824  seam_not_found=false;
2825  rc = false;
2826  }
2827  } else {
2828  seam_not_found=false;
2829  rc = false;
2830  }
2831  }
2832  }
2833  }
2834  }
2835  }
2836 
2837  return rc;
2838 }
2839 
2840 
2841 /*
2842  * Find where 2D trim curve crosses closed seam of surface UV
2843  */
2844 bool
2845 FindTrimSeamCrossing(const ON_BrepTrim &trim,double t0,double t1,double &seam_t,ON_2dPoint &from,ON_2dPoint &to,double tol)
2846 {
2847  bool rc = true;
2848  const ON_Surface *surf = trim.SurfaceOf();
2849 
2850  // quick bail out is surface not closed
2851  if (surf->IsClosed(0) || surf->IsClosed(1)) {
2852  ON_2dPoint p0 = trim.PointAt(t0);
2853  ON_2dPoint p1 = trim.PointAt(t1);
2854  ON_Interval dom[2];
2855  dom[0] = surf->Domain(0);
2856  dom[1] = surf->Domain(1);
2857 
2858  p0 = UnwrapUVPoint(surf,p0);
2859  p1 = UnwrapUVPoint(surf,p1);
2860 
2861  int udir=0;
2862  int vdir=0;
2863  if (ConsecutivePointsCrossClosedSeam(surf,p0,p1,udir,vdir,tol)) {
2864  ON_2dPoint p;
2865  //lets check to see if p0 || p1 are already on a seam
2866  int seam0=0;
2867  if ((seam0 = IsAtSeam(surf,p0, tol)) > 0) {
2868  ForceToClosestSeam(surf, p0, tol);
2869  }
2870  int seam1 = 0;
2871  if ((seam1 = IsAtSeam(surf,p1, tol)) > 0) {
2872  ForceToClosestSeam(surf, p1, tol);
2873  }
2874  if (seam0 > 0 ) {
2875  if (seam1 > 0) { // both p0 & p1 on seam shouldn't happen report error and return false
2876  rc = false;
2877  } else { // just p0 on seam
2878  from = to = p0;
2879  seam_t = t0;
2880  SwapUVSeamPoint(surf, to);
2881  }
2882  } else if (seam1 > 0) { // only p1 on seam
2883  from = to = p1;
2884  seam_t = t1;
2885  SwapUVSeamPoint(surf, from);
2886  } else { // crosses the seam somewhere in between the two points
2887  bool seam_not_found = true;
2888  while (seam_not_found) {
2889  double d0 = DistToNearestClosedSeam(surf,p0);
2890  double d1 = DistToNearestClosedSeam(surf,p1);
2891  if ((d0 > tol) && (d1 > tol)) {
2892  double t = t0 + (t1 - t0)*(d0/(d0+d1));
2893  int seam;
2894  p = trim.PointAt(t);
2895  if ((seam=IsAtSeam(surf,p, tol)) > 0) {
2896  ForceToClosestSeam(surf, p, tol);
2897  from = to = p;
2898  seam_t = t;
2899 
2900  if (p0.DistanceTo(p) < p1.DistanceTo(p)) {
2901  SwapUVSeamPoint(surf, to);
2902  } else {
2903  SwapUVSeamPoint(surf, from);
2904  }
2905  seam_not_found=false;
2906  rc = true;
2907  } else {
2908  if (ConsecutivePointsCrossClosedSeam(surf,p0,p,udir,vdir,tol)) {
2909  p1 = p;
2910  t1 = t;
2911  } else if (ConsecutivePointsCrossClosedSeam(surf,p,p1,udir,vdir,tol)) {
2912  p0 = p;
2913  t0=t;
2914  } else {
2915  seam_not_found=false;
2916  rc = false;
2917  }
2918  }
2919  } else {
2920  seam_not_found=false;
2921  rc = false;
2922  }
2923  }
2924  }
2925  }
2926  }
2927 
2928  return rc;
2929 }
2930 
2931 
2932 void
2934  double t,
2935  double s,
2936  double same_point_tol,
2937  double within_distance_tol)
2938 {
2939  if (!data)
2940  return;
2941 
2942  if (!data->surf || !data->curve)
2943  return;
2944 
2945  const ON_Curve* curve= data->curve;
2946  const ON_Surface *surf = data->surf;
2947 
2948  ON_2dPointArray *samples= new ON_2dPointArray();
2949  size_t numKnots = curve->SpanCount();
2950  double *knots = new double[numKnots+1];
2951 
2952  curve->GetSpanVector(knots);
2953 
2954  size_t istart = 0;
2955  while ((istart < (numKnots+1)) && (t >= knots[istart]))
2956  istart++;
2957 
2958  if (istart > 0) {
2959  knots[--istart] = t;
2960  }
2961 
2962  size_t istop = numKnots;
2963  while ((istop > 0) && (s <= knots[istop]))
2964  istop--;
2965 
2966  if (istop < numKnots) {
2967  knots[++istop] = s;
2968  }
2969 
2970  size_t degree = curve->Degree();
2971  size_t samplesperknotinterval=18*degree;
2972 
2973  ON_2dPoint pt;
2974  ON_2dPoint prev_pt;
2975  double prev_t = knots[istart];
2976  double offset = 0.0;
2977  double delta;
2978  for (size_t i=istart; i<istop; i++) {
2979  delta = (knots[i+1] - knots[i])/(double)samplesperknotinterval;
2980  if (i <= numKnots/2) {
2981  offset = PBC_FROM_OFFSET;
2982  } else {
2983  offset = -PBC_FROM_OFFSET;
2984  }
2985  for (size_t j=0; j<=samplesperknotinterval; j++) {
2986  if ((j == samplesperknotinterval) && (i < istop - 1))
2987  continue;
2988 
2989  double curr_t = knots[i]+j*delta;
2990  if (curr_t < (s-t)/2.0) {
2991  offset = PBC_FROM_OFFSET;
2992  } else {
2993  offset = -PBC_FROM_OFFSET;
2994  }
2995  ON_3dPoint p = curve->PointAt(curr_t);
2996  ON_3dPoint p3d = ON_3dPoint::UnsetPoint;
2997  double distance;
2998  if (surface_GetClosestPoint3dFirstOrder(surf,p,pt,p3d,distance,0,same_point_tol,within_distance_tol)) {
2999  if (IsAtSeam(surf,pt,PBC_SEAM_TOL) > 0) {
3000  ForceToClosestSeam(surf, pt, PBC_SEAM_TOL);
3001  }
3002  if ((i == istart) && (j == 0)) {
3003  // first point just append and set reference in prev_pt
3004  samples->Append(pt);
3005  prev_pt = pt;
3006  prev_t = curr_t;
3007  continue;
3008  }
3009  int udir= 0;
3010  int vdir= 0;
3011  if (ConsecutivePointsCrossClosedSeam(surf,pt,prev_pt,udir,vdir,PBC_SEAM_TOL)) {
3012  int pt_seam = surf->IsAtSeam(pt.x,pt.y);
3013  int prev_pt_seam = surf->IsAtSeam(prev_pt.x,prev_pt.y);
3014  if ( pt_seam > 0) {
3015  if ((prev_pt_seam > 0) && (samples->Count() == 1)) {
3016  samples->Empty();
3017  SwapUVSeamPoint(surf, prev_pt,pt_seam);
3018  samples->Append(prev_pt);
3019  } else {
3020  if (pt_seam == 3) {
3021  if (prev_pt_seam == 1) {
3022  pt.x = prev_pt.x;
3023  if (ConsecutivePointsCrossClosedSeam(surf,pt,prev_pt,udir,vdir,PBC_SEAM_TOL)) {
3024  SwapUVSeamPoint(surf, pt, 2);
3025  }
3026  } else if (prev_pt_seam == 2) {
3027  pt.y = prev_pt.y;
3028  if (ConsecutivePointsCrossClosedSeam(surf,pt,prev_pt,udir,vdir,PBC_SEAM_TOL)) {
3029  SwapUVSeamPoint(surf, pt, 1);
3030  }
3031  }
3032  } else {
3033  SwapUVSeamPoint(surf, pt, prev_pt_seam);
3034  }
3035  }
3036  } else if (prev_pt_seam > 0) {
3037  if (samples->Count() == 1) {
3038  samples->Empty();
3039  SwapUVSeamPoint(surf, prev_pt);
3040  samples->Append(prev_pt);
3041  }
3042  } else if (data->curve->IsClosed()) {
3043  ON_2dPoint from,to;
3044  double seam_t;
3045  if (Find3DCurveSeamCrossing(*data,prev_t,curr_t,offset,seam_t,from,to,PBC_TOL,same_point_tol,within_distance_tol)) {
3046  samples->Append(from);
3047  data->segments.push_back(samples);
3048  samples= new ON_2dPointArray();
3049  samples->Append(to);
3050  prev_pt = to;
3051  prev_t = seam_t;
3052  } else {
3053  std::cout << "Can not find seam crossing...." << std::endl;
3054  }
3055  }
3056  }
3057  samples->Append(pt);
3058 
3059  prev_pt = pt;
3060  prev_t = curr_t;
3061  }
3062  }
3063  }
3064  delete [] knots;
3065 
3066  if (samples != NULL) {
3067  data->segments.push_back(samples);
3068 
3069  int numsegs = data->segments.size();
3070 
3071  if (numsegs > 1) {
3072  if (curve->IsClosed()) {
3073  ON_2dPointArray *reordered_samples= new ON_2dPointArray();
3074  // must have walked over seam but have closed curve so reorder stitching
3075  int seg = 0;
3076  for (std::list<ON_2dPointArray *>::reverse_iterator rit=data->segments.rbegin(); rit!=data->segments.rend(); ++seg) {
3077  samples = *rit;
3078  if (seg < numsegs-1) { // since end points should be repeated
3079  reordered_samples->Append(samples->Count()-1,(const ON_2dPoint *)samples->Array());
3080  } else {
3081  reordered_samples->Append(samples->Count(),(const ON_2dPoint *)samples->Array());
3082  }
3083  data->segments.erase((++rit).base());
3084  rit = data->segments.rbegin();
3085  delete samples;
3086  }
3087  data->segments.clear();
3088  data->segments.push_back(reordered_samples);
3089 
3090  } else {
3091  //punt for now
3092  }
3093  }
3094  }
3095 
3096  return;
3097 }
3098 
3099 
3100 PBCData *
3101 pullback_samples(const ON_Surface* surf,
3102  const ON_Curve* curve,
3103  double tolerance,
3104  double flatness,
3105  double same_point_tol,
3106  double within_distance_tol)
3107 {
3108  if (!surf)
3109  return NULL;
3110 
3111  PBCData *data = new PBCData;
3112  data->tolerance = tolerance;
3113  data->flatness = flatness;
3114  data->curve = curve;
3115  data->surf = surf;
3116  data->surftree = NULL;
3117 
3118  double tmin, tmax;
3119  data->curve->GetDomain(&tmin, &tmax);
3120 
3121  if (surf->IsClosed(0) || surf->IsClosed(1)) {
3122  if ((tmin < 0.0) && (tmax > 0.0)) {
3123  ON_2dPoint uv = ON_2dPoint::UnsetPoint;
3124  ON_3dPoint p = curve->PointAt(0.0);
3125  ON_3dPoint p3d = ON_3dPoint::UnsetPoint;
3126  double distance;
3127  int quadrant = 0; // optional - 0 = default, 1 from NE quadrant, 2 from NW quadrant, 3 from SW quadrant, 4 from SE quadrant
3128  if (surface_GetClosestPoint3dFirstOrder(surf,p,uv,p3d,distance,quadrant,same_point_tol,within_distance_tol)) {
3129  if (IsAtSeam(surf, uv, PBC_SEAM_TOL) > 0) {
3130  ON_2dPointArray *samples1 = pullback_samples(data, tmin, 0.0,same_point_tol,within_distance_tol);
3131  ON_2dPointArray *samples2 = pullback_samples(data, 0.0, tmax,same_point_tol,within_distance_tol);
3132  if (samples1 != NULL) {
3133  data->segments.push_back(samples1);
3134  }
3135  if (samples2 != NULL) {
3136  data->segments.push_back(samples2);
3137  }
3138  } else {
3139  ON_2dPointArray *samples = pullback_samples(data, tmin, tmax,same_point_tol,within_distance_tol);
3140  if (samples != NULL) {
3141  data->segments.push_back(samples);
3142  }
3143  }
3144  } else {
3145  std::cerr << "pullback_samples:Error: cannot evaluate curve at parameter 0.0" << std::endl;
3146  delete data;
3147  return NULL;
3148  }
3149  } else {
3150  pullback_samples_from_closed_surface(data, tmin, tmax,same_point_tol,within_distance_tol);
3151  }
3152  } else {
3153  ON_2dPointArray *samples = pullback_samples(data, tmin, tmax,same_point_tol,within_distance_tol);
3154  if (samples != NULL) {
3155  data->segments.push_back(samples);
3156  }
3157  }
3158  return data;
3159 }
3160 
3161 
3162 ON_Curve*
3163 refit_edge(const ON_BrepEdge* edge, double UNUSED(tolerance))
3164 {
3165  double edge_tolerance = 0.01;
3166  ON_Brep *brep = edge->Brep();
3167 #ifdef SHOW_UNUSED
3168  ON_3dPoint start = edge->PointAtStart();
3169  ON_3dPoint end = edge->PointAtEnd();
3170 #endif
3171 
3172  ON_BrepTrim& trim1 = brep->m_T[edge->m_ti[0]];
3173  ON_BrepTrim& trim2 = brep->m_T[edge->m_ti[1]];
3174  ON_BrepFace *face1 = trim1.Face();
3175  ON_BrepFace *face2 = trim2.Face();
3176  const ON_Surface *surface1 = face1->SurfaceOf();
3177  const ON_Surface *surface2 = face2->SurfaceOf();
3178  bool removeTrimmed = false;
3179  brlcad::SurfaceTree *st1 = new brlcad::SurfaceTree(face1, removeTrimmed);
3180  brlcad::SurfaceTree *st2 = new brlcad::SurfaceTree(face2, removeTrimmed);
3181 
3182  ON_Curve *curve = brep->m_C3[edge->m_c3i];
3183  double t0, t1;
3184  curve->GetDomain(&t0, &t1);
3185  ON_Plane plane;
3186  curve->FrameAt(t0, plane);
3187 #ifdef SHOW_UNUSED
3188  ON_3dPoint origin = plane.Origin();
3189  ON_3dVector xaxis = plane.Xaxis();
3190  ON_3dVector yaxis = plane.Yaxis();
3191  ON_3dVector zaxis = plane.zaxis;
3192  ON_3dPoint px = origin + xaxis;
3193  ON_3dPoint py = origin + yaxis;
3194  ON_3dPoint pz = origin + zaxis;
3195 #endif
3196 
3197  int numKnots = curve->SpanCount();
3198  double *knots = new double[numKnots + 1];
3199  curve->GetSpanVector(knots);
3200 
3201  int samplesperknotinterval;
3202  int degree = curve->Degree();
3203 
3204  if (degree > 1) {
3205  samplesperknotinterval = 3 * degree;
3206  } else {
3207  samplesperknotinterval = 18 * degree;
3208  }
3209  ON_2dPoint pt;
3210  double t = 0.0;
3211  ON_3dPoint pointOnCurve;
3212  ON_3dPoint knudgedPointOnCurve;
3213  for (int i = 0; i <= numKnots; i++) {
3214  if (i <= numKnots / 2) {
3215  if (i > 0) {
3216  double delta = (knots[i] - knots[i - 1]) / (double) samplesperknotinterval;
3217  for (int j = 1; j < samplesperknotinterval; j++) {
3218  t = knots[i - 1] + j * delta;
3219  pointOnCurve = curve->PointAt(t);
3220  knudgedPointOnCurve = curve->PointAt(t + PBC_TOL);
3221 
3222  ON_3dPoint point = pointOnCurve;
3223  ON_3dPoint knudgepoint = knudgedPointOnCurve;
3224  ON_3dPoint ps1;
3225  ON_3dPoint ps2;
3226  bool found = false;
3227  double dist;
3228  while (!found) {
3229  ON_2dPoint uv;
3230  if (st1->getSurfacePoint((const ON_3dPoint&) point, uv, (const ON_3dPoint&) knudgepoint, edge_tolerance) > 0) {
3231  ps1 = surface1->PointAt(uv.x, uv.y);
3232  if (st2->getSurfacePoint((const ON_3dPoint&) point, uv, (const ON_3dPoint&) knudgepoint, edge_tolerance) > 0) {
3233  ps2 = surface2->PointAt(uv.x, uv.y);
3234  }
3235  }
3236  dist = ps1.DistanceTo(ps2);
3237  if (NEAR_ZERO(dist, PBC_TOL)) {
3238  point = ps1;
3239  found = true;
3240  } else {
3241  ON_3dVector v1 = ps1 - point;
3242  ON_3dVector v2 = ps2 - point;
3243  knudgepoint = point;
3244  ON_3dVector deltav = v1 + v2;
3245  if (NEAR_ZERO(deltav.Length(), PBC_TOL)) {
3246  found = true; // as close as we are going to get
3247  } else {
3248  point = point + v1 + v2;
3249  }
3250  }
3251  }
3252  }
3253  }
3254  t = knots[i];
3255  pointOnCurve = curve->PointAt(t);
3256  knudgedPointOnCurve = curve->PointAt(t + PBC_TOL);
3257  ON_3dPoint point = pointOnCurve;
3258  ON_3dPoint knudgepoint = knudgedPointOnCurve;
3259  ON_3dPoint ps1;
3260  ON_3dPoint ps2;
3261  bool found = false;
3262  double dist;
3263 
3264  while (!found) {
3265  ON_2dPoint uv;
3266  if (st1->getSurfacePoint((const ON_3dPoint&) point, uv, (const ON_3dPoint&) knudgepoint, edge_tolerance) > 0) {
3267  ps1 = surface1->PointAt(uv.x, uv.y);
3268  if (st2->getSurfacePoint((const ON_3dPoint&) point, uv, (const ON_3dPoint&) knudgepoint, edge_tolerance) > 0) {
3269  ps2 = surface2->PointAt(uv.x, uv.y);
3270  }
3271  }
3272  dist = ps1.DistanceTo(ps2);
3273  if (NEAR_ZERO(dist, PBC_TOL)) {
3274  point = ps1;
3275  found = true;
3276  } else {
3277  ON_3dVector v1 = ps1 - point;
3278  ON_3dVector v2 = ps2 - point;
3279  knudgepoint = point;
3280  ON_3dVector deltav = v1 + v2;
3281  if (NEAR_ZERO(deltav.Length(), PBC_TOL)) {
3282  found = true; // as close as we are going to get
3283  } else {
3284  point = point + v1 + v2;
3285  }
3286  }
3287  }
3288  } else {
3289  if (i > 0) {
3290  double delta = (knots[i] - knots[i - 1]) / (double) samplesperknotinterval;
3291  for (int j = 1; j < samplesperknotinterval; j++) {
3292  t = knots[i - 1] + j * delta;
3293  pointOnCurve = curve->PointAt(t);
3294  knudgedPointOnCurve = curve->PointAt(t - PBC_TOL);
3295 
3296  ON_3dPoint point = pointOnCurve;
3297  ON_3dPoint knudgepoint = knudgedPointOnCurve;
3298  ON_3dPoint ps1;
3299  ON_3dPoint ps2;
3300  bool found = false;
3301  double dist;
3302 
3303  while (!found) {
3304  ON_2dPoint uv;
3305  if (st1->getSurfacePoint((const ON_3dPoint&) point, uv, (const ON_3dPoint&) knudgepoint, edge_tolerance) > 0) {
3306  ps1 = surface1->PointAt(uv.x, uv.y);
3307  if (st2->getSurfacePoint((const ON_3dPoint&) point, uv, (const ON_3dPoint&) knudgepoint, edge_tolerance) > 0) {
3308  ps2 = surface2->PointAt(uv.x, uv.y);
3309  }
3310  }
3311  dist = ps1.DistanceTo(ps2);
3312  if (NEAR_ZERO(dist, PBC_TOL)) {
3313  point = ps1;
3314  found = true;
3315  } else {
3316  ON_3dVector v1 = ps1 - point;
3317  ON_3dVector v2 = ps2 - point;
3318  knudgepoint = point;
3319  ON_3dVector deltav = v1 + v2;
3320  if (NEAR_ZERO(deltav.Length(), PBC_TOL)) {
3321  found = true; // as close as we are going to get
3322  } else {
3323  point = point + v1 + v2;
3324  }
3325  }
3326  }
3327  }
3328  t = knots[i];
3329  pointOnCurve = curve->PointAt(t);
3330  knudgedPointOnCurve = curve->PointAt(t - PBC_TOL);
3331  ON_3dPoint point = pointOnCurve;
3332  ON_3dPoint knudgepoint = knudgedPointOnCurve;
3333  ON_3dPoint ps1;
3334  ON_3dPoint ps2;
3335  bool found = false;
3336  double dist;
3337 
3338  while (!found) {
3339  ON_2dPoint uv;
3340  if (st1->getSurfacePoint((const ON_3dPoint&) point, uv, (const ON_3dPoint&) knudgepoint, edge_tolerance) > 0) {
3341  ps1 = surface1->PointAt(uv.x, uv.y);
3342  if (st2->getSurfacePoint((const ON_3dPoint&) point, uv, (const ON_3dPoint&) knudgepoint, edge_tolerance) > 0) {
3343  ps2 = surface2->PointAt(uv.x, uv.y);
3344  }
3345  }
3346  dist = ps1.DistanceTo(ps2);
3347  if (NEAR_ZERO(dist, PBC_TOL)) {
3348  point = ps1;
3349  found = true;
3350  } else {
3351  ON_3dVector v1 = ps1 - point;
3352  ON_3dVector v2 = ps2 - point;
3353  knudgepoint = point;
3354  ON_3dVector deltav = v1 + v2;
3355  if (NEAR_ZERO(deltav.Length(), PBC_TOL)) {
3356  found = true; // as close as we are going to get
3357  } else {
3358  point = point + v1 + v2;
3359  }
3360  }
3361  }
3362  }
3363  }
3364  }
3365  delete [] knots;
3366 
3367 
3368  return NULL;
3369 }
3370 
3371 
3372 bool
3373 has_singularity(const ON_Surface *surf)
3374 {
3375  bool ret = false;
3376  if (UNLIKELY(!surf)) return ret;
3377  // 0 = south, 1 = east, 2 = north, 3 = west
3378  for (int i = 0; i < 4; i++) {
3379  if (surf->IsSingular(i)) {
3380  /*
3381  switch (i) {
3382  case 0:
3383  std::cout << "Singular South" << std::endl;
3384  break;
3385  case 1:
3386  std::cout << "Singular East" << std::endl;
3387  break;
3388  case 2:
3389  std::cout << "Singular North" << std::endl;
3390  break;
3391  case 3:
3392  std::cout << "Singular West" << std::endl;
3393  }
3394  */
3395  ret = true;
3396  }
3397  }
3398  return ret;
3399 }
3400 
3401 
3402 bool is_closed(const ON_Surface *surf)
3403 {
3404  bool ret = false;
3405  if (UNLIKELY(!surf)) return ret;
3406  // dir 0 = "s", 1 = "t"
3407  for (int i = 0; i < 2; i++) {
3408  if (surf->IsClosed(i)) {
3409 // switch (i) {
3410 // case 0:
3411 // std::cout << "Closed in U" << std::endl;
3412 // break;
3413 // case 1:
3414 // std::cout << "Closed in V" << std::endl;
3415 // }
3416  ret = true;
3417  }
3418  }
3419  return ret;
3420 }
3421 
3422 
3423 bool
3424 check_pullback_closed(std::list<PBCData*> &pbcs)
3425 {
3426  std::list<PBCData*>::iterator d = pbcs.begin();
3427  if ((*d) == NULL || (*d)->surf == NULL)
3428  return false;
3429 
3430  const ON_Surface *surf = (*d)->surf;
3431 
3432  //TODO:
3433  // 0 = U, 1 = V
3434  if (surf->IsClosed(0) && surf->IsClosed(1)) {
3435  //TODO: need to check how torus UV looks to determine checks
3436  std::cerr << "Is this some kind of torus????" << std::endl;
3437  } else if (surf->IsClosed(0)) {
3438  //check_pullback_closed_U(pbcs);
3439  std::cout << "check closed in U" << std::endl;
3440  } else if (surf->IsClosed(1)) {
3441  //check_pullback_closed_V(pbcs);
3442  std::cout << "check closed in V" << std::endl;
3443  }
3444  return true;
3445 }
3446 
3447 
3448 bool
3449 check_pullback_singular_east(std::list<PBCData*> &pbcs)
3450 {
3451  std::list<PBCData *>::iterator cs = pbcs.begin();
3452  if ((*cs) == NULL || (*cs)->surf == NULL)
3453  return false;
3454 
3455  const ON_Surface *surf = (*cs)->surf;
3456  double umin, umax;
3457  ON_2dPoint *prev = NULL;
3458 
3459  surf->GetDomain(0, &umin, &umax);
3460  std::cout << "Umax: " << umax << std::endl;
3461  while (cs != pbcs.end()) {
3462  PBCData *data = (*cs);
3463  std::list<ON_2dPointArray *>::iterator si = data->segments.begin();
3464  int segcnt = 0;
3465  while (si != data->segments.end()) {
3466  ON_2dPointArray *samples = (*si);
3467  std::cerr << std::endl << "Segment:" << ++segcnt << std::endl;
3468  if (true) {
3469  int ilast = samples->Count() - 1;
3470  std::cerr << std::endl << 0 << "- " << (*samples)[0].x << ", " << (*samples)[0].y << std::endl;
3471  std::cerr << ilast << "- " << (*samples)[ilast].x << ", " << (*samples)[ilast].y << std::endl;
3472  } else {
3473  for (int i = 0; i < samples->Count(); i++) {
3474  if (NEAR_EQUAL((*samples)[i].x, umax, PBC_TOL)) {
3475  if (prev != NULL) {
3476  std::cerr << "prev - " << prev->x << ", " << prev->y << std::endl;
3477  }
3478  std::cerr << i << "- " << (*samples)[i].x << ", " << (*samples)[i].y << std::endl << std::endl;
3479  }
3480  prev = &(*samples)[i];
3481  }
3482  }
3483  si ++;
3484  }
3485  cs++;
3486  }
3487  return true;
3488 }
3489 
3490 
3491 bool
3492 check_pullback_singular(std::list<PBCData*> &pbcs)
3493 {
3494  std::list<PBCData*>::iterator d = pbcs.begin();
3495  if ((*d) == NULL || (*d)->surf == NULL)
3496  return false;
3497 
3498  const ON_Surface *surf = (*d)->surf;
3499  int cnt = 0;
3500 
3501  for (int i = 0; i < 4; i++) {
3502  if (surf->IsSingular(i)) {
3503  cnt++;
3504  }
3505  }
3506  if (cnt > 2) {
3507  //TODO: I don't think this makes sense but check out
3508  std::cerr << "Is this some kind of sickness????" << std::endl;
3509  return false;
3510  } else if (cnt == 2) {
3511  if (surf->IsSingular(0) && surf->IsSingular(2)) {
3512  std::cout << "check singular North-South" << std::endl;
3513  } else if (surf->IsSingular(1) && surf->IsSingular(2)) {
3514  std::cout << "check singular East-West" << std::endl;
3515  } else {
3516  //TODO: I don't think this makes sense but check out
3517  std::cerr << "Is this some kind of sickness????" << std::endl;
3518  return false;
3519  }
3520  } else {
3521  if (surf->IsSingular(0)) {
3522  std::cout << "check singular South" << std::endl;
3523  } else if (surf->IsSingular(1)) {
3524  std::cout << "check singular East" << std::endl;
3525  if (check_pullback_singular_east(pbcs)) {
3526  return true;
3527  }
3528  } else if (surf->IsSingular(2)) {
3529  std::cout << "check singular North" << std::endl;
3530  } else if (surf->IsSingular(3)) {
3531  std::cout << "check singular West" << std::endl;
3532  }
3533  }
3534  return true;
3535 }
3536 
3537 
3538 #ifdef _DEBUG_TESTING_
3539 void
3540 print_pullback_data(std::string str, std::list<PBCData*> &pbcs, bool justendpoints)
3541 {
3542  std::list<PBCData*>::iterator cs = pbcs.begin();
3543  int trimcnt = 0;
3544  if (justendpoints) {
3545  // print out endpoints before
3546  std::cerr << "EndPoints " << str << ":" << std::endl;
3547  while (cs != pbcs.end()) {
3548  PBCData *data = (*cs);
3549  if (!data || !data->surf)
3550  continue;
3551 
3552  const ON_Surface *surf = data->surf;
3553  std::list<ON_2dPointArray *>::iterator si = data->segments.begin();
3554  int segcnt = 0;
3555  while (si != data->segments.end()) {
3556  ON_2dPointArray *samples = (*si);
3557  std::cerr << std::endl << " Segment:" << ++segcnt << std::endl;
3558  int ilast = samples->Count() - 1;
3559  std::cerr << " T:" << ++trimcnt << std::endl;
3560  int i = 0;
3561  int singularity = IsAtSingularity(surf, (*samples)[i], PBC_SEAM_TOL);
3562  int seam = IsAtSeam(surf, (*samples)[i], PBC_SEAM_TOL);
3563  std::cerr << "--------";
3564  if ((seam > 0) && (singularity >= 0)) {
3565  std::cerr << " S/S " << (*samples)[i].x << ", " << (*samples)[i].y;
3566  } else if (seam > 0) {
3567  std::cerr << " Seam " << (*samples)[i].x << ", " << (*samples)[i].y;
3568  } else if (singularity >= 0) {
3569  std::cerr << " Sing " << (*samples)[i].x << ", " << (*samples)[i].y;
3570  } else {
3571  std::cerr << " " << (*samples)[i].x << ", " << (*samples)[i].y;
3572  }
3573  ON_3dPoint p = surf->PointAt((*samples)[i].x, (*samples)[i].y);
3574  std::cerr << " (" << p.x << ", " << p.y << ", " << p.z << ") " << std::endl;
3575 
3576  i = ilast;
3577  singularity = IsAtSingularity(surf, (*samples)[i], PBC_SEAM_TOL);
3578  seam = IsAtSeam(surf, (*samples)[i], PBC_SEAM_TOL);
3579  std::cerr << " ";
3580  if ((seam > 0) && (singularity >= 0)) {
3581  std::cerr << " S/S " << (*samples)[i].x << ", " << (*samples)[i].y << std::endl;
3582  } else if (seam > 0) {
3583  std::cerr << " Seam " << (*samples)[i].x << ", " << (*samples)[i].y << std::endl;
3584  } else if (singularity >= 0) {
3585  std::cerr << " Sing " << (*samples)[i].x << ", " << (*samples)[i].y << std::endl;
3586  } else {
3587  std::cerr << " " << (*samples)[i].x << ", " << (*samples)[i].y << std::endl;
3588  }
3589  p = surf->PointAt((*samples)[i].x, (*samples)[i].y);
3590  std::cerr << " (" << p.x << ", " << p.y << ", " << p.z << ") " << std::endl;
3591  si++;
3592  }
3593  cs++;
3594  }
3595  } else {
3596  // print out all points
3597  trimcnt = 0;
3598  cs = pbcs.begin();
3599  std::cerr << str << ":" << std::endl;
3600  while (cs != pbcs.end()) {
3601  PBCData *data = (*cs);
3602  if (!data || !data->surf)
3603  continue;
3604 
3605  const ON_Surface *surf = data->surf;
3606  std::list<ON_2dPointArray *>::iterator si = data->segments.begin();
3607  int segcnt = 0;
3608  std::cerr << "2d surface domain: " << std::endl;
3609  std::cerr << "in rpp rpp" << surf->Domain(0).m_t[0] << " " << surf->Domain(0).m_t[1] << " " << surf->Domain(1).m_t[0] << " " << surf->Domain(1).m_t[1] << " 0.0 0.01" << std::endl;
3610  while (si != data->segments.end()) {
3611  ON_2dPointArray *samples = (*si);
3612  std::cerr << std::endl << " Segment:" << ++segcnt << std::endl;
3613  std::cerr << " T:" << ++trimcnt << std::endl;
3614  for (int i = 0; i < samples->Count(); i++) {
3615  int singularity = IsAtSingularity(surf, (*samples)[i], PBC_SEAM_TOL);
3616  int seam = IsAtSeam(surf, (*samples)[i], PBC_SEAM_TOL);
3618  std::cerr << "in pt_" << _debug_point_count_++ << " sph " << (*samples)[i].x << " " << (*samples)[i].y << " 0.0 0.1000" << std::endl;
3619  } else if (_debug_print_mged3d_points_) {
3620  ON_3dPoint p = surf->PointAt((*samples)[i].x, (*samples)[i].y);
3621  std::cerr << "in pt_" << _debug_point_count_++ << " sph " << p.x << " " << p.y << " " << p.z << " 0.1000" << std::endl;
3622  } else {
3623  if (i == 0) {
3624  std::cerr << "--------";
3625  } else {
3626  std::cerr << " ";
3627  }
3628  if ((seam > 0) && (singularity >= 0)) {
3629  std::cerr << " S/S " << (*samples)[i].x << ", " << (*samples)[i].y;
3630  } else if (seam > 0) {
3631  std::cerr << " Seam " << (*samples)[i].x << ", " << (*samples)[i].y;
3632  } else if (singularity >= 0) {
3633  std::cerr << " Sing " << (*samples)[i].x << ", " << (*samples)[i].y;
3634  } else {
3635  std::cerr << " " << (*samples)[i].x << ", " << (*samples)[i].y;
3636  }
3637  ON_3dPoint p = surf->PointAt((*samples)[i].x, (*samples)[i].y);
3638  std::cerr << " (" << p.x << ", " << p.y << ", " << p.z << ") " << std::endl;
3639  }
3640  }
3641  si++;
3642  }
3643  cs++;
3644  }
3645  }
3646  /////
3647 }
3648 #endif
3649 
3650 
3651 bool
3652 resolve_seam_segment_from_prev(const ON_Surface *surface, ON_2dPointArray &segment, ON_2dPoint *prev = NULL)
3653 {
3654  bool complete = false;
3655  double umin, umax, umid;
3656  double vmin, vmax, vmid;
3657 
3658  surface->GetDomain(0, &umin, &umax);
3659  surface->GetDomain(1, &vmin, &vmax);
3660  umid = (umin + umax) / 2.0;
3661  vmid = (vmin + vmax) / 2.0;
3662 
3663  for (int i = 0; i < segment.Count(); i++) {
3664  int singularity = IsAtSingularity(surface, segment[i], PBC_SEAM_TOL);
3665  int seam = IsAtSeam(surface, segment[i], PBC_SEAM_TOL);
3666  if ((seam > 0)) {
3667  if (prev != NULL) {
3668  //std::cerr << " at seam " << seam << " but has prev" << std::endl;
3669  //std::cerr << " prev: " << prev->x << ", " << prev->y << std::endl;
3670  //std::cerr << " curr: " << data->samples[i].x << ", " << data->samples[i].y << std::endl;
3671  switch (seam) {
3672  case 1: //east/west
3673  if (prev->x < umid) {
3674  segment[i].x = umin;
3675  } else {
3676  segment[i].x = umax;
3677  }
3678  break;
3679  case 2: //north/south
3680  if (prev->y < vmid) {
3681  segment[i].y = vmin;
3682  } else {
3683  segment[i].y = vmax;
3684  }
3685  break;
3686  case 3: //both
3687  if (prev->x < umid) {
3688  segment[i].x = umin;
3689  } else {
3690  segment[i].x = umax;
3691  }
3692  if (prev->y < vmid) {
3693  segment[i].y = vmin;
3694  } else {
3695  segment[i].y = vmax;
3696  }
3697  }
3698  prev = &segment[i];
3699  } else {
3700  //std::cerr << " at seam and no prev" << std::endl;
3701  complete = false;
3702  }
3703  } else {
3704  if (singularity < 0) {
3705  prev = &segment[i];
3706  } else {
3707  prev = NULL;
3708  }
3709  }
3710  }
3711  return complete;
3712 }
3713 
3714 
3715 bool
3716 resolve_seam_segment_from_next(const ON_Surface *surface, ON_2dPointArray &segment, ON_2dPoint *next = NULL)
3717 {
3718  bool complete = false;
3719  double umin, umax, umid;
3720  double vmin, vmax, vmid;
3721 
3722  surface->GetDomain(0, &umin, &umax);
3723  surface->GetDomain(1, &vmin, &vmax);
3724  umid = (umin + umax) / 2.0;
3725  vmid = (vmin + vmax) / 2.0;
3726 
3727  if (next != NULL) {
3728  complete = true;
3729  for (int i = segment.Count() - 1; i >= 0; i--) {
3730  int singularity = IsAtSingularity(surface, segment[i], PBC_SEAM_TOL);
3731  int seam = IsAtSeam(surface, segment[i], PBC_SEAM_TOL);
3732 
3733  if ((seam > 0)) {
3734  if (next != NULL) {
3735  switch (seam) {
3736  case 1: //east/west
3737  if (next->x < umid) {
3738  segment[i].x = umin;
3739  } else {
3740  segment[i].x = umax;
3741  }
3742  break;
3743  case 2: //north/south
3744  if (next->y < vmid) {
3745  segment[i].y = vmin;
3746  } else {
3747  segment[i].y = vmax;
3748  }
3749  break;
3750  case 3: //both
3751  if (next->x < umid) {
3752  segment[i].x = umin;
3753  } else {
3754  segment[i].x = umax;
3755  }
3756  if (next->y < vmid) {
3757  segment[i].y = vmin;
3758  } else {
3759  segment[i].y = vmax;
3760  }
3761  }
3762  next = &segment[i];
3763  } else {
3764  //std::cerr << " at seam and no prev" << std::endl;
3765  complete = false;
3766  }
3767  } else {
3768  if (singularity < 0) {
3769  next = &segment[i];
3770  } else {
3771  next = NULL;
3772  }
3773  }
3774  }
3775  }
3776  return complete;
3777 }
3778 
3779 
3780 bool
3781 resolve_seam_segment(const ON_Surface *surface, ON_2dPointArray &segment, bool &UNUSED(u_resolved), bool &UNUSED(v_resolved))
3782 {
3783  ON_2dPoint *prev = NULL;
3784  bool complete = false;
3785  double umin, umax, umid;
3786  double vmin, vmax, vmid;
3787  int prev_seam = 0;
3788 
3789  surface->GetDomain(0, &umin, &umax);
3790  surface->GetDomain(1, &vmin, &vmax);
3791  umid = (umin + umax) / 2.0;
3792  vmid = (vmin + vmax) / 2.0;
3793 
3794  for (int i = 0; i < segment.Count(); i++) {
3795  int singularity = IsAtSingularity(surface, segment[i], PBC_SEAM_TOL);
3796  int seam = IsAtSeam(surface, segment[i], PBC_SEAM_TOL);
3797 
3798  if ((seam > 0)) {
3799  if (prev != NULL) {
3800  //std::cerr << " at seam " << seam << " but has prev" << std::endl;
3801  //std::cerr << " prev: " << prev->x << ", " << prev->y << std::endl;
3802  //std::cerr << " curr: " << data->samples[i].x << ", " << data->samples[i].y << std::endl;
3803  switch (seam) {
3804  case 1: //east/west
3805  if (prev->x < umid) {
3806  segment[i].x = umin;
3807  } else {
3808  segment[i].x = umax;
3809  }
3810  break;
3811  case 2: //north/south
3812  if (prev->y < vmid) {
3813  segment[i].y = vmin;
3814  } else {
3815  segment[i].y = vmax;
3816  }
3817  break;
3818  case 3: //both
3819  if (prev->x < umid) {
3820  segment[i].x = umin;
3821  } else {
3822  segment[i].x = umax;
3823  }
3824  if (prev->y < vmid) {
3825  segment[i].y = vmin;
3826  } else {
3827  segment[i].y = vmax;
3828  }
3829  }
3830  prev = &segment[i];
3831  } else {
3832  //std::cerr << " at seam and no prev" << std::endl;
3833  complete = false;
3834  }
3835  prev_seam = seam;
3836  } else {
3837  if (singularity < 0) {
3838  prev = &segment[i];
3839  prev_seam = 0;
3840  } else {
3841  prev = NULL;
3842  }
3843  }
3844  }
3845  prev_seam = 0;
3846  if ((!complete) && (prev != NULL)) {
3847  complete = true;
3848  for (int i = segment.Count() - 2; i >= 0; i--) {
3849  int singularity = IsAtSingularity(surface, segment[i], PBC_SEAM_TOL);
3850  int seam = IsAtSeam(surface, segment[i], PBC_SEAM_TOL);
3851  if ((seam > 0)) {
3852  if (prev != NULL) {
3853  //std::cerr << " at seam " << seam << " but has prev" << std::endl;
3854  //std::cerr << " prev: " << prev->x << ", " << prev->y << std::endl;
3855  //std::cerr << " curr: " << data->samples[i].x << ", " << data->samples[i].y << std::endl;
3856  switch (seam) {
3857  case 1: //east/west
3858  if (prev->x < umid) {
3859  segment[i].x = umin;
3860  } else {
3861  segment[i].x = umax;
3862  }
3863  break;
3864  case 2: //north/south
3865  if (prev->y < vmid) {
3866  segment[i].y = vmin;
3867  } else {
3868  segment[i].y = vmax;
3869  }
3870  break;
3871  case 3: //both
3872  if (prev->x < umid) {
3873  segment[i].x = umin;
3874  } else {
3875  segment[i].x = umax;
3876  }
3877  if (prev->y < vmid) {
3878  segment[i].y = vmin;
3879  } else {
3880  segment[i].y = vmax;
3881  }
3882  }
3883  prev = &segment[i];
3884  } else {
3885  //std::cerr << " at seam and no prev" << std::endl;
3886  complete = false;
3887  }
3888  } else {
3889  if (singularity < 0) {
3890  prev = &segment[i];
3891  } else {
3892  prev = NULL;
3893  }
3894  }
3895  }
3896  }
3897  return complete;
3898 }
3899 
3900 
3901 /*
3902  * number_of_seam_crossings
3903  */
3904 int
3905 number_of_seam_crossings(std::list<PBCData*> &pbcs)
3906 {
3907  int rc = 0;
3908  std::list<PBCData*>::iterator cs;
3909 
3910  cs = pbcs.begin();
3911  while (cs != pbcs.end()) {
3912  PBCData *data = (*cs);
3913  if (!data || !data->surf)
3914  continue;
3915 
3916  const ON_Surface *surf = data->surf;
3917  std::list<ON_2dPointArray *>::iterator si = data->segments.begin();
3918  ON_2dPoint *pt = NULL;
3919  ON_2dPoint *prev_pt = NULL;
3920  ON_2dPoint *next_pt = NULL;
3921  while (si != data->segments.end()) {
3922  ON_2dPointArray *samples = (*si);
3923  for (int i = 0; i < samples->Count(); i++) {
3924  pt = &(*samples)[i];
3925  if (!IsAtSeam(surf,*pt,PBC_SEAM_TOL)) {
3926  if (prev_pt == NULL) {
3927  prev_pt = pt;
3928  } else {
3929  next_pt = pt;
3930  }
3931  int udir=0;
3932  int vdir=0;
3933  if (next_pt != NULL) {
3934  if (ConsecutivePointsCrossClosedSeam(surf,*prev_pt,*next_pt,udir,vdir,PBC_SEAM_TOL)) {
3935  rc++;
3936  }
3937  prev_pt = next_pt;
3938  next_pt = NULL;
3939  }
3940  }
3941  }
3942  if (si != data->segments.end())
3943  si++;
3944  }
3945 
3946  if (cs != pbcs.end())
3947  cs++;
3948  }
3949 
3950  return rc;
3951 }
3952 
3953 
3954 /*
3955  * if current and previous point on seam make sure they are on same seam
3956  */
3957 bool
3958 check_for_points_on_same_seam(std::list<PBCData*> &pbcs)
3959 {
3960 
3961  std::list<PBCData*>::iterator cs = pbcs.begin();
3962  ON_2dPoint *prev_pt = NULL;
3963  int prev_seam = 0;
3964  while ( cs != pbcs.end()) {
3965  PBCData *data = (*cs);
3966  const ON_Surface *surf = data->surf;
3967  std::list<ON_2dPointArray *>::iterator seg = data->segments.begin();
3968  while (seg != data->segments.end()) {
3969  ON_2dPointArray *points = (*seg);
3970  for (int i=0; i < points->Count(); i++) {
3971  ON_2dPoint *pt = points->At(i);
3972  int seam = IsAtSeam(surf,*pt,PBC_SEAM_TOL);
3973  if (seam > 0) {
3974  if (prev_seam > 0) {
3975  if ((seam == 1) && ((prev_seam % 2) == 1)) {
3976  pt->x = prev_pt->x;
3977  } else if ((seam == 2) && (prev_seam > 1)) {
3978  pt->y = prev_pt->y;
3979  } else if (seam == 3) {
3980  if ((prev_seam % 2) == 1) {
3981  pt->x = prev_pt->x;
3982  }
3983  if (prev_seam > 1) {
3984  pt->y = prev_pt->y;
3985  }
3986  }
3987  }
3988  prev_seam = seam;
3989  prev_pt = pt;
3990  }
3991  }
3992  seg++;
3993  }
3994  cs++;
3995  }
3996  return true;
3997 }
3998 
3999 
4000 /*
4001  * extend_pullback_at_shared_3D_curve_seam
4002  */
4003 bool
4005 {
4006  const ON_Curve *next_curve = NULL;
4007  std::set<const ON_Curve *> set;
4008  std::map<const ON_Curve *,int> map;
4009  std::list<PBCData*>::iterator cs = pbcs.begin();
4010 
4011  while ( cs != pbcs.end()) {
4012  PBCData *data = (*cs++);
4013  const ON_Curve *curve = data->curve;
4014  const ON_Surface *surf = data->surf;
4015 
4016  if (cs != pbcs.end()) {
4017  PBCData *nextdata = (*cs);
4018  next_curve = nextdata->curve;
4019  }
4020 
4021  if (curve == next_curve) {
4022  std::cerr << "Consecutive seam usage" << std::endl;
4023  //find which direction we need to extend
4024  if (surf->IsClosed(0) && !surf->IsClosed(1)) {
4025  double length = surf->Domain(0).Length();
4026  std::list<ON_2dPointArray *>::iterator seg = data->segments.begin();
4027  while (seg != data->segments.end()) {
4028  ON_2dPointArray *points = (*seg);
4029  for (int i=0; i < points->Count(); i++) {
4030  points->At(i)->x = points->At(i)->x + length;
4031  }
4032  seg++;
4033  }
4034  } else if (!surf->IsClosed(0) && surf->IsClosed(1)) {
4035  double length = surf->Domain(1).Length();
4036  std::list<ON_2dPointArray *>::iterator seg = data->segments.begin();
4037  while (seg != data->segments.end()) {
4038  ON_2dPointArray *points = (*seg);
4039  for (int i=0; i < points->Count(); i++) {
4040  points->At(i)->y = points->At(i)->y + length;
4041  }
4042  seg++;
4043  }
4044  } else {
4045  std::cerr << "both directions" << std::endl;
4046  }
4047  }
4048  next_curve = NULL;
4049  }
4050  return true;
4051 }
4052 
4053 
4054 /*
4055  * shift_closed_curve_split_over_seam
4056  */
4057 bool
4059 {
4060  if (pbcs.size() == 1) { // single curve for this loop
4061  std::list<PBCData*>::iterator cs;
4062 
4063  PBCData *data = pbcs.front();
4064  if (!data || !data->surf)
4065  return false;
4066 
4067  const ON_Surface *surf = data->surf;
4068  ON_Interval udom = surf->Domain(0);
4069  ON_Interval vdom = surf->Domain(1);
4070  std::list<ON_2dPointArray *>::iterator si = data->segments.begin();
4071  ON_2dPoint pt;
4072  ON_2dPoint prev_pt;
4073  if (data->curve->IsClosed()) {
4074  int numseamcrossings = number_of_seam_crossings(pbcs);
4075  if (numseamcrossings == 1) {
4076  ON_2dPointArray part1,part2;
4077  ON_2dPointArray* curr_point_array = &part2;
4078  while (si != data->segments.end()) {
4079  ON_2dPointArray *samples = (*si);
4080  for (int i = 0; i < samples->Count(); i++) {
4081  pt = (*samples)[i];
4082  if (i == 0) {
4083  prev_pt = pt;
4084  curr_point_array->Append(pt);
4085  continue;
4086  }
4087  int udir= 0;
4088  int vdir= 0;
4089  if (ConsecutivePointsCrossClosedSeam(surf,pt,prev_pt,udir,vdir,PBC_SEAM_TOL)) {
4090  if (surf->IsAtSeam(pt.x,pt.y) > 0) {
4091  SwapUVSeamPoint(surf, pt);
4092  curr_point_array->Append(pt);
4093  curr_point_array = &part1;
4094  SwapUVSeamPoint(surf, pt);
4095  } else if (surf->IsAtSeam(prev_pt.x,prev_pt.y) > 0) {
4096  SwapUVSeamPoint(surf, prev_pt);
4097  curr_point_array->Append(prev_pt);
4098  } else {
4099  std::cerr << "shift_single_curve_loop_straddled_over_seam(): Error expecting to see seam in sample points" << std::endl;
4100  }
4101  }
4102  curr_point_array->Append(pt);
4103  prev_pt = pt;
4104  }
4105  samples->Empty();
4106  samples->Append(part1.Count(),part1.Array());
4107  samples->Append(part2.Count(),part2.Array());
4108  if (si != data->segments.end())
4109  si++;
4110  }
4111  }
4112  }
4113  }
4114  return true;
4115 }
4116 
4117 
4118 /*
4119  * extend_over_seam_crossings - all incoming points assumed to be within UV bounds
4120  */
4121 bool
4122 extend_over_seam_crossings(std::list<PBCData*> &pbcs)
4123 {
4124  std::list<PBCData*>::iterator cs;
4125  ON_2dPoint *pt = NULL;
4126  ON_2dPoint *prev_pt = NULL;
4127 
4128  ///// Loop through and fix any seam ambiguities
4129  cs = pbcs.begin();
4130  while (cs != pbcs.end()) {
4131  PBCData *data = (*cs);
4132  if (!data || !data->surf)
4133  continue;
4134 
4135  std::list<ON_2dPointArray *>::iterator si = data->segments.begin();
4136  while (si != data->segments.end()) {
4137  ON_2dPointArray *samples = (*si);
4138  for (int i = 0; i < samples->Count(); i++) {
4139  pt = &(*samples)[i];
4140 
4141  if (prev_pt != NULL) {
4142  GetClosestExtendedPoint(data->surf,*pt,*prev_pt,PBC_SEAM_TOL);
4143  }
4144  prev_pt = pt;
4145  }
4146  if (si != data->segments.end())
4147  si++;
4148  }
4149  if (cs != pbcs.end())
4150  cs++;
4151  }
4152 
4153  return true;
4154 }
4155 
4156 
4157 /*
4158  * run through curve loop to determine correct start/end
4159  * points resolving ambiguities when point lies on a seam or
4160  * singularity
4161  */
4162 bool
4163 resolve_pullback_seams(std::list<PBCData*> &pbcs)
4164 {
4165  std::list<PBCData*>::iterator cs;
4166 
4167  ///// Loop through and fix any seam ambiguities
4168  ON_2dPoint *prev = NULL;
4169  ON_2dPoint *next = NULL;
4170  bool u_resolved = false;
4171  bool v_resolved = false;
4172  cs = pbcs.begin();
4173  while (cs != pbcs.end()) {
4174  PBCData *data = (*cs);
4175  if (!data || !data->surf)
4176  continue;
4177 
4178  const ON_Surface *surf = data->surf;
4179  double umin, umax;
4180  double vmin, vmax;
4181  surf->GetDomain(0, &umin, &umax);
4182  surf->GetDomain(1, &vmin, &vmax);
4183 
4184  std::list<ON_2dPointArray *>::iterator si = data->segments.begin();
4185  while (si != data->segments.end()) {
4186  ON_2dPointArray *samples = (*si);
4187  if (resolve_seam_segment(surf, *samples,u_resolved,v_resolved)) {
4188  // Found a starting point
4189  //1) walk back up with resolved next point
4190  next = (*samples).First();
4191  std::list<PBCData*>::reverse_iterator rcs(cs);
4192  rcs--;
4193  std::list<ON_2dPointArray *>::reverse_iterator rsi(si);
4194  while (rcs != pbcs.rend()) {
4195  PBCData *rdata = (*rcs);
4196  while (rsi != rdata->segments.rend()) {
4197  ON_2dPointArray *rsamples = (*rsi);
4198  // first try and resolve on own merits
4199  if (!resolve_seam_segment(surf, *rsamples,u_resolved,v_resolved)) {
4200  resolve_seam_segment_from_next(surf, *rsamples, next);
4201  }
4202  next = (*rsamples).First();
4203  rsi++;
4204  }
4205  rcs++;
4206  if (rcs != pbcs.rend()) {
4207  rdata = (*rcs);
4208  rsi = rdata->segments.rbegin();
4209  }
4210  }
4211 
4212  //2) walk rest of way down with resolved prev point
4213  if (samples->Count() > 1)
4214  prev = &(*samples)[samples->Count() - 1];
4215  else
4216  prev = NULL;
4217  si++;
4218  std::list<PBCData*>::iterator current(cs);
4219  while (cs != pbcs.end()) {
4220  while (si != data->segments.end()) {
4221  samples = (*si);
4222  // first try and resolve on own merits
4223  if (!resolve_seam_segment(surf, *samples,u_resolved,v_resolved)) {
4224  resolve_seam_segment_from_prev(surf, *samples, prev);
4225  }
4226  if (samples->Count() > 1)
4227  prev = &(*samples)[samples->Count() - 1];
4228  else
4229  prev = NULL;
4230  si++;
4231  }
4232  cs++;
4233  if (cs != pbcs.end()) {
4234  data = (*cs);
4235  si = data->segments.begin();
4236  }
4237  }
4238  // make sure to wrap back around with previous
4239  cs = pbcs.begin();
4240  data = (*cs);
4241  si = data->segments.begin();
4242  while ((cs != pbcs.end()) && (cs != current)) {
4243  while (si != data->segments.end()) {
4244  samples = (*si);
4245  // first try and resolve on own merits
4246  if (!resolve_seam_segment(surf, *samples,u_resolved,v_resolved)) {
4247  resolve_seam_segment_from_prev(surf, *samples, prev);
4248  }
4249  if (samples->Count() > 1)
4250  prev = &(*samples)[samples->Count() - 1];
4251  else
4252  prev = NULL;
4253  si++;
4254  }
4255  cs++;
4256  if (cs != pbcs.end()) {
4257  data = (*cs);
4258  si = data->segments.begin();
4259  }
4260  }
4261  }
4262  if (si != data->segments.end())
4263  si++;
4264  }
4265  if (cs != pbcs.end())
4266  cs++;
4267  }
4268  return true;
4269 }
4270 
4271 
4272 /*
4273  * run through curve loop to determine correct start/end
4274  * points resolving ambiguities when point lies on a seam or
4275  * singularity
4276  */
4277 bool
4278 resolve_pullback_singularities(std::list<PBCData*> &pbcs)
4279 {
4280  std::list<PBCData*>::iterator cs = pbcs.begin();
4281 
4282  ///// Loop through and fix any seam ambiguities
4283  ON_2dPoint *prev = NULL;
4284  bool complete = false;
4285  int checkcnt = 0;
4286 
4287  prev = NULL;
4288  complete = false;
4289  checkcnt = 0;
4290  while (!complete && (checkcnt < 2)) {
4291  cs = pbcs.begin();
4292  complete = true;
4293  checkcnt++;
4294  //std::cerr << "Checkcnt - " << checkcnt << std::endl;
4295  while (cs != pbcs.end()) {
4296  int singularity;
4297  prev = NULL;
4298  PBCData *data = (*cs);
4299  if (!data || !data->surf)
4300  continue;
4301 
4302  const ON_Surface *surf = data->surf;
4303  std::list<ON_2dPointArray *>::iterator si = data->segments.begin();
4304  while (si != data->segments.end()) {
4305  ON_2dPointArray *samples = (*si);
4306  for (int i = 0; i < samples->Count(); i++) {
4307  // 0 = south, 1 = east, 2 = north, 3 = west
4308  if ((singularity = IsAtSingularity(surf, (*samples)[i], PBC_SEAM_TOL)) >= 0) {
4309  if (prev != NULL) {
4310  //std::cerr << " at singularity " << singularity << " but has prev" << std::endl;
4311  //std::cerr << " prev: " << prev->x << ", " << prev->y << std::endl;
4312  //std::cerr << " curr: " << data->samples[i].x << ", " << data->samples[i].y << std::endl;
4313  switch (singularity) {
4314  case 0: //south
4315  (*samples)[i].x = prev->x;
4316  (*samples)[i].y = surf->Domain(1)[0];
4317  break;
4318  case 1: //east
4319  (*samples)[i].y = prev->y;
4320  (*samples)[i].x = surf->Domain(0)[1];
4321  break;
4322  case 2: //north
4323  (*samples)[i].x = prev->x;
4324  (*samples)[i].y = surf->Domain(1)[1];
4325  break;
4326  case 3: //west
4327  (*samples)[i].y = prev->y;
4328  (*samples)[i].x = surf->Domain(0)[0];
4329  }
4330  prev = NULL;
4331  //std::cerr << " curr now: " << data->samples[i].x << ", " << data->samples[i].y << std::endl;
4332  } else {
4333  //std::cerr << " at singularity " << singularity << " and no prev" << std::endl;
4334  //std::cerr << " curr: " << data->samples[i].x << ", " << data->samples[i].y << std::endl;
4335  complete = false;
4336  }
4337  } else {
4338  prev = &(*samples)[i];
4339  }
4340  }
4341  if (!complete) {
4342  //std::cerr << "Lets work backward:" << std::endl;
4343  for (int i = samples->Count() - 2; i >= 0; i--) {
4344  // 0 = south, 1 = east, 2 = north, 3 = west
4345  if ((singularity = IsAtSingularity(surf, (*samples)[i], PBC_SEAM_TOL)) >= 0) {
4346  if (prev != NULL) {
4347  //std::cerr << " at singularity " << singularity << " but has prev" << std::endl;
4348  //std::cerr << " prev: " << prev->x << ", " << prev->y << std::endl;
4349  //std::cerr << " curr: " << data->samples[i].x << ", " << data->samples[i].y << std::endl;
4350  switch (singularity) {
4351  case 0: //south
4352  (*samples)[i].x = prev->x;
4353  (*samples)[i].y = surf->Domain(1)[0];
4354  break;
4355  case 1: //east
4356  (*samples)[i].y = prev->y;
4357  (*samples)[i].x = surf->Domain(0)[1];
4358  break;
4359  case 2: //north
4360  (*samples)[i].x = prev->x;
4361  (*samples)[i].y = surf->Domain(1)[1];
4362  break;
4363  case 3: //west
4364  (*samples)[i].y = prev->y;
4365  (*samples)[i].x = surf->Domain(0)[0];
4366  }
4367  prev = NULL;
4368  //std::cerr << " curr now: " << data->samples[i].x << ", " << data->samples[i].y << std::endl;
4369  } else {
4370  //std::cerr << " at singularity " << singularity << " and no prev" << std::endl;
4371  //std::cerr << " curr: " << data->samples[i].x << ", " << data->samples[i].y << std::endl;
4372  complete = false;
4373  }
4374  } else {
4375  prev = &(*samples)[i];
4376  }
4377  }
4378  }
4379  si++;
4380  }
4381  cs++;
4382  }
4383  }
4384 
4385  return true;
4386 }
4387 
4388 
4389 void
4391 {
4392  std::list<PBCData*>::iterator cs = pbcs.begin();
4393  while (cs != pbcs.end()) {
4394  PBCData *data = (*cs);
4395  std::list<ON_2dPointArray *>::iterator si = data->segments.begin();
4396  while (si != data->segments.end()) {
4397  ON_2dPointArray *samples = (*si);
4398  if (samples->Count() == 0) {
4399  si = data->segments.erase(si);
4400  } else {
4401  for (int i = 0; i < samples->Count() - 1; i++) {
4402  while ((i < (samples->Count() - 1)) && (*samples)[i].DistanceTo((*samples)[i + 1]) < 1e-9) {
4403  samples->Remove(i + 1);
4404  }
4405  }
4406  si++;
4407  }
4408  }
4409  if (data->segments.empty()) {
4410  cs = pbcs.erase(cs);
4411  } else {
4412  cs++;
4413  }
4414  }
4415 }
4416 
4417 
4418 bool
4419 check_pullback_data(std::list<PBCData*> &pbcs)
4420 {
4421  std::list<PBCData*>::iterator d = pbcs.begin();
4422 
4423  if (d == pbcs.end() || (*d) == NULL || (*d)->surf == NULL)
4424  return false;
4425 
4426  const ON_Surface *surf = (*d)->surf;
4427  bool singular = has_singularity(surf);
4428  bool closed = is_closed(surf);
4429 
4430  if (singular) {
4431  if (!resolve_pullback_singularities(pbcs)) {
4432  std::cerr << "Error: Can not resolve singular ambiguities." << std::endl;
4433  }
4434  }
4435 
4436  if (closed) {
4437  // check for same 3D curve use
4438  if (!check_for_points_on_same_seam(pbcs)) {
4439  std::cerr << "Error: Can not extend pullback at shared 3D curve seam." << std::endl;
4440  return false;
4441  }
4442  // check for same 3D curve use
4444  std::cerr << "Error: Can not extend pullback at shared 3D curve seam." << std::endl;
4445  return false;
4446  }
4448  std::cerr << "Error: Can not resolve seam ambiguities." << std::endl;
4449  return false;
4450  }
4451  if (!resolve_pullback_seams(pbcs)) {
4452  std::cerr << "Error: Can not resolve seam ambiguities." << std::endl;
4453  return false;
4454  }
4455  if (!extend_over_seam_crossings(pbcs)) {
4456  std::cerr << "Error: Can not resolve seam ambiguities." << std::endl;
4457  return false;
4458  }
4459  }
4460 
4461  // consecutive duplicates within segment will cause problems in curve fit
4463 
4464  return true;
4465 }
4466 
4467 
4468 int
4469 check_pullback_singularity_bridge(const ON_Surface *surf, const ON_2dPoint &p1, const ON_2dPoint &p2)
4470 {
4471  if (has_singularity(surf)) {
4472  int is, js;
4473  if (((is = IsAtSingularity(surf, p1, PBC_SEAM_TOL)) >= 0) && ((js = IsAtSingularity(surf, p2, PBC_SEAM_TOL)) >= 0)) {
4474  //create new singular trim
4475  if (is == js) {
4476  return is;
4477  }
4478  }
4479  }
4480  return -1;
4481 }
4482 
4483 
4484 int
4485 check_pullback_seam_bridge(const ON_Surface *surf, const ON_2dPoint &p1, const ON_2dPoint &p2)
4486 {
4487  if (is_closed(surf)) {
4488  int is, js;
4489  if (((is = IsAtSeam(surf, p1, PBC_SEAM_TOL)) > 0) && ((js = IsAtSeam(surf, p2, PBC_SEAM_TOL)) > 0)) {
4490  //create new seam trim
4491  if (is == js) {
4492  // need to check if seam 3d points are equal
4493  double endpoint_distance = p1.DistanceTo(p2);
4494  double t0, t1;
4495 
4496  int dir = is - 1;
4497  surf->GetDomain(dir, &t0, &t1);
4498  if (endpoint_distance > 0.5 * (t1 - t0)) {
4499  return is;
4500  }
4501  }
4502  }
4503  }
4504  return -1;
4505 }
4506 
4507 
4508 ON_Curve*
4509 pullback_curve(const brlcad::SurfaceTree* surfacetree,
4510  const ON_Surface* surf,
4511  const ON_Curve* curve,
4512  double tolerance,
4513  double flatness)
4514 {
4515  PBCData data;
4516  data.tolerance = tolerance;
4517  data.flatness = flatness;
4518  data.curve = curve;
4519  data.surf = surf;
4520  data.surftree = (brlcad::SurfaceTree*)surfacetree;
4521  ON_2dPointArray samples;
4522  data.segments.push_back(&samples);
4523 
4524  // Step 1 - adaptively sample the curve
4525  double tmin, tmax;
4526  data.curve->GetDomain(&tmin, &tmax);
4527  ON_2dPoint& start = samples.AppendNew(); // new point is added to samples and returned
4528  if (!toUV(data, start, tmin, 0.001)) {
4529  return NULL; // fails if first point is out of tolerance!
4530  }
4531 
4532  ON_2dPoint uv;
4533  ON_3dPoint p = curve->PointAt(tmin);
4534  ON_3dPoint from = curve->PointAt(tmin + 0.0001);
4535  brlcad::SurfaceTree *st = (brlcad::SurfaceTree *)surfacetree;
4536  if (st->getSurfacePoint((const ON_3dPoint&)p, uv, (const ON_3dPoint&)from) < 0) {
4537  std::cerr << "Error: Can not get surface point." << std::endl;
4538  }
4539 
4540  ON_2dPoint p1, p2;
4541 
4542 #ifdef SHOW_UNUSED
4543  if (!data.surf)
4544  return NULL;
4545 
4546  const ON_Surface *surf = data.surf;
4547 #endif
4548 
4549  if (toUV(data, p1, tmin, PBC_TOL) && toUV(data, p2, tmax, -PBC_TOL)) {
4550 #ifdef SHOW_UNUSED
4551  ON_3dPoint a = surf->PointAt(p1.x, p1.y);
4552  ON_3dPoint b = surf->PointAt(p2.x, p2.y);
4553 #endif
4554 
4555  p = curve->PointAt(tmax);
4556  from = curve->PointAt(tmax - 0.0001);
4557  if (st->getSurfacePoint((const ON_3dPoint&)p, uv, (const ON_3dPoint&)from) < 0) {
4558  std::cerr << "Error: Can not get surface point." << std::endl;
4559  }
4560 
4561  if (!sample(data, tmin, tmax, p1, p2)) {
4562  return NULL;
4563  }
4564 
4565  for (int i = 0; i < samples.Count(); i++) {
4566  std::cerr << samples[i].x << ", " << samples[i].y << std::endl;
4567  }
4568  } else {
4569  return NULL;
4570  }
4571 
4572  return interpolateCurve(samples);
4573 }
4574 
4575 
4576 ON_Curve*
4578  const brlcad::SurfaceTree* surfacetree,
4579  const ON_Surface* surf,
4580  const ON_Curve* curve,
4581  double tolerance,
4582  double flatness)
4583 {
4584  PBCData data;
4585  data.tolerance = tolerance;
4586  data.flatness = flatness;
4587  data.curve = curve;
4588  data.surf = surf;
4589  data.surftree = (brlcad::SurfaceTree*)surfacetree;
4590  ON_2dPointArray samples;
4591  data.segments.push_back(&samples);
4592 
4593  // Step 1 - adaptively sample the curve
4594  double tmin, tmax;
4595  data.curve->GetDomain(&tmin, &tmax);
4596  ON_2dPoint& start = samples.AppendNew(); // new point is added to samples and returned
4597  if (!toUV(data, start, tmin, 0.001)) {
4598  return NULL; // fails if first point is out of tolerance!
4599  }
4600 
4601  ON_2dPoint p1, p2;
4602 
4603  if (toUV(data, p1, tmin, PBC_TOL) && toUV(data, p2, tmax, -PBC_TOL)) {
4604  if (!sample(data, tmin, tmax, p1, p2)) {
4605  return NULL;
4606  }
4607 
4608  for (int i = 0; i < samples.Count(); i++) {
4609  if (seam_dir == NORTH_SEAM) {
4610  samples[i].y = 1.0;
4611  } else if (seam_dir == EAST_SEAM) {
4612  samples[i].x = 1.0;
4613  } else if (seam_dir == SOUTH_SEAM) {
4614  samples[i].y = 0.0;
4615  } else if (seam_dir == WEST_SEAM) {
4616  samples[i].x = 0.0;
4617  }
4618  std::cerr << samples[i].x << ", " << samples[i].y << std::endl;
4619  }
4620  } else {
4621  return NULL;
4622  }
4623 
4624  return interpolateCurve(samples);
4625 }
4626 
4627 
4628 // Local Variables:
4629 // tab-width: 8
4630 // mode: C++
4631 // c-basic-offset: 4
4632 // indent-tabs-mode: t
4633 // c-file-style: "stroustrup"
4634 // End:
4635 // ex: shiftwidth=4 tabstop=8
ustring wrap
double surface_GetOptimalNormalVSplit(const ON_Surface *surf, const ON_Interval &u_interval, const ON_Interval &v_interval, double tol)
ON_Curve * interpolateCurve(ON_2dPointArray &samples)
ON_2dPoint UnwrapUVPoint(const ON_Surface *surf, const ON_2dPoint &pt, double tol)
double surface_GetClosestPoint3dFirstOrderSubdivision(const ON_Surface *surf, const ON_3dPoint &p, const ON_Interval &u_interval, double u, const ON_Interval &v_interval, double v, double current_closest_dist, ON_2dPoint &p2d, ON_3dPoint &p3d, double same_point_tol, double within_distance_tol, int level)
ON_Curve * refit_edge(const ON_BrepEdge *edge, double tolerance)
std::vector< double > knots
if lu s
Definition: nmg_mod.c:3860
#define N
Definition: randmt.c:39
void bu_semaphore_acquire(unsigned int i)
Definition: semaphore.c:180
int check_pullback_seam_bridge(const ON_Surface *surf, const ON_2dPoint &p1, const ON_2dPoint &p2)
Definition: raytrace.h:368
#define st
bool check_for_points_on_same_seam(std::list< PBCData * > &pbcs)
#define PBC_SEAM_TOL
Definition: PullbackCurve.h:69
bool isFlat(const ON_2dPoint &p1, const ON_2dPoint &m, const ON_2dPoint &p2, double flatness)
bool Find3DCurveSeamCrossing(PBCData &data, double t0, double t1, double offset, double &seam_t, ON_2dPoint &from, ON_2dPoint &to, double tol, double same_point_tol, double within_distance_tol)
Header file for the BRL-CAD common definitions.
void brep_r(const ON_Surface *surf, plane_ray &pr, pt2d_t uv, ON_3dPoint &pt, ON_3dVector &su, ON_3dVector &sv, pt2d_t R)
ON_BOOL32 surface_GetIntervalMinMaxDistance(const ON_3dPoint &p, ON_BoundingBox &bbox, double &min_distance, double &max_distance)
ON_NurbsCurve * newNURBSCurve(BSpline &spline, int dimension=3)
double surface_GetClosestPoint3dFirstOrderByRange(const ON_Surface *surf, const ON_3dPoint &p, const ON_Interval &u_range, const ON_Interval &v_range, double current_closest_dist, ON_2dPoint &p2d, ON_3dPoint &p3d, double same_point_tol, double within_distance_tol, int level)
int IsAtSeam(const ON_Surface *surf, int dir, double u, double v, double tol)
void SwapUVSeamPoint(const ON_Surface *surf, ON_2dPoint &p, int hint)
#define MAX_FASTF
Definition: defines.h:340
bool check_pullback_data(std::list< PBCData * > &pbcs)
int getKnotInterval(BSpline &bspline, double u)
ON_2dPointArray * pullback_samples(PBCData *data, double t, double s, double same_point_tol, double within_distance_tol)
COMPLEX data[64]
Definition: fftest.c:34
bool check_pullback_singular(std::list< PBCData * > &pbcs)
ON_NurbsCurve * interpolateLocalCubicCurve(ON_2dPointArray &Q)
bool extend_pullback_at_shared_3D_curve_seam(std::list< PBCData * > &pbcs)
#define BREP_INTERSECTION_ROOT_EPSILON
Definition: brep.h:66
void brep_get_plane_ray(ON_Ray &r, plane_ray &pr)
double surface_GetOptimalNormalUSplit(const ON_Surface *surf, const ON_Interval &u_interval, const ON_Interval &v_interval, double tol)
bool resolve_seam_segment_from_next(const ON_Surface *surface, ON_2dPointArray &segment, ON_2dPoint *next=NULL)
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
bool _debug_print_mged3d_points_
ON_Curve * pullback_curve(const brlcad::SurfaceTree *surfacetree, const ON_Surface *surf, const ON_Curve *curve, double tolerance, double flatness)
bool FindTrimSeamCrossing(const ON_BrepTrim &trim, double t0, double t1, double &seam_t, ON_2dPoint &from, ON_2dPoint &to, double tol)
ON_BOOL32 surface_EvNormal(const ON_Surface *surf, double s, double t, ON_3dPoint &point, ON_3dVector &normal, int side, int *hint)
Coord * point
Definition: chull3d.cpp:52
#define PBC_FROM_OFFSET
Definition: PullbackCurve.h:68
void print_pullback_data(std::string str, std::list< PBCData * > &pbcs, bool justendpoints)
int bu_parallel_id(void)
Definition: parallel.c:152
ON_BOOL32 surface_GetBoundingBox(const ON_Surface *surf, const ON_Interval &u_interval, const ON_Interval &v_interval, ON_BoundingBox &bbox, ON_BOOL32 bGrowBox)
#define BU_SEM_MALLOC
Definition: parallel.h:183
#define UNUSED(parameter)
Definition: common.h:239
bool sample(PBCData &data, double t1, double t2, const ON_2dPoint &p1, const ON_2dPoint &p2)
#define RANGE_HI
goto out
Definition: nmg_mod.c:3846
void brep_newton_iterate(plane_ray &pr, pt2d_t R, ON_3dVector &su, ON_3dVector &sv, pt2d_t uv, pt2d_t out_uv)
seam_direction
Definition: PullbackCurve.h:59
bool _debug_print_mged2d_points_
double drand48(void)
ON_Curve * pullback_seam_curve(enum seam_direction seam_dir, const brlcad::SurfaceTree *surfacetree, const ON_Surface *surf, const ON_Curve *curve, double tolerance, double flatness)
void bu_semaphore_release(unsigned int i)
Definition: semaphore.c:218
bool extend_over_seam_crossings(std::list< PBCData * > &pbcs)
ustring alpha
bool surface_GetClosestPoint3dFirstOrder(const ON_Surface *surf, const ON_3dPoint &p, ON_2dPoint &p2d, ON_3dPoint &p3d, double &current_distance, int quadrant, double same_point_tol, double within_distance_tol)
#define PBC_TOL
Definition: PullbackCurve.h:67
bool resolve_seam_segment(const ON_Surface *surface, ON_2dPointArray &segment, bool &u_resolved, bool &v_resolved)
Definition: clone.c:124
int check_pullback_singularity_bridge(const ON_Surface *surf, const ON_2dPoint &p1, const ON_2dPoint &p2)
#define MAX_PSW
Definition: parallel.h:48
bool check_pullback_closed(std::list< PBCData * > &pbcs)
HIDDEN const double cs
Definition: sh_prj.c:617
void remove_consecutive_intersegment_duplicates(std::list< PBCData * > &pbcs)
bool trim_GetClosestPoint3dFirstOrder(const ON_BrepTrim &trim, const ON_3dPoint &p, ON_2dPoint &p2d, double &t, double &distance, const ON_Interval *interval, double same_point_tol, double within_distance_tol)
bool resolve_seam_segment_from_prev(const ON_Surface *surface, ON_2dPointArray &segment, ON_2dPoint *prev=NULL)
bool resolve_pullback_singularities(std::list< PBCData * > &pbcs)
double randomPointFromRange(PBCData &data, ON_2dPoint &out, double lo, double hi)
void utah_ray_planes(const ON_Ray &r, ON_3dVector &p1, double &p1d, ON_3dVector &p2, double &p2d)
bool shift_single_curve_loop_straddled_over_seam(std::list< PBCData * > &pbcs)
void GetClosestExtendedPoint(const ON_Surface *surf, ON_2dPoint &pt, ON_2dPoint &prev_pt, double tol)
ON_2dPointArray controls
#define RANGE_LO
bool ConsecutivePointsCrossClosedSeam(const ON_Surface *surf, const ON_2dPoint &pt, const ON_2dPoint &prev_pt, int &udir, int &vdir, double tol)
void ForceToClosestSeam(const ON_Surface *surf, ON_2dPoint &pt, double tol)
ON_BOOL32 face_GetBoundingBox(const ON_BrepFace &face, ON_BoundingBox &bbox, ON_BOOL32 bGrowBox)
bool toUV(PBCData &data, ON_2dPoint &out_pt, double t, double knudge=0.0, double within_distance_tol=BREP_EDGE_MISS_TOLERANCE)
#define A
Definition: msr.c:51
int _debug_point_count_
double DistToNearestClosedSeam(const ON_Surface *surf, const ON_2dPoint &pt)
#define Q
Definition: msr.c:54
bool resolve_pullback_seams(std::list< PBCData * > &pbcs)
bool has_singularity(const ON_Surface *surf)
int number_of_seam_crossings(std::list< PBCData * > &pbcs)
void pullback_samples_from_closed_surface(PBCData *data, double t, double s, double same_point_tol, double within_distance_tol)
bool is_closed(const ON_Surface *surf)
int IsAtSingularity(const ON_Surface *surf, double u, double v, double tol)
HIDDEN const point_t delta
Definition: sh_prj.c:618
std::vector< double > params
ON_BOOL32 GetDomainSplits(const ON_Surface *surf, const ON_Interval &u_interval, const ON_Interval &v_interval, ON_Interval domSplits[2][2])
double fastf_t
Definition: defines.h:300
enum seam_direction seam_direction(ON_2dPoint uv1, ON_2dPoint uv2)
void generateKnots(BSpline &bspline)
bool check_pullback_singular_east(std::list< PBCData * > &pbcs)
bool _debug_print_
#define UNLIKELY(expression)
Definition: common.h:282