xmsgeom  1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
geoms.cpp
Go to the documentation of this file.
1 //------------------------------------------------------------------------------
8 //------------------------------------------------------------------------------
9 
10 // 1. Precompiled header
11 
12 // 2. My own header
13 #include <xmsgeom/geometry/geoms.h>
14 
15 // 3. Standard library headers
16 #include <cfloat>
17 #include <algorithm>
18 
19 // 4. External library headers
20 
21 // 5. Shared code headers
22 #include <xmscore/math/math.h> // GT_TOL, LT_TOL
23 #include <xmscore/misc/XmError.h>
24 #include <xmscore/points/pt.h> // Pt3d
25 #include <xmscore/stl/vector.h>
26 #include <xmsgeom/geometry/GmPolygon.h>
27 #include <xmscore/misc/xmstype.h> // XM_ZERO_TOL
28 #include <xmsgeom/triangulate/triangles.h>
29 #include <xmscore/misc/XmConst.h>
30 
31 // 6. Non-shared code headers
32 
33 //----- Namespace declaration --------------------------------------------------
34 
35 namespace xms
36 {
37 //----- Constants / Enumerations -----------------------------------------------
38 
39 //----- Forward declarations ---------------------------------------------------
40 
41 //----- External globals -------------------------------------------------------
42 
43 //----- Internal function prototypes -------------------------------------------
44 
45 //------------------------------------------------------------------------------
51 //------------------------------------------------------------------------------
52 bool gmPointInOrOnBox2d(const Pt3d& a_bMin, const Pt3d& a_bMax, const Pt3d& a_pt)
53 {
54  if (a_pt.x < a_bMin.x || a_pt.y < a_bMin.y || a_pt.x > a_bMax.x || a_pt.y > a_bMax.y)
55  return false;
56  return true;
57 } // gmBoxesOverlap2d
58 //------------------------------------------------------------------------------
65 //------------------------------------------------------------------------------
66 bool gmBoxesOverlap2d(const Pt3d& a_b1Min,
67  const Pt3d& a_b1Max,
68  const Pt3d& a_b2Min,
69  const Pt3d& a_b2Max)
70 {
71  if (a_b1Max.x < a_b2Min.x)
72  return false;
73  if (a_b1Min.x > a_b2Max.x)
74  return false;
75  if (a_b1Max.y < a_b2Min.y)
76  return false;
77  if (a_b1Min.y > a_b2Max.y)
78  return false;
79  return true;
80 } // gmBoxesOverlap2d
81 //------------------------------------------------------------------------------
91 //------------------------------------------------------------------------------
92 void gmCalculateNormalizedPlaneCoefficients(const Pt3d& p1,
93  const Pt3d& p2,
94  const Pt3d& p3,
95  double* a,
96  double* b,
97  double* c,
98  double* d)
99 {
100  // call the other version
101  gmCalculateNormalizedPlaneCoefficients(&p1, &p2, &p3, a, b, c, d);
102 } // gmCalculateNormalizedPlaneCoefficients
103 //------------------------------------------------------------------------------
113 //------------------------------------------------------------------------------
114 void gmCalculateNormalizedPlaneCoefficients(const Pt3d* p1,
115  const Pt3d* p2,
116  const Pt3d* p3,
117  double* a,
118  double* b,
119  double* c,
120  double* d)
121 {
122  double x1(p1->x), y1(p1->y), z1(p1->z);
123  double x2(p2->x), y2(p2->y), z2(p2->z);
124  double x3(p3->x), y3(p3->y), z3(p3->z);
125  *a = (y1 * (z2 - z3) + y2 * (z3 - z1) + y3 * (z1 - z2));
126  *b = (z1 * (x2 - x3) + z2 * (x3 - x1) + z3 * (x1 - x2));
127  *c = (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2));
128  double mag(sqrt((*a) * (*a) + (*b) * (*b) + (*c) * (*c)));
129  *a /= mag;
130  *b /= mag;
131  *c /= mag;
132  *d = -(*a) * x1 - (*b) * y1 - (*c) * z1;
133 } // gmCalculateNormalizedPlaneCoefficients
134 //------------------------------------------------------------------------------
139 //------------------------------------------------------------------------------
140 double gmMiddleZ(const VecPt3d& a_points)
141 {
142  // Get an average z value
143  double zmin = XM_DBL_HIGHEST;
144  double zmax = XM_DBL_LOWEST;
145  for (size_t i = 0; i < a_points.size(); ++i)
146  {
147  double z = a_points[i].z;
148  zmin = std::min(zmin, z);
149  zmax = std::max(zmax, z);
150  }
151  return (zmin + zmax) / 2.0;
152 } // gmMiddleZ
153 //------------------------------------------------------------------------------
158 //------------------------------------------------------------------------------
159 PtInOutOrOn_enum gmPtInCircumcircle(const Pt3d& pt, Pt3d circumcirclePts[3])
160 {
161  double xc, yc, r2;
162 
163  if (!gmCircumcircleWithTol(&circumcirclePts[0], &circumcirclePts[1], &circumcirclePts[2], &xc,
164  &yc, &r2, gmXyTol()))
165  {
166  return PT_ERROR;
167  }
168  /* compute distance from (xc,yc) to pt squared */
169  double delta = sqrt(r2) - sqrt(sqr(pt.x - xc) + sqr(pt.y - yc));
170  if (fabs(delta) > gmXyTol())
171  {
172  if (delta > gmXyTol())
173  {
174  return PT_IN;
175  }
176  else
177  {
178  return PT_OUT;
179  }
180  }
181  return PT_ON;
182 } // gmPtInCircumcircle
183 //------------------------------------------------------------------------------
188 //------------------------------------------------------------------------------
189 double gmXyDistanceSquared(const Pt3d& pt1, const Pt3d& pt2)
190 {
191  return (sqr(pt1.x - pt2.x) + sqr(pt1.y - pt2.y));
192 } // gmXyDistanceSqared
193 //------------------------------------------------------------------------------
204 //------------------------------------------------------------------------------
205 bool gmCircumcircleWithTol(const Pt3d* pt1,
206  const Pt3d* pt2,
207  const Pt3d* pt3,
208  double* xc,
209  double* yc,
210  double* r2,
211  double tol)
212 {
213  bool ok = true;
214  double det11, det12, det13, det21, det22, det23;
215  double determinate;
216  /* compute these */
217  det11 = pt1->x - pt2->x;
218  det12 = pt1->y - pt2->y;
219  det13 = det11 * (pt1->x + pt2->x) / 2.0 + det12 * (pt1->y + pt2->y) / 2.0;
220 
221  det21 = pt3->x - pt2->x;
222  det22 = pt3->y - pt2->y;
223  det23 = det21 * (pt3->x + pt2->x) / 2.0 + det22 * (pt3->y + pt2->y) / 2.0;
224  /* compute determinant */
225  determinate = det11 * det22 - det12 * det21;
226  /* zero determinate indicates collinear pts */
227  if (fabs(determinate) < tol)
228  {
229  ok = false;
230  determinate = tol;
231  /* the old subtri perturbed a vertex to
232  * give a non zero determinant as follows:
233 i = 0;
234 while (determinate == 0.0) {
235 i++;
236 Nx = ((double)i/1000.0)*(vtx1->y - vtx0->y);
237 Ny = ((double)i/1000.0)*(vtx0->x - vtx1->x);
238 x1new = vtx1->x + Nx;
239 y1new = vtx1->y + Ny;
240 det11 = x1new - vtx0->x;
241 det12 = y1new - vtx0->y;
242 det13 = det11*(x1new + vtx0->x)/2.0 + det12*(y1new + vtx0->y)/2.0;
243 determinate = det11 * det22 - det12 * det21;
244 }
245  * END OF COMMENTED CODE */
246  }
247  *xc = (det13 * det22 - det23 * det12) / determinate;
248  *yc = (det11 * det23 - det21 * det13) / determinate;
249  *r2 = sqr(pt1->x - *xc) + sqr(pt1->y - *yc);
250  return ok;
251 } // gmCircumcircleWithTol
252 //------------------------------------------------------------------------------
261 //------------------------------------------------------------------------------
262 int gmCartToBary(const Pt3d* cart, const Pt3d* orig, double coef[6], int dir, Pt3d* bary)
263 {
264  double x, y, z;
265  x = cart->x - orig->x;
266  y = cart->y - orig->y;
267  z = cart->z - orig->z;
268  /* get magnitudes of the plane normal */
269  switch (dir)
270  {
271  case (0):
272  bary->x = coef[0] * y + coef[1] * z + coef[2];
273  bary->y = coef[3] * y + coef[4] * z + coef[5];
274  bary->z = 1.0 - bary->x - bary->y;
275  break;
276  case (1):
277  bary->x = coef[0] * z + coef[1] * x + coef[2];
278  bary->y = coef[3] * z + coef[4] * x + coef[5];
279  bary->z = 1.0 - bary->x - bary->y;
280  break;
281  case (2):
282  bary->x = coef[0] * x + coef[1] * y + coef[2];
283  bary->y = coef[3] * x + coef[4] * y + coef[5];
284  bary->z = 1.0 - bary->x - bary->y;
285  break;
286  }
287  return (XM_SUCCESS);
288 } // gmCartToBary
289 //------------------------------------------------------------------------------
301 //------------------------------------------------------------------------------
302 int gmBaryPrepare(const Pt3d* p1,
303  const Pt3d* p2,
304  const Pt3d* p3,
305  const Pt3d* norm,
306  Pt3d* orig,
307  double coef[6],
308  int* dir,
309  bool flag)
310 {
311  double x1, x2, x3, y1, y2, y3, z1, z2, z3, sizeinv;
312  /* compute the direction */
313  x1 = fabs(norm->x);
314  y1 = fabs(norm->y);
315  z1 = fabs(norm->z);
316  if (x1 > y1 && x1 > z1)
317  *dir = 0;
318  else if (y1 > z1)
319  *dir = 1;
320  else
321  *dir = 2;
322  /* get the origin */
323  if (flag)
324  {
325  orig->x = Mmin(p1->x, Mmin(p2->x, p3->x));
326  orig->y = Mmin(p1->y, Mmin(p2->y, p3->y));
327  orig->z = Mmin(p1->z, Mmin(p2->z, p3->z));
328  }
329  /* compute the coefficients */
330  x1 = p1->x - orig->x;
331  y1 = p1->y - orig->y;
332  z1 = p1->z - orig->z;
333  x2 = p2->x - orig->x;
334  y2 = p2->y - orig->y;
335  z2 = p2->z - orig->z;
336  x3 = p3->x - orig->x;
337  y3 = p3->y - orig->y;
338  z3 = p3->z - orig->z;
339  switch (*dir)
340  {
341  case (0):
342  sizeinv = 1.0 / ((y2 - y3) * (z1 - z3) - (z2 - z3) * (y1 - y3));
343  coef[0] = (z3 - z2) * sizeinv;
344  coef[1] = (y2 - y3) * sizeinv;
345  coef[2] = (y3 * z2 - z3 * y2) * sizeinv;
346  coef[3] = (z1 - z3) * sizeinv;
347  coef[4] = (y3 - y1) * sizeinv;
348  coef[5] = (y1 * z3 - z1 * y3) * sizeinv;
349  break;
350  case (1):
351  sizeinv = 1.0 / ((z2 - z3) * (x1 - x3) - (x2 - x3) * (z1 - z3));
352  coef[0] = (x3 - x2) * sizeinv;
353  coef[1] = (z2 - z3) * sizeinv;
354  coef[2] = (z3 * x2 - x3 * z2) * sizeinv;
355  coef[3] = (x1 - x3) * sizeinv;
356  coef[4] = (z3 - z1) * sizeinv;
357  coef[5] = (z1 * x3 - x1 * z3) * sizeinv;
358  break;
359  case (2):
360  sizeinv = 1.0 / ((x2 - x3) * (y1 - y3) - (y2 - y3) * (x1 - x3));
361  coef[0] = (y3 - y2) * sizeinv;
362  coef[1] = (x2 - x3) * sizeinv;
363  coef[2] = (x3 * y2 - y3 * x2) * sizeinv;
364  coef[3] = (y1 - y3) * sizeinv;
365  coef[4] = (x3 - x1) * sizeinv;
366  coef[5] = (x1 * y3 - y1 * x3) * sizeinv;
367  break;
368  }
369  return (XM_SUCCESS);
370 } // gmBaryPrepare
371 //------------------------------------------------------------------------------
379 //------------------------------------------------------------------------------
380 bool gmColinearWithTol(const Pt3d& p1, const Pt3d& p2, const Pt3d& p3, const double tol)
381 {
382  if (gmOnLineWithTol(p1, p2, p3.x, p3.y, tol))
383  return true;
384  else if (gmOnLineWithTol(p2, p3, p1.x, p1.y, tol))
385  return true;
386  else if (gmOnLineWithTol(p3, p1, p2.x, p2.y, tol))
387  return true;
388  else
389  return false;
390 } // gmColinearWithTol
391 //------------------------------------------------------------------------------
406 //------------------------------------------------------------------------------
407 bool gmIntersectLineSegmentsWithTol(const Pt3d& one1,
408  const Pt3d& one2,
409  const Pt3d& two1,
410  const Pt3d& two2,
411  double* xi /*=NULL*/,
412  double* yi /*=NULL*/,
413  double* zi1 /*=NULL*/,
414  double* zi2 /*=NULL*/,
415  double tol /*=0.0*/)
416 {
417  double minx1, minx2, miny1, miny2, maxx1, maxx2, maxy1, maxy2;
418  double dx1, dy1, dz1, dx2, dy2, dz2, lambda, mu, cross;
419  // do a trivial rejection
420  // RDJ - 4/20/2004 Do tests one at a time so if we fail on one we stop
421  minx1 = std::min(one1.x, one2.x);
422  maxx2 = std::max(two1.x, two2.x);
423  if (GT_TOL(minx1, maxx2, tol))
424  {
425  return false;
426  }
427  maxx1 = std::max(one1.x, one2.x);
428  minx2 = std::min(two1.x, two2.x);
429  if (LT_TOL(maxx1, minx2, tol))
430  {
431  return false;
432  }
433  miny1 = std::min(one1.y, one2.y);
434  maxy2 = std::max(two1.y, two2.y);
435  if (GT_TOL(miny1, maxy2, tol))
436  {
437  return false;
438  }
439  maxy1 = std::max(one1.y, one2.y);
440  miny2 = std::min(two1.y, two2.y);
441  if (LT_TOL(maxy1, miny2, tol))
442  {
443  return false;
444  }
445  // define the vectors
446  dx1 = one2.x - one1.x;
447  dy1 = one2.y - one1.y;
448  dz1 = one2.z - one1.z;
449  dx2 = two2.x - two1.x;
450  dy2 = two2.y - two1.y;
451  dz2 = two2.z - two1.z;
452  /* see if lines are parallel */
453  cross = (dx1 * dy2) -
454  (dy1 * dx2); // dy1/dx1 = dy2/dx2 lines have same slope
455  // This checks the line slopes. If equal then the lines are parallel or
456  // the same line. We assume they are parallel and exit.
457  // But really theta should be computed using vector magnitudes
458  if (EQ_TOL(cross, 0.0, tol))
459  return false;
460  // compute the value of lambda
461  lambda = (dy2 * (two1.x - one1.x) + dx2 * (one1.y - two1.y)) / cross;
462 
463  // There is some question as to what effect the tolerance should have on this
464  // function. It was decided that in the case where the segments do not
465  // intersect, but the minimum distance between the two segments is within
466  // tolerance, then the location of the minimum will be returned as the
467  // point of intersection.
468  // If the point of intersection is off the end of the first segment, then
469  // set it to be at the end and set the checkDistance flag.
470  bool checkDistance(false);
471  if (lambda < 0.0)
472  {
473  lambda = 0.0;
474  checkDistance = true;
475  }
476  else if (lambda > 1.0)
477  {
478  lambda = 1.0;
479  checkDistance = true;
480  }
481  // The intersection is inside of segment 1, so we need to check segment 2
482  // Compute mu from lambda and the two parametric segment equations
483  // two1 + mu (two2 - two1) = one1 + lambda ( one2 - one1)
484  // so
485  // mu = (one1 + lambda (one2-one1) - two1)/(two2-two1)
486  if (fabs(dx2) > fabs(dy2))
487  mu = (one1.x + lambda * dx1 - two1.x) / dx2;
488  else
489  mu = (one1.y + lambda * dy1 - two1.y) / dy2;
490 
491  // If the point of intersection is off the end of the second segment, then
492  // set it to be at the end and set the checkDistance flag.
493  if (mu < 0.0)
494  {
495  mu = 0.0;
496  checkDistance = true;
497  }
498  else if (mu > 1.0)
499  {
500  mu = 1.0;
501  checkDistance = true;
502  }
503  // if checkDistance flag is true check distance between mu/lambda
504  // positions (nearest points). If it is not within tolerance, return false.
505  Pt3d lambdapos = one1 + Pt3d(lambda * dx1, lambda * dy1, lambda * dz1);
506  Pt3d mupos = two1 + Pt3d(mu * dx2, mu * dy2, mu * dz2);
507  if (checkDistance)
508  {
509  if (!gmEqualPointsXY(lambdapos, mupos, tol))
510  {
511  return false;
512  }
513  }
514 
515  if (gmColinearWithTol(one1, one2, lambdapos, tol / 2) &&
516  gmColinearWithTol(two1, two2, lambdapos, tol / 2))
517  {
518  if (xi)
519  {
520  *xi = lambdapos.x;
521  }
522  if (yi)
523  {
524  *yi = lambdapos.y;
525  }
526  }
527  else
528  {
529  if (xi)
530  {
531  *xi = mupos.x;
532  }
533  if (yi)
534  {
535  *yi = mupos.y;
536  }
537  }
538  if (zi2)
539  *zi2 = mupos.z;
540  if (zi1)
541  *zi1 = lambdapos.z;
542  return true;
543 } // gmIntersectLineSegmentsWithTol
544 //------------------------------------------------------------------------------
550 //------------------------------------------------------------------------------
551 bool gmCounterClockwiseTri(const Pt3d& vtx0, const Pt3d& vtx1, const Pt3d& vtx2)
552 {
553  double triarea1 = trArea(vtx0, vtx1, vtx2);
554  return (triarea1 > 0.0);
555 } // gmCounterClockwiseTri
556 //------------------------------------------------------------------------------
565 //------------------------------------------------------------------------------
566 double gmCross2D(const double& dx1, const double& dy1, const double& dx2, const double& dy2)
567 {
568  return (dx1 * dy2) - (dx2 * dy1);
569 } // gmCross2D
570 //------------------------------------------------------------------------------
578 //------------------------------------------------------------------------------
579 double gmAngleBetween2DVectors(double dxp, double dyp, double dxn, double dyn)
580 {
581  double magn, magp;
582 
583  magn = sqrt(sqr(dxn) + sqr(dyn));
584  magp = sqrt(sqr(dxp) + sqr(dyp));
585  return gmAngleBetween2DVectors(dxp, dyp, dxn, dyn, magn, magp);
586 } // gmAngleBetween2DVectors
587 //----- OVERLOAD ---------------------------------------------------------------
597 //----- OVERLOAD ---------------------------------------------------------------
598 double gmAngleBetween2DVectors(double dxp,
599  double dyp,
600  double dxn,
601  double dyn,
602  double a_magn,
603  double a_magp)
604 {
605  double theangle, cosign;
606 
607  if (a_magn == 0.0 || a_magp == 0.0)
608  return (0.0);
609  cosign = (dxn * dxp + dyn * dyp) / (a_magn * a_magp);
610  if (cosign > 1.0)
611  cosign = 1.0;
612  if (cosign < -1.0)
613  cosign = -1.0;
614  theangle = acos(cosign);
615  if (theangle == 0.0)
616  {
617  if ((dxn * dxp) + (dyn * dyp) < 0.0)
618  theangle = XM_PI;
619  }
620  else if (gmCross2D(dxp, dyp, dxn, dyn) < 0.0)
621  theangle = 2 * XM_PI - theangle;
622  return theangle;
623 } // gmAngleBetween2DVectors
624 //------------------------------------------------------------------------------
630 //------------------------------------------------------------------------------
631 double gmAngleBetweenEdges(const Pt3d& p1, const Pt3d& p2, const Pt3d& p3)
632 {
633  double dxp, dyp, dxn, dyn;
634  /* compute the vectors */
635  dxp = p1.x - p2.x;
636  dyp = p1.y - p2.y;
637  dxn = p3.x - p2.x;
638  dyn = p3.y - p2.y;
639  return gmAngleBetween2DVectors(dxp, dyp, dxn, dyn);
640 } // gmAngleBetweenEdges
641 //------------------------------------------------------------------------------
647 //------------------------------------------------------------------------------
648 double gmAngleBetweenEdges(const Pt2d& p1, const Pt2d& p2, const Pt2d& p3)
649 {
650  double dxp, dyp, dxn, dyn;
651  /* compute the vectors */
652  dxp = p1.x - p2.x;
653  dyp = p1.y - p2.y;
654  dxn = p3.x - p2.x;
655  dyn = p3.y - p2.y;
656  return gmAngleBetween2DVectors(dxp, dyp, dxn, dyn);
657 } // gmAngleBetweenEdges
658 //------------------------------------------------------------------------------
672 //------------------------------------------------------------------------------
673 double gmComputeDeviationInDirection(const Pt3d& a_p0, const Pt3d& a_p1, const Pt3d& a_p2)
674 {
675  double x1, y1, x2, y2, m, cosine_dev;
676  x1 = a_p1.x - a_p0.x;
677  y1 = a_p1.y - a_p0.y;
678  x2 = a_p2.x - a_p1.x;
679  y2 = a_p2.y - a_p1.y;
680  m = sqrt(x1 * x1 + y1 * y1) * sqrt(x2 * x2 + y2 * y2);
681  if (m > 0.0)
682  {
683  cosine_dev = Mmin((x1 * x2 + y1 * y2) / (m), 1.0);
684  cosine_dev = Mmax(cosine_dev, -1.0);
685  return acos(cosine_dev);
686  }
687  else
688  {
689  return XM_PI;
690  }
691 } // gmComputeDeviationInDirection
692 //------------------------------------------------------------------------------
707 //------------------------------------------------------------------------------
708 bool gmOnLineWithTol(const Pt3d& p1,
709  const Pt3d& p2,
710  const double x,
711  const double y,
712  const double tol)
713 {
714  // compute vector components
715  double dx = p2.x - p1.x;
716  double dy = p2.y - p1.y;
717  double mag = sqrt(sqr(dx) + sqr(dy));
718  // check for extremely small segment
719  if (mag <= tol)
720  return gmEqualPointsXY(p1.x, p1.y, x, y);
721  else
722  {
723  double a = -dy / mag;
724  double b = dx / mag;
725  double c = -a * p2.x - b * p2.y;
726  // compute distance from line to (x,y)
727  double d = a * x + b * y + c;
728  return fabs(d) <= tol;
729  }
730 } // gmOnLineWithTol
731 //------------------------------------------------------------------------------
742 //------------------------------------------------------------------------------
743 bool gmOnLineAndBetweenEndpointsWithTol(const Pt3d& a_pt1,
744  const Pt3d& a_pt2,
745  const double a_x,
746  const double a_y,
747  double a_tol)
748 {
749  if ((Mmin(a_pt1.x, a_pt2.x) - a_tol <= a_x && Mmax(a_pt1.x, a_pt2.x) + a_tol >= a_x) &&
750  (Mmin(a_pt1.y, a_pt2.y) - a_tol <= a_y && Mmax(a_pt1.y, a_pt2.y) + a_tol >= a_y))
751  return gmOnLineWithTol(a_pt1, a_pt2, a_x, a_y, a_tol) == true;
752  else
753  return false;
754 } // gmOnLineAndBetweenEndpointsWithTol
755 //------------------------------------------------------------------------------
762 //------------------------------------------------------------------------------
763 void gmAddToExtents(const Pt3d& a_pt, Pt3d& a_min, Pt3d& a_max)
764 {
765  a_min.x = std::min(a_pt.x, a_min.x);
766  a_min.y = std::min(a_pt.y, a_min.y);
767  a_min.z = std::min(a_pt.z, a_min.z);
768  a_max.x = std::max(a_pt.x, a_max.x);
769  a_max.y = std::max(a_pt.y, a_max.y);
770  a_max.z = std::max(a_pt.z, a_max.z);
771 } // gmAddToExtents
772 //------------------------------------------------------------------------------
777 //------------------------------------------------------------------------------
778 void gmAddToExtents(const Pt3d& a_pt, Pt2d& a_min, Pt2d& a_max)
779 {
780  a_min.x = std::min(a_pt.x, a_min.x);
781  a_min.y = std::min(a_pt.y, a_min.y);
782  a_max.x = std::max(a_pt.x, a_max.x);
783  a_max.y = std::max(a_pt.y, a_max.y);
784 } // gmAddToExtents
785 //------------------------------------------------------------------------------
790 //------------------------------------------------------------------------------
791 void gmAddToExtents(const Pt2d& a_pt, Pt3d& a_min, Pt3d& a_max)
792 {
793  a_min.x = std::min(a_pt.x, a_min.x);
794  a_min.y = std::min(a_pt.y, a_min.y);
795  a_max.x = std::max(a_pt.x, a_max.x);
796  a_max.y = std::max(a_pt.y, a_max.y);
797 } // gmAddToExtents
798 //------------------------------------------------------------------------------
807 //------------------------------------------------------------------------------
808 double gmXyDistance(double x1, double y1, double x2, double y2)
809 {
810  return sqrt(sqr(x1 - x2) + sqr(y1 - y2));
811 } // gmXyDistance
812 //------------------------------------------------------------------------------
817 //------------------------------------------------------------------------------
818 double gmXyDistance(const Pt3d& pt1, const Pt3d& pt2)
819 {
820  return sqrt(sqr(pt1.x - pt2.x) + sqr(pt1.y - pt2.y));
821 } // gmXyDistance
822 //------------------------------------------------------------------------------
838 //------------------------------------------------------------------------------
839 Turn_enum gmTurn(const Pt3d& a_v1, const Pt3d& a_v2, const Pt3d& a_v3, double a_angtol)
840 {
841  // compute sin T = (v3-v2)x(v1-v2)/(d12*d32)
842  double dx32 = a_v3.x - a_v2.x;
843  double dy32 = a_v3.y - a_v2.y;
844  double dx12 = a_v1.x - a_v2.x;
845  double dy12 = a_v1.y - a_v2.y;
846  double d32 = sqrt(dx32 * dx32 + dy32 * dy32);
847  double d12 = sqrt(dx12 * dx12 + dy12 * dy12);
848  double mag = d12 * d32;
849  double sint = (dx32 * dy12 - dx12 * dy32) / mag;
850 
851  // for 99.999 > T > 0.1 deg - sin T > 0.0017
852  if (sint > a_angtol)
853  return TURN_LEFT;
854  else if (sint < -a_angtol)
855  return TURN_RIGHT;
856 
857  // compute cos T = (v3-v2)DOT(v1-v2)/(d12*d32)
858  // Near colinear case, Cosine should be near 1 or -1,
859  // -1 indicates 180 deg
860  double cost = (dx12 * dx32 + dy12 * dy32) / mag;
861  if (cost < 0.0)
862  return TURN_COLINEAR_180;
863  return TURN_COLINEAR_0;
864 } // gmTurn
865 //------------------------------------------------------------------------------
870 //------------------------------------------------------------------------------
871 Pt3d gmComputeCentroid(const VecPt3d& a_points)
872 {
873  Pt3d centroid;
874  size_t size = a_points.size();
875  for (size_t i = 0; i < size; ++i)
876  centroid += a_points[i];
877  return (centroid / (double)size);
878 } // gmComputeCentroid
879 //------------------------------------------------------------------------------
883 //------------------------------------------------------------------------------
884 Pt3d gmComputePolygonCentroid(const VecPt3d& pts)
885 {
886  Pt3d centroid;
887  if (pts.empty())
888  return centroid;
889  // get offset to use in calculation below to fix precision issues
890  double xMax = XM_DBL_LOWEST, yMax = XM_DBL_LOWEST, xMin = XM_DBL_HIGHEST, yMin = XM_DBL_HIGHEST;
891  size_t i = 0;
892  for (i = 0; i < pts.size(); ++i)
893  {
894  double x = pts[i].x;
895  double y = pts[i].y;
896  xMax = (x > xMax) ? x : xMax;
897  yMax = (y > yMax) ? y : yMax;
898  xMin = (x < xMin) ? x : xMin;
899  yMin = (y < yMin) ? y : yMin;
900  }
901  double xOffset = (xMax + xMin) / 2.0;
902  double yOffset = (yMax + yMin) / 2.0;
903  // For all vertices except last
904  double signedArea = 0.0;
905  for (i = 0; i < pts.size() - 1; ++i)
906  {
907  double x0 = pts[i].x - xOffset;
908  double y0 = pts[i].y - yOffset;
909  double x1 = pts[i + 1].x - xOffset;
910  double y1 = pts[i + 1].y - yOffset;
911  double a = x0 * y1 - x1 * y0;
912  signedArea += a;
913  centroid.x += (x0 + x1) * a;
914  centroid.y += (y0 + y1) * a;
915  }
916  // Do last vertex
917  double x0 = pts[i].x - xOffset;
918  double y0 = pts[i].y - yOffset;
919  double x1 = pts[0].x - xOffset;
920  double y1 = pts[0].y - yOffset;
921  double a = x0 * y1 - x1 * y0;
922  signedArea += a;
923  centroid.x += (x0 + x1) * a;
924  centroid.y += (y0 + y1) * a;
925  signedArea *= 0.5;
926  centroid.x /= (6.0 * signedArea);
927  centroid.y /= (6.0 * signedArea);
928  centroid.x += xOffset;
929  centroid.y += yOffset;
930  return centroid;
931 } // gmComputePolygonCentroid
932 //------------------------------------------------------------------------------
941 //------------------------------------------------------------------------------
942 bool gmLinesIntersect(const Pt3d& one1, const Pt3d& one2, const Pt3d& two1, const Pt3d& two2)
943 {
944  // do a trivial rejection
945  double min = std::min(one1.x, one2.x);
946  double max = std::max(two1.x, two2.x);
947  const double tol(XM_ZERO_TOL);
948  if (GT_TOL(min, max, tol))
949  return false;
950  max = std::max(one1.x, one2.x);
951  min = std::min(two1.x, two2.x);
952  if (LT_TOL(max, min, tol))
953  return false;
954  min = std::min(one1.y, one2.y);
955  max = std::max(two1.y, two2.y);
956  if (GT_TOL(min, max, tol))
957  return false;
958  max = std::max(one1.y, one2.y);
959  min = std::min(two1.y, two2.y);
960  if (LT_TOL(max, min, tol))
961  return false;
962  // define the vectors
963  Pt2d d1(one2 - one1), d2(two2 - two1);
964  // see if lines are parallel
965  // SMS bug fix: replaced tolerance check with a non-zero check.
966  // Sometimes we have very small tolerances
967  double cross = (d2.x * d1.y) - (d2.y * d1.x);
968  if (cross == 0.0)
969  return false;
970  // if the lines intersect s should be between 0.0 and 1.0
971  double s = ((d1.x * (two1.y - one1.y)) + (d1.y * (one1.x - two1.x))) / cross;
972  if (LT_TOL(s, 0.0, tol) || GT_TOL(s, 1.0, tol))
973  return false;
974  // if the lines intersect t should be between -1.0 and 0.0
975  double t = ((d2.x * (one1.y - two1.y)) + (d2.y * (two1.x - one1.x))) / cross;
976  if (LT_TOL(t, -1.0, tol) || GT_TOL(t, 0.0, tol))
977  return false;
978  return true;
979 } // gmLinesIntersect
980 //------------------------------------------------------------------------------
993 //------------------------------------------------------------------------------
994 template <typename T>
995 int gmPointInPolygon2D(const T* a_verts,
996  const size_t a_num,
997  const double a_x,
998  const double a_y,
999  const double a_tol)
1000 {
1001  if (a_verts && a_num)
1002  {
1003  int nmleft = 0, nmrght = 0;
1004  double dy2 = fabs(a_verts[0].y - a_y);
1005  double diff;
1006  for (size_t i = 0; i < a_num; i++)
1007  {
1008  size_t i2 = 0;
1009  if (i != a_num - 1)
1010  i2 = i + 1;
1011  double dy1 = dy2;
1012  dy2 = fabs(a_verts[i2].y - a_y);
1013  if (dy1 <= a_tol)
1014  {
1015  if (dy2 <= a_tol)
1016  {
1017  // Case #1 - edge is on the same y value as the point
1018  // - see if the point is on the edge
1019  if (((a_verts[i].x >= a_x) && (a_verts[i2].x <= a_x)) ||
1020  ((a_verts[i].x <= a_x) && (a_verts[i2].x >= a_x)))
1021  return (0);
1022  }
1023  else
1024  {
1025  // Case #2 - first vertex has same y value as point
1026  // - see if the point actually coincides with
1027  // the ith vertex
1028  diff = a_verts[i].x - a_x;
1029  if (fabs(diff) <= a_tol)
1030  return (0);
1031  else if (a_verts[i].y < a_verts[i2].y)
1032  {
1033  // edge ascends, classify as right or left
1034  if (diff > 0)
1035  ++nmrght;
1036  else
1037  ++nmleft;
1038  }
1039  }
1040  }
1041  else if (dy2 <= a_tol)
1042  {
1043  // Case #3 - 2nd vertex of edge lies on dividing plane
1044  // - see if the point actually coincides with
1045  // the i2-th vertex
1046  diff = a_verts[i2].x - a_x;
1047  if (fabs(diff) <= a_tol)
1048  return (0);
1049  else if (a_verts[i2].y < a_verts[i].y)
1050  {
1051  // if edge descends classify as right or left
1052  if (diff > 0)
1053  ++nmrght;
1054  else
1055  ++nmleft;
1056  }
1057  }
1058  else if (((a_verts[i].y < a_y) && (a_verts[i2].y > a_y)) ||
1059  ((a_verts[i].y > a_y) && (a_verts[i2].y < a_y)))
1060  {
1061  // Case #4 - edge cleanly intersects dividing plane
1062  // - flag edge as left, right or on edge
1063  double val = a_verts[i].x + (a_verts[i2].x - a_verts[i].x) * (a_y - a_verts[i].y) /
1064  (a_verts[i2].y - a_verts[i].y);
1065  diff = val - a_x;
1066  if (fabs(diff) <= a_tol)
1067  return (0);
1068  else if (diff > 0)
1069  ++nmrght;
1070  else
1071  ++nmleft;
1072  }
1073  }
1074  nmleft = nmleft % 2;
1075  nmrght = nmrght % 2;
1076  if (nmleft != nmrght)
1077  return (-1); // this should never happen actually
1078  else if (nmleft == 1)
1079  return (1);
1080  else
1081  return (-1);
1082  }
1083  return (-1);
1084 } // gmPointInPolygon2D
1085 //------------------------------------------------------------------------------
1092 //------------------------------------------------------------------------------
1093 int gmPointInPolygon2D(const Pt3d* a_verts, size_t a_num, double a_x, double a_y)
1094 {
1095  return gmPointInPolygon2D(a_verts, a_num, a_x, a_y, gmXyTol());
1096 } // gmPointInPolygon2D
1097 //------------------------------------------------------------------------------
1103 //------------------------------------------------------------------------------
1104 int gmPointInPolygon2D(const Pt3d* a_verts, size_t a_num, Pt3d a_pt)
1105 {
1106  return gmPointInPolygon2D(a_verts, a_num, a_pt.x, a_pt.y, gmXyTol());
1107 } // gmPointInPolygon2D
1108 //------------------------------------------------------------------------------
1114 //------------------------------------------------------------------------------
1115 int gmPointInPolygon2D(const Pt2i* a_verts, size_t a_num, Pt2d a_pt)
1116 {
1117  return gmPointInPolygon2D(a_verts, a_num, a_pt.x, a_pt.y, gmXyTol());
1118 } // gmPointInPolygon2D
1119 //------------------------------------------------------------------------------
1125 //------------------------------------------------------------------------------
1126 int gmPointInPolygon2D(const Pt2i* a_verts, size_t a_num, Pt3d a_pt)
1127 {
1128  return gmPointInPolygon2D(a_verts, a_num, a_pt.x, a_pt.y, gmXyTol());
1129 } // gmPointInPolygon2D
1130 //------------------------------------------------------------------------------
1136 //------------------------------------------------------------------------------
1137 int gmPointInPolygon2D(const Pt2i* a_verts, size_t a_num, Pt2i a_pt)
1138 {
1139  return gmPointInPolygon2D(a_verts, a_num, a_pt.x, a_pt.y, gmXyTol());
1140 } // gmPointInPolygon2D
1141 //------------------------------------------------------------------------------
1146 //------------------------------------------------------------------------------
1147 int gmPointInPolygon2D(const VecPt3d& a_verts, const Pt3d& a_pt)
1148 {
1149  return gmPointInPolygon2D(&a_verts[0], a_verts.size(), a_pt.x, a_pt.y, gmXyTol());
1150 } // gmPointInPolygon2D
1151 //------------------------------------------------------------------------------
1159 //------------------------------------------------------------------------------
1160 template int gmPointInPolygon2D<Pt3d>(const Pt3d* a_verts,
1161  const size_t a_num,
1162  const double a_x,
1163  const double a_y,
1164  const double a_tol);
1165 //------------------------------------------------------------------------------
1173 //------------------------------------------------------------------------------
1174 template int gmPointInPolygon2D<Pt2d>(const Pt2d* a_verts,
1175  const size_t a_num,
1176  const double a_x,
1177  const double a_y,
1178  const double a_tol);
1179 //------------------------------------------------------------------------------
1187 //------------------------------------------------------------------------------
1188 template int gmPointInPolygon2D<Pt2i>(const Pt2i* a_verts,
1189  const size_t a_num,
1190  const double a_x,
1191  const double a_y,
1192  const double a_tol);
1193 //------------------------------------------------------------------------------
1199 //------------------------------------------------------------------------------
1200 double gmComputeXyTol(const Pt3d& a_mn, const Pt3d& a_mx)
1201 {
1202  double d = gmXyDistance(a_mn, a_mx);
1203  double const kFactor = 1e-9;
1204  double xytol = d * kFactor;
1205  if (xytol < kFactor)
1206  {
1207  xytol = kFactor;
1208  }
1209  return xytol;
1210 } // gmComputeXyTol
1211 //------------------------------------------------------------------------------
1216 //------------------------------------------------------------------------------
1217 double gmXyTol(bool a_set /*false*/, double a_value /*1e-9*/)
1218 {
1219  static double xytol = a_value;
1220  if (a_set)
1221  xytol = a_value;
1222  return xytol;
1223 } // gmXyTol
1224 //------------------------------------------------------------------------------
1229 //------------------------------------------------------------------------------
1230 double gmZTol(bool a_set, double a_value)
1231 {
1232  static double ztol = a_value;
1233  if (a_set)
1234  ztol = a_value;
1235  return ztol;
1236 } // gmZTol
1237 //------------------------------------------------------------------------------
1245 //------------------------------------------------------------------------------
1246 bool gmEqualPointsXY(double x1, double y1, double x2, double y2, double tolerance)
1247 {
1248  double dx = fabs(x1 - x2);
1249  double dy = fabs(y1 - y2);
1250  if (dx > tolerance || dy > tolerance)
1251  return false;
1252  else if (sqrt(dx * dx + dy * dy) <= tolerance)
1253  return true;
1254  else
1255  return false;
1256 } // gmEqualPointsXY
1257 //------------------------------------------------------------------------------
1264 //------------------------------------------------------------------------------
1265 bool gmEqualPointsXY(double x1, double y1, double x2, double y2)
1266 {
1267  return gmEqualPointsXY(x1, y1, x2, y2, gmXyTol());
1268 } // gmEqualPointsXY
1269 //------------------------------------------------------------------------------
1275 //------------------------------------------------------------------------------
1276 bool gmEqualPointsXY(const Pt2d& a_pt1, const Pt2d& a_pt2, double tol)
1277 {
1278  return gmEqualPointsXY(a_pt1.x, a_pt1.y, a_pt2.x, a_pt2.y, tol);
1279 } // gmEqualPointsXY
1280 //------------------------------------------------------------------------------
1286 //------------------------------------------------------------------------------
1287 bool gmEqualPointsXY(const Pt3d& a_pt1, const Pt3d& a_pt2, double tol)
1288 {
1289  return gmEqualPointsXY(a_pt1.x, a_pt1.y, a_pt2.x, a_pt2.y, tol);
1290 } // gmEqualPointsXY
1291 //------------------------------------------------------------------------------
1296 //------------------------------------------------------------------------------
1297 bool gmEqualPointsXY(const Pt2i& point1, const Pt2i& point2)
1298 {
1299  if (point1.x == point2.x && point1.y == point2.y)
1300  return true;
1301  return false;
1302 } // gmEqualPointsXY
1303 //------------------------------------------------------------------------------
1313 //------------------------------------------------------------------------------
1314 bool gmEqualPointsXYZ(double x1,
1315  double y1,
1316  double z1,
1317  double x2,
1318  double y2,
1319  double z2,
1320  double tolerance)
1321 {
1322  if ((fabs(x1 - x2) <= tolerance) && (fabs(y1 - y2) <= tolerance) && (fabs(z1 - z2) <= tolerance))
1323  return true;
1324  return false;
1325 } // gmEqualPointsXYZ
1326 //------------------------------------------------------------------------------
1335 //------------------------------------------------------------------------------
1336 bool gmEqualPointsXYZ(double x1, double y1, double z1, double x2, double y2, double z2)
1337 {
1338  return gmEqualPointsXYZ(x1, y1, z1, x2, y2, z2, gmXyTol());
1339 } // gmEqualPointsXYZ
1340 //------------------------------------------------------------------------------
1346 //------------------------------------------------------------------------------
1347 bool gmEqualPointsXYZ(const Pt3d& pt1, const Pt3d& pt2, double tol)
1348 {
1349  return gmEqualPointsXYZ(pt1.x, pt1.y, pt1.z, pt2.x, pt2.y, pt2.z, tol);
1350 } // gmEqualPointsXYZ
1351 //------------------------------------------------------------------------------
1360 //------------------------------------------------------------------------------
1361 bool gmPointInTriangleWithTol(const Pt3d* p1,
1362  const Pt3d* p2,
1363  const Pt3d* p3,
1364  double x,
1365  double y,
1366  double tol)
1367 {
1368  if (gmInsideOrOnLineWithTol(p1, p2, p3, x, y, tol))
1369  if (gmInsideOrOnLineWithTol(p2, p3, p1, x, y, tol))
1370  if (gmInsideOrOnLineWithTol(p3, p1, p2, x, y, tol))
1371  return true;
1372  return false;
1373 } // gmPointInTriangleWithTol
1374 //------------------------------------------------------------------------------
1387 //------------------------------------------------------------------------------
1388 bool gmInsideOrOnLineWithTol(const Pt3d* p1,
1389  const Pt3d* p2,
1390  const Pt3d* inpoint,
1391  const double x,
1392  const double y,
1393  const double tol,
1394  double* dist /*= NULL */)
1395 {
1396  double dx, dy, a, b, c, mag;
1397  double d1, d2;
1398  /* make sure line is not a point */
1399  dx = p2->x - p1->x;
1400  dy = p2->y - p1->y;
1401  mag = sqrt(sqr(dx) + sqr(dy));
1402  if (mag <= tol)
1403  {
1404  return gmEqualPointsXY(p1->x, p1->y, x, y);
1405  }
1406  else
1407  {
1408  a = -dy / mag;
1409  b = dx / mag;
1410  c = -a * p1->x - b * p1->y;
1411  /* compute distance from line to (x,y) */
1412  d1 = a * x + b * y + c;
1413  /* compute distance from line to opposite p */
1414  d2 = a * inpoint->x + b * inpoint->y + c;
1415  // if inpoint is really on, we can't make a judgment on the point
1416  if (Mfabs(d2) <= tol)
1417  {
1418  // AKZ - Since "inpoint" is usually the other point in a triangle
1419  // this usually indicates a poorly formed triangle.
1420  // You should determine how this triangle was formed.
1421  // XM_ASSERT(0);
1422  return false;
1423  }
1424  /* return true if point on line */
1425  if (Mfabs(d1) <= tol)
1426  {
1427  if (dist)
1428  {
1429  if (d1 * d2 < 0.0)
1430  {
1431  *dist = Mfabs(d1);
1432  }
1433  else
1434  {
1435  *dist = -Mfabs(d1);
1436  }
1437  }
1438  return true;
1439  }
1440  /* return true if d1 and d2 have same sign */
1441  if ((d1 < 0.0) && (d2 < 0.0))
1442  {
1443  if (dist)
1444  {
1445  *dist = -Mfabs(d1);
1446  }
1447  return true;
1448  }
1449  else if ((d1 > 0.0) && (d2 > 0.0))
1450  {
1451  if (dist)
1452  {
1453  *dist = -Mfabs(d1);
1454  }
1455  return true;
1456  }
1457  else
1458  {
1459  if (dist)
1460  {
1461  *dist = Mfabs(d1);
1462  }
1463  return false;
1464  }
1465  }
1466 } // gmInsideOrOnLineWithTol
1467 //------------------------------------------------------------------------------
1474 //------------------------------------------------------------------------------
1475 double gmPolygonArea(const Pt3d* pts, size_t npoints)
1476 {
1477  size_t id;
1478  double area = 0.0;
1479 
1480  /* original method with precision errors
1481  for (id = 0; id < npoints; id++)
1482  {
1483  if (id != (npoints - 1))
1484  {
1485  area += (pts[id].x * pts[id + 1].y);
1486  area -= (pts[id].y * pts[id + 1].x);
1487  }
1488  else
1489  {
1490  area += (pts[id].x * pts[0].y);
1491  area -= (pts[id].y * pts[0].x);
1492  }
1493  }
1494  area /= 2.0;
1495  */
1496 
1497  // AKZ 2/15/2018
1498  // I changed the implementation to translate the polygon
1499  // so that the first point is at the origin
1500  // Reduces round off error due to large coordinates
1501  // Reduces the number of computations because the first and last
1502  // computations in the loop would be 0.0
1503  VecDbl x, y;
1504  double x0 = pts[0].x;
1505  double y0 = pts[0].y;
1506  for (id = 1; id < npoints; id++)
1507  {
1508  x.push_back((pts[id].x - x0));
1509  y.push_back((pts[id].y - y0));
1510  }
1511  for (id = 0; id < npoints - 2; id++)
1512  {
1513  area += (x[id] * y[id + 1]);
1514  area -= (y[id] * x[id + 1]);
1515  }
1516  area /= 2.0;
1517 
1518  return (area);
1519 } // gmPolygonArea
1520 //------------------------------------------------------------------------------
1525 //------------------------------------------------------------------------------
1526 VecPt3d gmArrayToVecPt3d(double* a_array, int a_size)
1527 {
1528  VecPt3d v(a_size / 2);
1529  for (int i = 0; i < a_size; i += 2)
1530  {
1531  v[i / 2].x = a_array[i];
1532  v[i / 2].y = a_array[i + 1];
1533  }
1534  return v;
1535 } // gmArrayToVecPt3d
1536 //------------------------------------------------------------------------------
1541 //------------------------------------------------------------------------------
1542 void gmEnvelopeOfPts(const VecPt3d& a_pts, Pt3d& a_min, Pt3d& a_max)
1543 {
1544  a_min = a_max = Pt3d();
1545  XM_ENSURE_TRUE(!a_pts.empty());
1546  a_min = a_max = a_pts.front();
1547  for (size_t i = 0; i < a_pts.size(); ++i)
1548  {
1549  if (a_pts[i].x < a_min.x)
1550  a_min.x = a_pts[i].x;
1551  if (a_pts[i].y < a_min.y)
1552  a_min.y = a_pts[i].y;
1553  if (a_pts[i].z < a_min.z)
1554  a_min.z = a_pts[i].z;
1555  if (a_pts[i].x > a_max.x)
1556  a_max.x = a_pts[i].x;
1557  if (a_pts[i].y > a_max.y)
1558  a_max.y = a_pts[i].y;
1559  if (a_pts[i].z > a_max.z)
1560  a_max.z = a_pts[i].z;
1561  }
1562 } // gmEnvelopeOfPts
1563 //------------------------------------------------------------------------------
1570 //------------------------------------------------------------------------------
1571 void gmOrderPointsCounterclockwise(const VecPt3d& a_pts, VecInt& a_ccwOrder, int a_startindex)
1572 {
1573  int numnodes = (int)a_pts.size();
1574 
1575  // compute centroid of points
1576  Pt2d center;
1577  double sumx = 0.0;
1578  double sumy = 0.0;
1579  for (int i = 0; i < numnodes; i++)
1580  {
1581  sumx += a_pts[i].x;
1582  sumy += a_pts[i].y;
1583  }
1584 
1585  center.x = sumx / numnodes;
1586  center.y = sumy / numnodes;
1587 
1588  // compute polar angle for each node about the centroid of the element
1589  // along with the original point index
1590  std::vector<std::pair<float, int> > angleIndex(numnodes);
1591  for (int i = 0; i < numnodes; i++)
1592  {
1593  float angle = (float)(atan2(a_pts[i].y - center.y, a_pts[i].x - center.x));
1594  angleIndex[i] = std::pair<float, int>(angle, i);
1595  }
1596 
1597  std::stable_sort(angleIndex.begin(), angleIndex.end());
1598 
1599  // find where starting index is in sorted array
1600  auto startIter = angleIndex.begin();
1601  while (startIter != angleIndex.end() && startIter->second != a_startindex)
1602  {
1603  ++startIter;
1604  }
1605  if (startIter == angleIndex.end())
1606  startIter = angleIndex.begin();
1607 
1608  // assign the result preserving the order with a_startindex the first node
1609  a_ccwOrder.resize(numnodes, 0);
1610  size_t j = 0;
1611  for (auto iter = startIter; iter != angleIndex.end(); ++iter)
1612  {
1613  a_ccwOrder[j] = iter->second;
1614  j++;
1615  }
1616  for (auto iter = angleIndex.begin(); iter != startIter; ++iter)
1617  {
1618  a_ccwOrder[j] = iter->second;
1619  j++;
1620  }
1621 } // gmOrderPointsCounterclockwise
1622 //------------------------------------------------------------------------------
1625 //------------------------------------------------------------------------------
1626 void gmOrderPointsCounterclockwise(VecPt3d& a_pts)
1627 {
1628  if (a_pts.empty())
1629  return;
1630 
1631  VecInt ccwOrder;
1632  VecPt3d pts(a_pts);
1633  gmOrderPointsCounterclockwise(pts, ccwOrder);
1634  for (size_t i = 0; i < ccwOrder.size(); ++i)
1635  {
1636  a_pts[i] = pts[ccwOrder[i]];
1637  }
1638 } // gmOrderPointsCounterclockwise
1639 //------------------------------------------------------------------------------
1649 //------------------------------------------------------------------------------
1650 double gmFindClosestPtOnSegment(const Pt3d& a_pt1,
1651  const Pt3d& a_pt2,
1652  const Pt3d& a_pt,
1653  Pt3d& a_newpt,
1654  const double a_tol)
1655 {
1656  double t = gmPtDistanceAlongSegment(a_pt1, a_pt2, a_pt, a_tol);
1657  if (t < 0.0)
1658  {
1659  a_newpt.x = a_pt1.x;
1660  a_newpt.y = a_pt1.y;
1661  }
1662  else if (t > 1.0)
1663  {
1664  a_newpt.x = a_pt2.x;
1665  a_newpt.y = a_pt2.y;
1666  }
1667  else
1668  {
1669  double dx = a_pt2.x - a_pt1.x;
1670  double dy = a_pt2.y - a_pt1.y;
1671  a_newpt.x = a_pt1.x + dx * t;
1672  a_newpt.y = a_pt1.y + dy * t;
1673  }
1674  return t;
1675 } // gmFindClosestPtOnSegment
1676 //------------------------------------------------------------------------------
1684 //------------------------------------------------------------------------------
1685 double gmPtDistanceAlongSegment(const Pt3d& a_pt1,
1686  const Pt3d& a_pt2,
1687  const Pt3d& a_pt,
1688  const double a_tol)
1689 {
1690  double dx, dy, t;
1691 
1692  dx = a_pt2.x - a_pt1.x;
1693  dy = a_pt2.y - a_pt1.y;
1694 
1695  if ((dx == 0.0 && dy == 0.0) || (sqrt(dx * dx + dy * dy) <= a_tol))
1696  {
1697  t = -1.0;
1698  }
1699  else
1700  {
1701  t = ((a_pt.x - a_pt1.x) * dx + (a_pt.y - a_pt1.y) * dy) / (dx * dx + dy * dy);
1702  }
1703 
1704  return t;
1705 } // gmPtDistanceAlongSegment
1706 //------------------------------------------------------------------------------
1717 //------------------------------------------------------------------------------
1718 bool gmInsideOfLineWithTol(const Pt3d& a_vertex1,
1719  const Pt3d& a_vertex2,
1720  const Pt3d& a_oppositevertex,
1721  const double a_x,
1722  const double a_y,
1723  const double a_tol)
1724 {
1725  double dx = a_vertex2.x - a_vertex1.x;
1726  double dy = a_vertex2.y - a_vertex1.y;
1727  double mag = sqrt(sqr(dx) + sqr(dy));
1728  if (mag <= a_tol)
1729  {
1730  return gmEqualPointsXY(a_vertex1.x, a_vertex1.y, a_x, a_y);
1731  }
1732  else
1733  {
1734  double a = -dy / mag;
1735  double b = dx / mag;
1736  double c = -a * a_vertex1.x - b * a_vertex1.y;
1737  /* compute the distance from the line to the Point */
1738  double d1 = a * a_x + b * a_y + c;
1739  /* compute the distance from the line to the opposite vertex */
1740  double d2 = a * a_oppositevertex.x + b * a_oppositevertex.y + c;
1741  if (fabs(d1) <= a_tol)
1742  return false;
1743  else if ((d1 < 0.0) && (d2 < 0.0))
1744  return true;
1745  else if ((d1 > 0.0) && (d2 > 0.0))
1746  return true;
1747  else
1748  return false;
1749  }
1750 } // gmInsideOfLine
1751 //------------------------------------------------------------------------------
1756 //------------------------------------------------------------------------------
1757 void gmExtents2D(const VecPt3d& a_points, Pt2d& a_min, Pt2d& a_max)
1758 {
1759  XM_ENSURE_TRUE_VOID_NO_ASSERT(!a_points.empty());
1760 
1761  a_min.x = a_max.x = a_points.at(0).x;
1762  a_min.y = a_max.y = a_points.at(0).y;
1763  for (int i = 1; i < (int)a_points.size(); i++)
1764  {
1765  gmAddToExtents(a_points[i], a_min, a_max);
1766  }
1767 } // gmExtents2D
1768 //----- OVERLOAD ---------------------------------------------------------------
1773 //----- OVERLOAD ---------------------------------------------------------------
1774 void gmExtents2D(const VecPt3d& a_points, Pt3d& a_min, Pt3d& a_max)
1775 {
1776  XM_ENSURE_TRUE_VOID_NO_ASSERT(!a_points.empty());
1777 
1778  a_min.x = a_max.x = a_points.at(0).x;
1779  a_min.y = a_max.y = a_points.at(0).y;
1780  a_min.z = a_max.z = 0.0;
1781  for (int i = 1; i < (int)a_points.size(); i++)
1782  {
1783  gmAddToExtents((Pt2d)a_points[i], a_min, a_max);
1784  }
1785 } // gmExtents2D
1786 //------------------------------------------------------------------------------
1791 //------------------------------------------------------------------------------
1792 void gmExtents3D(const VecPt3d& a_points, Pt3d& a_min, Pt3d& a_max)
1793 {
1794  XM_ENSURE_TRUE_VOID_NO_ASSERT(!a_points.empty());
1795 
1796  a_min.x = a_max.x = a_points.at(0).x;
1797  a_min.y = a_max.y = a_points.at(0).y;
1798  a_min.z = a_max.z = a_points.at(0).z;
1799  for (int i = 1; i < (int)a_points.size(); i++)
1800  {
1801  gmAddToExtents(a_points[i], a_min, a_max);
1802  }
1803 } // gmExtents3D
1804 //------------------------------------------------------------------------------
1809 //------------------------------------------------------------------------------
1810 double gmPerpendicularAngle(const Pt3d& a_pt1, const Pt3d& a_pt2)
1811 {
1812  double hypot;
1813  double deltax, deltay, arad, theangle;
1814 
1815  deltax = a_pt1.x - a_pt2.x;
1816  deltay = a_pt1.y - a_pt2.y;
1817  hypot = sqrt(sqr(deltax) + sqr(deltay));
1818  arad = deltax / hypot;
1819  if (arad > .9999)
1820  arad = 1.0;
1821  else if (arad < -.9999)
1822  arad = -1.0;
1823  if (deltay >= 0.0)
1824  theangle = acos(arad);
1825  else
1826  theangle = 2 * XM_PI - acos(arad);
1827  return (theangle - (XM_PI / 2));
1828 } // gmPerpendicularAngle
1829 //------------------------------------------------------------------------------
1836 //------------------------------------------------------------------------------
1837 double gmBisectingAngle(const Pt3d& a_p1, const Pt3d& a_p2, const Pt3d& a_p3)
1838 {
1839  double dxn, dyn, dxp, dyp, magn, magp, angletoedge1, theanglebetween;
1840  double cosign;
1841 
1842  dxp = a_p1.x - a_p2.x;
1843  dyp = a_p1.y - a_p2.y;
1844  dxn = a_p3.x - a_p2.x;
1845  dyn = a_p3.y - a_p2.y;
1846  angletoedge1 = atan2(dyp, dxp);
1847  magn = sqrt(sqr(dxn) + sqr(dyn));
1848  magp = sqrt(sqr(dxp) + sqr(dyp));
1849  cosign = (dxn * dxp + dyn * dyp) / (magn * magp);
1850  if (cosign > .99999)
1851  cosign = 1.0;
1852  if (cosign < -.99999)
1853  cosign = -1.0;
1854  theanglebetween = acos(cosign);
1855  if (theanglebetween == 0.0)
1856  {
1857  if ((dxn * dxp) + (dyn * dyp) < 0.0)
1858  theanglebetween = XM_PI;
1859  }
1860  else if (gmCross2D(dxp, dyp, dxn, dyn) < 0.0)
1861  theanglebetween = 2 * XM_PI - theanglebetween;
1862  return angletoedge1 + theanglebetween / 2.0;
1863 } // gmBisectingAngle
1864 //------------------------------------------------------------------------------
1876 //------------------------------------------------------------------------------
1877 void gmComponentMagnitudes(double* a_x, double* a_y, double* a_mag, double* a_dir, bool a_tomagdir)
1878 {
1879  double rads;
1880 
1881  if (a_tomagdir)
1882  { // convert (x,y) to (mag,dir)
1883  if (fabs(*a_x) < XM_ZERO_TOL && fabs(*a_y) < XM_ZERO_TOL)
1884  {
1885  *a_mag = 0.0;
1886  *a_dir = 0.0;
1887  }
1888  else
1889  {
1890  if (*a_x == 0)
1891  *a_x = XM_ZERO_TOL;
1892  *a_mag = sqrt(sqr(*a_x) + sqr(*a_y));
1893  *a_dir = (atan(*a_y / *a_x)) * (180 / XM_PI);
1894  if (*a_x < 0)
1895  (*a_dir) += 180;
1896  if (*a_dir < 0)
1897  (*a_dir) += 360;
1898  }
1899  }
1900  else
1901  { // convert (mag,dir) to (x,y)
1902  rads = *a_dir * (XM_PI / 180);
1903  *a_x = cos(rads) * *a_mag;
1904  *a_y = sin(rads) * *a_mag;
1905  if (fabs(*a_x) < XM_ZERO_TOL)
1906  *a_x = 0;
1907  if (fabs(*a_y) < XM_ZERO_TOL)
1908  *a_y = 0;
1909  }
1910 } // gmComponentMagnitudes
1911 //------------------------------------------------------------------------------
1917 //------------------------------------------------------------------------------
1918 Pt3d gmCreateVector(const Pt3d& a_p1, const Pt3d& a_p2)
1919 {
1920  Pt3d vector;
1921  vector.x = a_p2.x - a_p1.x;
1922  vector.y = a_p2.y - a_p1.y;
1923  vector.z = a_p2.z - a_p1.z;
1924  return vector;
1925 } // gmCreateVector
1926 //------------------------------------------------------------------------------
1932 //------------------------------------------------------------------------------
1933 double gmConvertAngleToBetween0And360(double a_angle, bool a_InDegrees /*= true*/
1934 )
1935 {
1936  double ang = a_angle;
1937 
1938  if (!a_InDegrees)
1939  {
1940  ang *= (180.0 / XM_PI);
1941  }
1942 
1943 #if BOOST_OS_WINDOWS
1944  while (LT_TOL(ang, 0.0, DBL_EPSILON) && _finite(ang))
1945  {
1946 #else
1947  while (LT_TOL(ang, 0.0, DBL_EPSILON) && std::isfinite(ang))
1948  {
1949 #endif
1950  ang += 360.0;
1951  }
1952 #if BOOST_OS_WINDOWS
1953  while (GTEQ_TOL(ang, 360.0, DBL_EPSILON) && _finite(ang))
1954  {
1955 #else
1956  while (GTEQ_TOL(ang, 360.0, DBL_EPSILON) && std::isfinite(ang))
1957  {
1958 #endif
1959  ang -= 360.0;
1960  }
1961 
1962  return ang;
1963 } // gmConvertAngleToBetween0And360
1964 //------------------------------------------------------------------------------
1969 //------------------------------------------------------------------------------
1970 void gmCross3D(const Pt3d& a_vec1, const Pt3d& a_vec2, Pt3d* a_vec3)
1971 {
1972  a_vec3->x = a_vec1.y * a_vec2.z - a_vec1.z * a_vec2.y;
1973  a_vec3->y = a_vec1.z * a_vec2.x - a_vec1.x * a_vec2.z;
1974  a_vec3->z = a_vec1.x * a_vec2.y - a_vec1.y * a_vec2.x;
1975 } // gmCross3D
1976 //------------------------------------------------------------------------------
1983 //------------------------------------------------------------------------------
1984 inline double gmDot3D(const Pt3d& a_vec1, const Pt3d& a_vec2)
1985 {
1986  return (a_vec1.x * a_vec2.x + a_vec1.y * a_vec2.y + a_vec1.z * a_vec2.z);
1987 } // gmDot3D
1988 //------------------------------------------------------------------------------
2002 //
2003 // Copyright 2001, softSurfer (www.softsurfer.com)
2004 // This code may be freely used and modified for any purpose providing that
2005 // this copyright notice is included with it. SoftSurfer makes no warranty
2006 // for this code, and cannot be held liable for any real or imagined damage
2007 // resulting from its use. Users of this code must verify correctness for
2008 // their application.
2009 // http://geometryalgorithms.com/Archive/algorithm_0105/algorithm_0105.htm
2010 //------------------------------------------------------------------------------
2011 int gmIntersectTriangleAndLineSegment(const Pt3d& a_pt1,
2012  const Pt3d& a_pt2,
2013  const Pt3d& a_t0,
2014  const Pt3d& a_t1,
2015  const Pt3d& a_t2,
2016  Pt3d& a_IntersectPt)
2017 {
2018  double a, b, r;
2019  Pt3d dir, n, NullVector(0.0, 0.0, 0.0), u, v, w, w0;
2020 
2021  // get the triangle edge vectors and plane normal
2022  u = a_t1 - a_t0;
2023  v = a_t2 - a_t0;
2024  gmCross3D(u, v, &n);
2025 
2026  // check if the triangle is degenerate
2027  if (n == NullVector)
2028  {
2029  return -1;
2030  }
2031 
2032  // get the ray direction vector
2033  dir = a_pt2 - a_pt1;
2034  w0 = a_pt1 - a_t0;
2035  a = -gmDot3D(n, w0);
2036  b = gmDot3D(n, dir);
2037 
2038  // see if ray is parallel to the triangle
2039  if (fabs(b) < XM_ZERO_TOL)
2040  {
2041  // see if ray lies in triangle plane
2042  if (a == 0)
2043  {
2044  return 2;
2045  }
2046  // else ray is disjoint from the triangle plane
2047  else
2048  {
2049  return 0;
2050  }
2051  }
2052 
2053  // get the intersection point or ray with triangle plane
2054  r = a / b;
2055 
2056  // see if there is an intersection
2057  // if (r < 0.0 || r > 1.0) {
2058  if (r < -FLT_EPSILON || r > 1.0 + FLT_EPSILON)
2059  {
2060  return 0;
2061  }
2062 
2063  // intersect point of ray and plane
2064  a_IntersectPt = a_pt1 + dir * r;
2065 
2066  // see if the intersection is inside of the triangle
2067  double D, uu, uv, vv, wu, wv;
2068 
2069  uu = gmDot3D(u, u);
2070  uv = gmDot3D(u, v);
2071  vv = gmDot3D(v, v);
2072  w = a_IntersectPt - a_t0;
2073  wu = gmDot3D(w, u);
2074  wv = gmDot3D(w, v);
2075  D = uv * uv - uu * vv;
2076 
2077  // get the test parametric coords
2078  double s, t;
2079 
2080  s = (uv * wv - vv * wu) / D;
2081  if (s < 0.0 || s > 1.0)
2082  {
2083  // the intersect point is outside the triangle
2084  return 0;
2085  }
2086  t = (uv * wu - uu * wv) / D;
2087  if (t < 0.0 || (s + t) > 1.0)
2088  {
2089  // the intersect point is outside the triangle
2090  return 0;
2091  }
2092 
2093  // the intersect point is inside the triangle
2094  return 1;
2095 } // gmIntersectTriangleAndLineSegment
2096 /***********************************************************************
2097 * FUNCTION gm2DDistanceToLineWithTol
2098 * PURPOSE return the xy distance from (x,y) to line through (pt1, pt2)
2099 * NOTES this doesn't return fabs of the dist!! Differs from
2100 gm2DDistanceToLineSegmentWithTol in that the line is not
2101 treated like a line segment and the point on the line closest
2102 to xy may be outside the 2 points defining the line.
2103 ***********************************************************************/
2104 //------------------------------------------------------------------------------
2112 //------------------------------------------------------------------------------
2113 double gm2DDistanceToLineWithTol(const Pt3d* a_pt1,
2114  const Pt3d* a_pt2,
2115  double a_x,
2116  double a_y,
2117  double a_tol)
2118 {
2119  double a1, b1, c, mag;
2120  double dist;
2121  // see if the (x,y) is on infinite line
2122  a1 = a_pt2->x - a_pt1->x;
2123  b1 = a_pt2->y - a_pt1->y;
2124  mag = sqrt(a1 * a1 + b1 * b1);
2125  // handle case of line segment with length < tol distance to either point (pt1)
2126  if (mag <= a_tol)
2127  {
2128  return sqrt(sqr(a_pt1->x - a_x) + sqr(a_pt1->y - a_y));
2129  }
2130 
2131  // compute line equation
2132  double a, b;
2133  a = -b1 / mag;
2134  b = a1 / mag;
2135  c = -a * a_pt1->x - b * a_pt1->y;
2136  // compute distance from the line to (x,y)
2137  dist = a * a_x + b * a_y + c;
2138 
2139  return dist;
2140 } // gm2DDistanceToLineWithTol
2141 
2142 } // namespace xms
2143 
2144 //----- Tests ------------------------------------------------------------------
2145 
2146 #ifdef CXX_TEST
2147 
2148 //#include <RunTest.h>
2149 
2150 #include <boost/timer/timer.hpp>
2151 
2152 #include <xmscore/testing/TestTools.h>
2153 #include <xmsgeom/geometry/GmPolygon.t.h>
2154 #include <xmsgeom/geometry/geoms.t.h>
2155 #include <xmscore/misc/XmLog.h>
2156 
2157 //----- Namespace declaration --------------------------------------------------
2158 
2159 // namespace xms {
2160 
2167 //------------------------------------------------------------------------------
2169 //------------------------------------------------------------------------------
2170 GmPointInPolyUnitTests::GmPointInPolyUnitTests()
2171 : m_outPoly()
2172 , m_inPoly()
2173 , m_timer()
2174 , m_count(0)
2175 , m_status(XM_NODATA)
2176 {
2177 } // GmPointInPolyUnitTests::GmPointInPolyUnitTests
2178 //------------------------------------------------------------------------------
2180 //------------------------------------------------------------------------------
2181 void GmPointInPolyUnitTests::DoTest()
2182 {
2183  Setup();
2184  CheckPoints();
2185  GetResults();
2186 } // GmPointInPolyUnitTests::DoTest
2187 //------------------------------------------------------------------------------
2189 //------------------------------------------------------------------------------
2190 void GmPointInPolyUnitTests::Setup()
2191 {
2192  GmPolygonUnitTests::SetUpPolyWithHole(m_outPoly, m_inPoly);
2193 } // GmPointInPolyUnitTests::Setup
2194 //------------------------------------------------------------------------------
2196 //------------------------------------------------------------------------------
2197 void GmPointInPolyUnitTests::CheckPoints()
2198 {
2199  xms::Pt3d pt;
2200  const double xStart = -5;
2201  const double xEnd = 55;
2202  const double yStart = 35;
2203  const double yEnd = -5;
2204  const double kIncrement = 1e-1;
2205  m_timer.start();
2206  for (double y = yStart; y >= yEnd; y -= kIncrement)
2207  {
2208  for (double x = xStart; x <= xEnd; x += kIncrement)
2209  {
2210  pt.Set(x, y, 0);
2211  CheckPoint(pt);
2212  ++m_count;
2213  }
2214  }
2215 } // GmPointInPolyUnitTests::CheckPoints
2216 //--------------------------------------------------------------------------
2218 //--------------------------------------------------------------------------
2219 void GmPointInPolyUnitTests::GetResults()
2220 {
2221  // Using elapsed user + system time, not wall clock time, so this should be
2222  // similar on other computers. "For a production process, the wall clock
2223  // time may be what is most important. To study the efficiency of code,
2224  // total CPU time (user + system) is often a much better measure."
2225  // http://www.boost.org/doc/libs/1_59_0/libs/timer/doc/cpu_timers.html
2226 
2227  // Get elapsed user + system time, not wall clock time
2228  boost::timer::cpu_times const elapsed_times(m_timer.elapsed());
2229  boost::timer::nanosecond_type time = elapsed_times.system + elapsed_times.user;
2230 
2231  // convert nanoseconds to seconds
2232  const double NANO_PER_SEC = 1e9;
2233  double seconds = time / NANO_PER_SEC;
2234 
2235  // Make sure we checked a bunch of points
2236  int const kCountTotal = 240000;
2237  TS_ASSERT_EQUALS(m_count, kCountTotal);
2238 
2239  // For speed tests we don't check correctness of every point. We just
2240  // make sure m_status has changed.
2241  TS_ASSERT(m_status != XM_NODATA);
2242 
2243  // Check that the time is below a max time and write the time to the log
2244  // so we can check this in release.
2245  TS_ASSERT(seconds < MaxTime());
2246  XM_LOG(xmlog::debug, std::to_string((long long)seconds));
2247 } // GmPointInPolyUnitTests::GetResults
2248 
2251 class GmPointInPolyTester_gmPointInPolygon2D : public GmPointInPolyUnitTests
2252 {
2253 public:
2254  //----------------------------------------------------------------------------
2256  //----------------------------------------------------------------------------
2257  GmPointInPolyTester_gmPointInPolygon2D()
2258  : GmPointInPolyUnitTests()
2259  {
2260  }
2261 
2262 private:
2263  //----------------------------------------------------------------------------
2266  //----------------------------------------------------------------------------
2267  virtual void CheckPoint(const xms::Pt3d& a_pt) override
2268  {
2271  int outer = xms::gmPointInPolygon2D(&m_outPoly[0], m_outPoly.size(), a_pt);
2272  int inner = xms::gmPointInPolygon2D(&m_inPoly[0], m_inPoly.size(), a_pt);
2273  if (outer == 1 && inner == -1)
2274  {
2275  m_status = 1; // in
2276  }
2277  else if (outer == -1 || inner == 1)
2278  {
2279  m_status = 0; // completely out or in the hole
2280  }
2281  else if (outer == 0 || inner == 0)
2282  {
2283  m_status = 2; // on the outer poly or the inner poly
2284  }
2285  } // CheckPoint
2286  //----------------------------------------------------------------------------
2289  //----------------------------------------------------------------------------
2290  virtual double MaxTime() override
2291  {
2292  // It takes .09 - .17 seconds on my machine in release.
2293  return 0.5;
2294  // Hopefully big enough for the tests to pass on all machines but give us
2295  // an error if it starts going a lot slower for some reason.
2296  } // MaxTime
2297 
2298 private:
2300  XM_DISALLOW_COPY_AND_ASSIGN(GmPointInPolyTester_gmPointInPolygon2D);
2301 }; // GmPointInPolyTester_gmPointInPolygon2D
2302 
2307 //------------------------------------------------------------------------------
2310 //------------------------------------------------------------------------------
2311 // virtual
2312 // const CxxTest::TestGroup& GeomsXmsngUnitTests::group ()
2313 //{
2314 // return *CxxTest::TestGroup::GetGroup(CxxTest::TG_INTERMEDIATE);
2315 //} // GeomsXmsngUnitTests::group
2316 //------------------------------------------------------------------------------
2319 // 10 *------*
2320 // | |
2321 // | *
2322 // | |
2323 // 0 *------*
2324 // 0 10
2326 //------------------------------------------------------------------------------
2327 void GeomsXmsngUnitTests::test_gmComputeCentroid()
2328 {
2329  using xms::Pt3d;
2330 
2331  // 10 *------*
2332  // | |
2333  // | *
2334  // | |
2335  // 0 *------*
2336  // 0 10
2337 
2338  std::vector<Pt3d> points{{0, 0, 0}, {10, 0, 0}, {10, 5, 0}, {10, 10, 0}, {0, 10, 0}};
2339  const double kDelta = 1e-5;
2340  TS_ASSERT_DELTA_PT3D(xms::gmComputeCentroid(points), Pt3d(6, 5, 0), kDelta);
2341  // Notice the 6 above means it's not the centroid, but the average
2342 } // GeomsXmsngUnitTests::test_gmComputeCentroid
2343 //------------------------------------------------------------------------------
2346 // 10 *------*
2347 // | |
2348 // | *
2349 // | |
2350 // 0 *------*
2351 // 0 10
2353 //------------------------------------------------------------------------------
2354 void GeomsXmsngUnitTests::test_gmComputePolygonCentroid()
2355 {
2356  using xms::Pt3d;
2357 
2358  std::vector<Pt3d> points{{0, 0, 0}, {10, 0, 0}, {10, 5, 0}, {10, 10, 0}, {0, 10, 0}};
2359  const double kDelta = 1e-5;
2360  TS_ASSERT_DELTA_PT3D(xms::gmComputePolygonCentroid(points), Pt3d(5, 5, 0), kDelta);
2361 } // GeomsXmsngUnitTests::test_gmComputePolygonCentroid
2362 //------------------------------------------------------------------------------
2366 // 30 - 45--44 41--40 37--36 33--32 29------28
2367 // | | | | | | | | | | |
2368 // 25 - 46 43--42 39--38 35--34 31--30 26--27
2369 // | | |
2370 // 20 - 47---48 2---3---4---5---6---7---8 25--24
2371 // | | | | |
2372 // 15 - 50---49 1 18---17 14---13 10---9 22--23
2373 // | | | | | | | | |
2374 // 10 - 51 0---19 16---15 12---11 21--20
2375 // | | |
2376 // 5 - 52 2---3 6---7 10--11 14--15 18--19
2377 // | | | | | | | | | | |
2378 // 0 - 0---1 4---5 8---9 12--13 16--17
2379 //
2380 // |---|---|---|---|---|---|---|---|---|---|
2381 // 0 5 10 15 20 25 30 35 40 45 50 55
2383 //------------------------------------------------------------------------------
2384 void GeomsXmsngUnitTests::test_gmPointInPolygon2D_Speed()
2385 {
2386 #ifndef _DEBUG
2387  GmPointInPolyTester_gmPointInPolygon2D tester;
2388  tester.DoTest();
2389 #endif
2390 } // GeomsXmsngUnitTests::test_gmPointInPolygon2D_Speed
2391 
2396 //------------------------------------------------------------------------------
2399 //------------------------------------------------------------------------------
2400 #ifndef CXXTEST4
2401 const CxxTest::TestGroup& GeomsXmsngIntermediateTests::group()
2402 {
2403  return *CxxTest::TestGroup::GetGroup(CxxTest::TG_INTERMEDIATE);
2404  // return CxxTest::TestSuite::group();
2405 }
2406 #endif
2407 //------------------------------------------------------------------------------
2410 // 15 - 0 1 2 3 4
2411 // |
2412 // 10 - 5 6---7---8 9
2413 // | | |
2414 // 5 - 10 11 12 13 14
2415 // | | |
2416 // 0 - 15 16--17--18 19
2417 // |
2418 // -15- 20 21 22 23 24
2419 //
2420 // |---|---|---|---|
2421 // -5 0 5 10 15
2423 //------------------------------------------------------------------------------
2424 void GeomsXmsngIntermediateTests::test_gmPointInPolygon2D()
2425 {
2426  using xms::Pt3d;
2427  using xms::gmPointInPolygon2D;
2428 
2429  // CCW, first point not repeated. We also test it CW below.
2430  xms::VecPt3d poly{{0, 0, 0}, {5, 0, 0}, {10, 0, 0}, {10, 5, 0}, {10, 10, 0}, {0, 10, 0}};
2431  for (size_t i = 0; i < 2; ++i)
2432  {
2433  if (i == 1)
2434  {
2435  // Make it CW and recheck
2436  std::reverse(poly.begin(), poly.end());
2437  }
2438  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(-5, 15, 0)), -1); // 0
2439  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(0, 15, 0)), -1); // 1
2440  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(5, 15, 0)), -1); // 2
2441  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(10, 15, 0)), -1); // 3
2442  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(15, 15, 0)), -1); // 4
2443  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(-5, 10, 0)), -1); // 5
2444  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(0, 10, 0)), 0); // 6
2445  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(5, 10, 0)), 0); // 7
2446  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(10, 10, 0)), 0); // 8
2447  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(15, 10, 0)), -1); // 9
2448  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(-5, 5, 0)), -1); // 10
2449  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(0, 5, 0)), 0); // 11
2450  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(5, 5, 0)), 1); // 12
2451  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(10, 5, 0)), 0); // 13
2452  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(15, 5, 0)), -1); // 14
2453  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(-5, 0, 0)), -1); // 15
2454  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(0, 0, 0)), 0); // 16
2455  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(5, 0, 0)), 0); // 17
2456  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(10, 0, 0)), 0); // 18
2457  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(15, 0, 0)), -1); // 19
2458  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(-5, -5, 0)), -1); // 20
2459  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(0, -5, 0)), -1); // 21
2460  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(5, -5, 0)), -1); // 22
2461  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(10, -5, 0)), -1); // 23
2462  TS_ASSERT_EQUALS(gmPointInPolygon2D(&poly[0], poly.size(), Pt3d(15, -5, 0)), -1); // 24
2463  }
2464 } // GeomsXmsngIntermediateTests::test_gmPointInPolygon2D
2465 
2466  //} // namespace xms
2467 
2468 #endif
GmPointInPolyUnitTests::Setup
virtual void Setup()
Setup the polygons.
Definition: geoms.cpp:2190
GmPointInPolyTester_gmPointInPolygon2D
For testing point in polygon speed.
Definition: geoms.cpp:2251
TS_ASSERT_DELTA_PT3D
#define TS_ASSERT_DELTA_PT3D(a, b, delta)
XM_LOG
#define XM_LOG(A, B)
GeomsXmsngIntermediateTests::group
virtual const CxxTest::TestGroup & group()
Defines the test group.
Definition: geoms.cpp:2401
Pt3d
Pt3< double > Pt3d
GmPointInPolyUnitTests::GetResults
virtual void GetResults()
Get the timer results and do some assertions.
Definition: geoms.cpp:2219
xms::GTEQ_TOL
bool GTEQ_TOL(_T A, _U B, _V tol)
xms::XM_DBL_LOWEST
static const double XM_DBL_LOWEST
xms::EQ_TOL
bool EQ_TOL(const _T &A, const _U &B, const _V &tolerance)
triangles.h
Functions dealing with triangles.
GmPointInPolyUnitTests::MaxTime
virtual double MaxTime()=0
GmPointInPolyTester_gmPointInPolygon2D::MaxTime
virtual double MaxTime() override
Maximum time test should take.
Definition: geoms.cpp:2290
GmPointInPolyUnitTests::CheckPoint
virtual void CheckPoint(const xms::Pt3d &a_pt)=0
GmPointInPolyTester_gmPointInPolygon2D::CheckPoint
virtual void CheckPoint(const xms::Pt3d &a_pt) override
Determine if a point is in or out.
Definition: geoms.cpp:2267
xms::LT_TOL
bool LT_TOL(_T A, _U B, _V tolerance)
GmPolygonUnitTests::SetUpPolyWithHole
static void SetUpPolyWithHole(xms::VecPt3d &a_outPoly, xms::VecPt3d &a_inPolys)
Used in tests to create a polygon with lots of segments and 2 holes.
Definition: GmPolygon.cpp:625
GmPointInPolyUnitTests::m_timer
boost::timer::cpu_timer m_timer
Timer.
Definition: geoms.t.h:65
XM_ZERO_TOL
#define XM_ZERO_TOL
VecInt
std::vector< int > VecInt
GmPointInPolyUnitTests::m_status
int m_status
Status (in, out, on) of at least one pt.
Definition: geoms.t.h:67
GmPointInPolyUnitTests::DoTest
void DoTest()
Run the test. This is the main function to call.
Definition: geoms.cpp:2181
geoms.h
Functions dealing with geometry.
XmConst.h
VecDbl
std::vector< double > VecDbl
xms::Mmax
_T Mmax(_T a, _U b)
XM_ENSURE_TRUE
#define XM_ENSURE_TRUE(...)
XM_ENSURE_TRUE_VOID_NO_ASSERT
#define XM_ENSURE_TRUE_VOID_NO_ASSERT(...)
VecPt3d
std::vector< Pt3d > VecPt3d
GmPointInPolyUnitTests::GmPointInPolyUnitTests
GmPointInPolyUnitTests()
constructor
Definition: geoms.cpp:2170
XM_NODATA
#define XM_NODATA
xms::GT_TOL
bool GT_TOL(_T A, _U B, _V tolerance)
GeomsXmsngUnitTests::test_gmPointInPolygon2D_Speed
void test_gmPointInPolygon2D_Speed()
Test lots of points for timing purposes. Only in release, not debug.
Definition: geoms.cpp:2384
XM_DISALLOW_COPY_AND_ASSIGN
#define XM_DISALLOW_COPY_AND_ASSIGN(TypeName)
xms::Pt3::Set
void Set(T a_x, T a_y, T a_z)
Pt2i
Pt2< int > Pt2i
GeomsXmsngUnitTests::test_gmComputePolygonCentroid
void test_gmComputePolygonCentroid()
Definition: geoms.cpp:2354
xms::sqr
_T sqr(const _T x)
GmPointInPolyUnitTests::m_count
int m_count
Number of points checked.
Definition: geoms.t.h:66
GmPointInPolyUnitTests::m_inPoly
std::vector< xms::Pt3d > m_inPoly
Input polygon.
Definition: geoms.t.h:64
GmPointInPolyUnitTests::CheckPoints
virtual void CheckPoints()
Check a lot of points to see if they are in the polys.
Definition: geoms.cpp:2197
GmPointInPolyUnitTests
Used for speed tests of various point in poly functions / classes.
Definition: geoms.t.h:42
XM_PI
#define XM_PI
xms::Mfabs
_T Mfabs(_T a)
GmPolygon.t.h
XmError.h
XmLog.h
xms::Mmin
_T Mmin(_T a, _U b)
GeomsXmsngUnitTests::test_gmComputeCentroid
void test_gmComputeCentroid()
Definition: geoms.cpp:2327
GmPointInPolyUnitTests::m_outPoly
std::vector< xms::Pt3d > m_outPoly
Output polygon.
Definition: geoms.t.h:63
TestTools.h
GmPointInPolyTester_gmPointInPolygon2D::GmPointInPolyTester_gmPointInPolygon2D
GmPointInPolyTester_gmPointInPolygon2D()
Constructor.
Definition: geoms.cpp:2257
GmPolygon.h
xms::XM_DBL_HIGHEST
static const double XM_DBL_HIGHEST