Thea
AxisAlignedBox3.hpp
1 //============================================================================
2 //
3 // This file is part of the Thea toolkit.
4 //
5 // This software is distributed under the BSD license, as detailed in the
6 // accompanying LICENSE.txt file. Portions are derived from other works:
7 // their respective licenses and copyright information are reproduced in
8 // LICENSE.txt and/or in the relevant source files.
9 //
10 // Author: Siddhartha Chaudhuri
11 // First version: 2009
12 //
13 //============== License for line-AAB distance calculation ===================
14 //
15 // David Eberly, Geometric Tools, Redmond WA 98052
16 // Copyright (c) 1998-2017
17 // Distributed under the Boost Software License, Version 1.0.
18 // http://www.boost.org/LICENSE_1_0.txt
19 // http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
20 // File Version: 3.0.0 (2016/06/19)
21 //
22 //============================================================================
23 
24 #ifndef __Thea_AxisAlignedBox3_hpp__
25 #define __Thea_AxisAlignedBox3_hpp__
26 
27 #include "Common.hpp"
28 #include "AxisAlignedBoxN.hpp"
29 #include <array>
30 #include <utility>
31 
32 namespace Thea {
33 
34 namespace Internal {
35 
36 // From https://www.geometrictools.com/GTEngine/Include/Mathematics/GteDistLine3AlignedBox3.h
37 template <typename T>
38 class LineBoxDistance3
39 {
40  public:
41  typedef Vector<3, T> VectorT;
42 
43  struct Result
44  {
45  T sqrDistance;
46  T lineParameter;
47  };
48 
49  // Compute the distance and closest point between a line and an
50  // axis-aligned box whose center is the origin. On input, 'point' is the
51  // line origin and 'direction' is the line direction. On output, 'point'
52  // is the point on the box closest to the line. The 'direction' is
53  // non-const to allow transforming the problem into the first octant.
54  static void DoQuery(VectorT & point, VectorT & direction,
55  VectorT const & boxExtent, Result & result);
56 
57  private:
58  static void Face(int i0, int i1, int i2, VectorT & pnt,
59  VectorT const & dir, VectorT const & PmE,
60  VectorT const & boxExtent, Result & result);
61 
62  static void CaseNoZeros(VectorT & pnt, VectorT const & dir,
63  VectorT const & boxExtent, Result & result);
64 
65  static void Case0(int i0, int i1, int i2, VectorT & pnt,
66  VectorT const & dir, VectorT const & boxExtent,
67  Result & result);
68 
69  static void Case00(int i0, int i1, int i2, VectorT & pnt,
70  VectorT const & dir, VectorT const & boxExtent,
71  Result & result);
72 
73  static void Case000(VectorT & pnt, VectorT const & boxExtent,
74  Result & result);
75 
76 }; // class LineBoxDistance3
77 
78 } // namespace Internal
79 
81 template <typename T>
82 class /* THEA_API */ AxisAlignedBoxN<3, T> : public Internal::AxisAlignedBoxNBase<3, T>
83 {
84  private:
86 
87  public:
88  typedef typename BaseT::VectorT VectorT;
89 
91  AxisAlignedBoxN() : BaseT() {}
92 
94  AxisAlignedBoxN(VectorT const & v) : BaseT(v) {}
95 
97  AxisAlignedBoxN(VectorT const & lo_, VectorT const & hi_) : BaseT(lo_, hi_) {}
98 
104  template <typename TransformT> AxisAlignedBoxN transformAndBound(TransformT const & tr) const
105  {
106  VectorT const & lo_ = this->getLow();
107  VectorT const & hi_ = this->getHigh();
108 
109  // We sequentially compute the corners in the following order:
110  // 0, 6, 5, 1, 2, 4, 7, 3
111  // This sequence allows us to only change one member at a time to get at all corners. For each one, we transform it and
112  // merge the resulting point.
113 
114  // min min min
115  VectorT current_corner = lo_;
116  AxisAlignedBoxN result = AxisAlignedBoxN(tr * current_corner);
117 
118  // min min max
119  current_corner.z() = hi_.z();
120  result.merge(tr * current_corner);
121 
122  // min max max
123  current_corner.y() = hi_.y();
124  result.merge(tr * current_corner);
125 
126  // min max min
127  current_corner.z() = lo_.z();
128  result.merge(tr * current_corner);
129 
130  // max max min
131  current_corner.x() = hi_.x();
132  result.merge(tr * current_corner);
133 
134  // max max max
135  current_corner.z() = hi_.z();
136  result.merge(tr * current_corner);
137 
138  // max min max
139  current_corner.y() = lo_.y();
140  result.merge(tr * current_corner);
141 
142  // max min min
143  current_corner.z() = lo_.z();
144  result.merge(tr * current_corner);
145 
146  return result;
147  }
148 
159  std::pair<uintx, uintx> const & getEdge(intx i) const
160  {
161  static std::pair<uintx, uintx> const EDGES[12] = {
162  { 0b000, 0b100 },
163  { 0b000, 0b010 },
164  { 0b000, 0b001 },
165 
166  { 0b011, 0b111 },
167  { 0b011, 0b001 },
168  { 0b011, 0b010 },
169 
170  { 0b101, 0b100 },
171  { 0b101, 0b111 },
172  { 0b101, 0b001 },
173 
174  { 0b110, 0b111 },
175  { 0b110, 0b100 },
176  { 0b110, 0b010 }
177  };
178 
179  debugAssertM(i >= 0 && i < 12, "AxisAlignedBox3: Edge index out of bounds");
180 
181  return EDGES[i];
182  }
183 
194  std::array<uintx, 4> const & getFace(intx i) const
195  {
196  static std::array<uintx, 4> const FACES[6] = {
197  { 0, 4, 6, 2 }, // -X
198  { 1, 3, 7, 5 }, // +X
199 
200  { 0, 1, 5, 4 }, // -Y
201  { 2, 6, 7, 3 }, // +Y
202 
203  { 0, 2, 3, 1 }, // -Z
204  { 4, 5, 7, 6 } // +Z
205  };
206 
207  debugAssertM(i >= 0 && i < 6, "AxisAlignedBox3: Face index out of bounds");
208 
209  return FACES[i];
210  }
211 
219  VectorT const & getFaceNormal(intx i) const
220  {
221  // Synced with the face order in getFace()
222  static VectorT const NORMALS[6] = {
223  -VectorT::UnitX(),
224  VectorT::UnitX(),
225 
226  -VectorT::UnitY(),
227  VectorT::UnitY(),
228 
229  -VectorT::UnitZ(),
230  VectorT::UnitZ()
231  };
232 
233  debugAssertM(i >= 0 && i < 6, "AxisAlignedBox3: Face index out of bounds");
234 
235  return NORMALS[i];
236  }
237 
238  using BaseT::distance;
239  using BaseT::squaredDistance;
240 
241  T distance(LineN<3, T> const & line) const { return std::sqrt(squaredDistance(line)); }
242 
243  T squaredDistance(LineN<3, T> const & line, VectorT * this_pt = nullptr, VectorT * line_pt = nullptr) const
244  {
245  // Translate the line and box so that the box has center at the origin.
246  VectorT boxCenter = this->getCenter(), boxExtent = 0.5f * this->getExtent();
247  VectorT point = line.getPoint() - boxCenter;
248  VectorT direction = line.getDirection();
249 
250  typename Internal::LineBoxDistance3<T>::Result result;
251  Internal::LineBoxDistance3<T>::DoQuery(point, direction, boxExtent, result);
252 
253  if (this_pt) *this_pt = boxCenter + point;
254  if (line_pt) *line_pt = line.getPoint() + result.lineParameter * line.getDirection();
255 
256  return result.sqrDistance;
257  }
258 
259  T distance(LineSegmentN<3, T> const & seg) const { return std::sqrt(squaredDistance(seg)); }
260 
261  T squaredDistance(LineSegmentN<3, T> const & seg, VectorT * this_pt = nullptr, VectorT * seg_pt = nullptr) const
262  {
263  // Translate the line and box so that the box has center at the origin.
264  VectorT boxCenter = this->getCenter(), boxExtent = 0.5f * this->getExtent();
265  VectorT e0 = seg.getEndpoint(0);
266  VectorT point = e0 - boxCenter;
267 
268  // FIXME: Does DoQuery actually need a unit length direction vector?
269  T seg_len = seg.length();
270  VectorT seg_dir = (seg_len >= Math::eps<T>() ? seg.getDirection() / seg_len : VectorT::UnitX());
271  VectorT direction = seg_dir; // can be overwritten by DoQuery
272 
273  typename Internal::LineBoxDistance3<T>::Result result;
274  Internal::LineBoxDistance3<T>::DoQuery(point, direction, boxExtent, result);
275 
276  if (result.lineParameter >= 0 && result.lineParameter <= seg_len)
277  {
278  if (this_pt) *this_pt = boxCenter + point;
279  if (seg_pt) *seg_pt = e0 + result.lineParameter * seg_dir;
280  }
281  else
282  {
283  VectorT e1 = seg.getEndpoint(1);
284  VectorT c0 = this->closestPoint(e0);
285  VectorT c1 = this->closestPoint(e1);
286  Real sqdist0 = (c0 - e0).squaredNorm();
287  Real sqdist1 = (c1 - e1).squaredNorm();
288  if (sqdist0 < sqdist1)
289  {
290  if (this_pt) *this_pt = c0;
291  if (seg_pt) *seg_pt = e0;
292  result.sqrDistance = sqdist0;
293  }
294  else
295  {
296  if (this_pt) *this_pt = c1;
297  if (seg_pt) *seg_pt = e1;
298  result.sqrDistance = sqdist1;
299  }
300  }
301 
302  return result.sqrDistance;
303  }
304 
305  T distance(RayN<3, T> const & ray) const { return std::sqrt(squaredDistance(ray)); }
306 
307  T squaredDistance(RayN<3, T> const & ray, VectorT * this_pt = nullptr, VectorT * ray_pt = nullptr) const
308  {
309  // Translate the line and box so that the box has center at the origin.
310  VectorT boxCenter = this->getCenter(), boxExtent = 0.5f * this->getExtent();
311  VectorT point = ray.getOrigin() - boxCenter;
312  VectorT ray_dir = ray.getDirection().normalized(); // FIXME: Does DoQuery actually need a unit length direction vector?
313  VectorT direction = ray_dir; // can be overwritten by DoQuery
314 
315  typename Internal::LineBoxDistance3<T>::Result result;
316  Internal::LineBoxDistance3<T>::DoQuery(point, direction, boxExtent, result);
317 
318  if (result.lineParameter >= 0)
319  {
320  if (this_pt) *this_pt = boxCenter + point;
321  if (ray_pt) *ray_pt = ray.getOrigin() + result.lineParameter * ray_dir;
322  }
323  else
324  {
325  VectorT c = this->closestPoint(ray.getOrigin());
326  result.sqrDistance = (c - ray.getOrigin()).squaredNorm();
327 
328  if (this_pt) *this_pt = c;
329  if (ray_pt) *ray_pt = ray.getOrigin();
330  }
331 
332  return result.sqrDistance;
333  }
334 
335 }; // class AxisAlignedBoxN<3, T>
336 
337 #ifdef THEA_EXPORT_INSTANTIATION
338  template class THEA_API AxisAlignedBoxN<3, Real>;
339 #endif
340 
342 typedef AxisAlignedBoxN<3, Real> AxisAlignedBox3;
343 
344 namespace Internal {
345 
346 //=============================================================================================================================
347 //
348 // Line-to-axis-aligned-box distance calculations in 3D. Note that "extent" below refers to half the extent of the box above.
349 //
350 //=============================================================================================================================
351 
352 template <typename T>
353 void
354 LineBoxDistance3<T>::DoQuery(VectorT & point, VectorT & direction,
355  VectorT const & boxExtent, Result & result)
356 {
357  result.sqrDistance = (T)0;
358  result.lineParameter = (T)0;
359 
360  // Apply reflections so that direction vector has nonnegative components.
361  bool reflect[3];
362  for (int i = 0; i < 3; ++i)
363  {
364  if (direction[i] < (T)0)
365  {
366  point[i] = -point[i];
367  direction[i] = -direction[i];
368  reflect[i] = true;
369  }
370  else
371  {
372  reflect[i] = false;
373  }
374  }
375 
376  if (direction[0] > (T)0)
377  {
378  if (direction[1] > (T)0)
379  {
380  if (direction[2] > (T)0) // (+,+,+)
381  {
382  CaseNoZeros(point, direction, boxExtent, result);
383  }
384  else // (+,+,0)
385  {
386  Case0(0, 1, 2, point, direction, boxExtent, result);
387  }
388  }
389  else
390  {
391  if (direction[2] > (T)0) // (+,0,+)
392  {
393  Case0(0, 2, 1, point, direction, boxExtent, result);
394  }
395  else // (+,0,0)
396  {
397  Case00(0, 1, 2, point, direction, boxExtent, result);
398  }
399  }
400  }
401  else
402  {
403  if (direction[1] > (T)0)
404  {
405  if (direction[2] > (T)0) // (0,+,+)
406  {
407  Case0(1, 2, 0, point, direction, boxExtent, result);
408  }
409  else // (0,+,0)
410  {
411  Case00(1, 0, 2, point, direction, boxExtent, result);
412  }
413  }
414  else
415  {
416  if (direction[2] > (T)0) // (0,0,+)
417  {
418  Case00(2, 0, 1, point, direction, boxExtent, result);
419  }
420  else // (0,0,0)
421  {
422  Case000(point, boxExtent, result);
423  }
424  }
425  }
426 
427  // Undo the reflections applied previously.
428  for (int i = 0; i < 3; ++i)
429  {
430  if (reflect[i])
431  {
432  point[i] = -point[i];
433  }
434  }
435 }
436 
437 template <typename T>
438 void
439 LineBoxDistance3<T>::Face(int i0, int i1,
440  int i2, VectorT & pnt, VectorT const & dir,
441  VectorT const & PmE, VectorT const & boxExtent, Result & result)
442 {
443  VectorT PpE;
444  T lenSqr, inv, tmp, param, t, delta;
445 
446  PpE[i1] = pnt[i1] + boxExtent[i1];
447  PpE[i2] = pnt[i2] + boxExtent[i2];
448  if (dir[i0] * PpE[i1] >= dir[i1] * PmE[i0])
449  {
450  if (dir[i0] * PpE[i2] >= dir[i2] * PmE[i0])
451  {
452  // v[i1] >= -e[i1], v[i2] >= -e[i2] (distance = 0)
453  pnt[i0] = boxExtent[i0];
454  inv = ((T)1) / dir[i0];
455  pnt[i1] -= dir[i1] * PmE[i0] * inv;
456  pnt[i2] -= dir[i2] * PmE[i0] * inv;
457  result.lineParameter = -PmE[i0] * inv;
458  }
459  else
460  {
461  // v[i1] >= -e[i1], v[i2] < -e[i2]
462  lenSqr = dir[i0] * dir[i0] + dir[i2] * dir[i2];
463  tmp = lenSqr * PpE[i1] - dir[i1] * (dir[i0] * PmE[i0] +
464  dir[i2] * PpE[i2]);
465  if (tmp <= ((T)2)*lenSqr * boxExtent[i1])
466  {
467  t = tmp / lenSqr;
468  lenSqr += dir[i1] * dir[i1];
469  tmp = PpE[i1] - t;
470  delta = dir[i0] * PmE[i0] + dir[i1] * tmp + dir[i2] * PpE[i2];
471  param = -delta / lenSqr;
472  result.sqrDistance += PmE[i0] * PmE[i0] + tmp * tmp +
473  PpE[i2] * PpE[i2] + delta * param;
474 
475  result.lineParameter = param;
476  pnt[i0] = boxExtent[i0];
477  pnt[i1] = t - boxExtent[i1];
478  pnt[i2] = -boxExtent[i2];
479  }
480  else
481  {
482  lenSqr += dir[i1] * dir[i1];
483  delta = dir[i0] * PmE[i0] + dir[i1] * PmE[i1] + dir[i2] * PpE[i2];
484  param = -delta / lenSqr;
485  result.sqrDistance += PmE[i0] * PmE[i0] + PmE[i1] * PmE[i1]
486  + PpE[i2] * PpE[i2] + delta * param;
487 
488  result.lineParameter = param;
489  pnt[i0] = boxExtent[i0];
490  pnt[i1] = boxExtent[i1];
491  pnt[i2] = -boxExtent[i2];
492  }
493  }
494  }
495  else
496  {
497  if (dir[i0] * PpE[i2] >= dir[i2] * PmE[i0])
498  {
499  // v[i1] < -e[i1], v[i2] >= -e[i2]
500  lenSqr = dir[i0] * dir[i0] + dir[i1] * dir[i1];
501  tmp = lenSqr * PpE[i2] - dir[i2] * (dir[i0] * PmE[i0] +
502  dir[i1] * PpE[i1]);
503  if (tmp <= ((T)2)*lenSqr * boxExtent[i2])
504  {
505  t = tmp / lenSqr;
506  lenSqr += dir[i2] * dir[i2];
507  tmp = PpE[i2] - t;
508  delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * tmp;
509  param = -delta / lenSqr;
510  result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] +
511  tmp * tmp + delta * param;
512 
513  result.lineParameter = param;
514  pnt[i0] = boxExtent[i0];
515  pnt[i1] = -boxExtent[i1];
516  pnt[i2] = t - boxExtent[i2];
517  }
518  else
519  {
520  lenSqr += dir[i2] * dir[i2];
521  delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * PmE[i2];
522  param = -delta / lenSqr;
523  result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] +
524  PmE[i2] * PmE[i2] + delta * param;
525 
526  result.lineParameter = param;
527  pnt[i0] = boxExtent[i0];
528  pnt[i1] = -boxExtent[i1];
529  pnt[i2] = boxExtent[i2];
530  }
531  }
532  else
533  {
534  // v[i1] < -e[i1], v[i2] < -e[i2]
535  lenSqr = dir[i0] * dir[i0] + dir[i2] * dir[i2];
536  tmp = lenSqr * PpE[i1] - dir[i1] * (dir[i0] * PmE[i0] +
537  dir[i2] * PpE[i2]);
538  if (tmp >= (T)0)
539  {
540  // v[i1]-edge is closest
541  if (tmp <= ((T)2)*lenSqr * boxExtent[i1])
542  {
543  t = tmp / lenSqr;
544  lenSqr += dir[i1] * dir[i1];
545  tmp = PpE[i1] - t;
546  delta = dir[i0] * PmE[i0] + dir[i1] * tmp + dir[i2] * PpE[i2];
547  param = -delta / lenSqr;
548  result.sqrDistance += PmE[i0] * PmE[i0] + tmp * tmp +
549  PpE[i2] * PpE[i2] + delta * param;
550 
551  result.lineParameter = param;
552  pnt[i0] = boxExtent[i0];
553  pnt[i1] = t - boxExtent[i1];
554  pnt[i2] = -boxExtent[i2];
555  }
556  else
557  {
558  lenSqr += dir[i1] * dir[i1];
559  delta = dir[i0] * PmE[i0] + dir[i1] * PmE[i1]
560  + dir[i2] * PpE[i2];
561  param = -delta / lenSqr;
562  result.sqrDistance += PmE[i0] * PmE[i0] + PmE[i1] * PmE[i1]
563  + PpE[i2] * PpE[i2] + delta * param;
564 
565  result.lineParameter = param;
566  pnt[i0] = boxExtent[i0];
567  pnt[i1] = boxExtent[i1];
568  pnt[i2] = -boxExtent[i2];
569  }
570  return;
571  }
572 
573  lenSqr = dir[i0] * dir[i0] + dir[i1] * dir[i1];
574  tmp = lenSqr * PpE[i2] - dir[i2] * (dir[i0] * PmE[i0] +
575  dir[i1] * PpE[i1]);
576  if (tmp >= (T)0)
577  {
578  // v[i2]-edge is closest
579  if (tmp <= ((T)2)*lenSqr * boxExtent[i2])
580  {
581  t = tmp / lenSqr;
582  lenSqr += dir[i2] * dir[i2];
583  tmp = PpE[i2] - t;
584  delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * tmp;
585  param = -delta / lenSqr;
586  result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] +
587  tmp * tmp + delta * param;
588 
589  result.lineParameter = param;
590  pnt[i0] = boxExtent[i0];
591  pnt[i1] = -boxExtent[i1];
592  pnt[i2] = t - boxExtent[i2];
593  }
594  else
595  {
596  lenSqr += dir[i2] * dir[i2];
597  delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1]
598  + dir[i2] * PmE[i2];
599  param = -delta / lenSqr;
600  result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1]
601  + PmE[i2] * PmE[i2] + delta * param;
602 
603  result.lineParameter = param;
604  pnt[i0] = boxExtent[i0];
605  pnt[i1] = -boxExtent[i1];
606  pnt[i2] = boxExtent[i2];
607  }
608  return;
609  }
610 
611  // (v[i1],v[i2])-corner is closest
612  lenSqr += dir[i2] * dir[i2];
613  delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * PpE[i2];
614  param = -delta / lenSqr;
615  result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1]
616  + PpE[i2] * PpE[i2] + delta * param;
617 
618  result.lineParameter = param;
619  pnt[i0] = boxExtent[i0];
620  pnt[i1] = -boxExtent[i1];
621  pnt[i2] = -boxExtent[i2];
622  }
623  }
624 }
625 
626 template <typename T>
627 void
628 LineBoxDistance3<T>::CaseNoZeros(VectorT & pnt, VectorT const & dir,
629  VectorT const & boxExtent, Result & result)
630 {
631  VectorT PmE = pnt - boxExtent;
632  T prodDxPy = dir[0] * PmE[1];
633  T prodDyPx = dir[1] * PmE[0];
634  T prodDzPx, prodDxPz, prodDzPy, prodDyPz;
635 
636  if (prodDyPx >= prodDxPy)
637  {
638  prodDzPx = dir[2] * PmE[0];
639  prodDxPz = dir[0] * PmE[2];
640  if (prodDzPx >= prodDxPz)
641  {
642  // line intersects x = e0
643  Face(0, 1, 2, pnt, dir, PmE, boxExtent, result);
644  }
645  else
646  {
647  // line intersects z = e2
648  Face(2, 0, 1, pnt, dir, PmE, boxExtent, result);
649  }
650  }
651  else
652  {
653  prodDzPy = dir[2] * PmE[1];
654  prodDyPz = dir[1] * PmE[2];
655  if (prodDzPy >= prodDyPz)
656  {
657  // line intersects y = e1
658  Face(1, 2, 0, pnt, dir, PmE, boxExtent, result);
659  }
660  else
661  {
662  // line intersects z = e2
663  Face(2, 0, 1, pnt, dir, PmE, boxExtent, result);
664  }
665  }
666 }
667 
668 template <typename T>
669 void
670 LineBoxDistance3<T>::Case0(int i0, int i1,
671  int i2, VectorT & pnt, VectorT const & dir,
672  VectorT const & boxExtent, Result & result)
673 {
674  T PmE0 = pnt[i0] - boxExtent[i0];
675  T PmE1 = pnt[i1] - boxExtent[i1];
676  T prod0 = dir[i1] * PmE0;
677  T prod1 = dir[i0] * PmE1;
678  T delta, invLSqr, inv;
679 
680  if (prod0 >= prod1)
681  {
682  // line intersects P[i0] = e[i0]
683  pnt[i0] = boxExtent[i0];
684 
685  T PpE1 = pnt[i1] + boxExtent[i1];
686  delta = prod0 - dir[i0] * PpE1;
687  if (delta >= (T)0)
688  {
689  invLSqr = ((T)1) / (dir[i0] * dir[i0] + dir[i1] * dir[i1]);
690  result.sqrDistance += delta * delta * invLSqr;
691  pnt[i1] = -boxExtent[i1];
692  result.lineParameter = -(dir[i0] * PmE0 + dir[i1] * PpE1) * invLSqr;
693  }
694  else
695  {
696  inv = ((T)1) / dir[i0];
697  pnt[i1] -= prod0 * inv;
698  result.lineParameter = -PmE0 * inv;
699  }
700  }
701  else
702  {
703  // line intersects P[i1] = e[i1]
704  pnt[i1] = boxExtent[i1];
705 
706  T PpE0 = pnt[i0] + boxExtent[i0];
707  delta = prod1 - dir[i1] * PpE0;
708  if (delta >= (T)0)
709  {
710  invLSqr = ((T)1) / (dir[i0] * dir[i0] + dir[i1] * dir[i1]);
711  result.sqrDistance += delta * delta * invLSqr;
712  pnt[i0] = -boxExtent[i0];
713  result.lineParameter = -(dir[i0] * PpE0 + dir[i1] * PmE1) * invLSqr;
714  }
715  else
716  {
717  inv = ((T)1) / dir[i1];
718  pnt[i0] -= prod1 * inv;
719  result.lineParameter = -PmE1 * inv;
720  }
721  }
722 
723  if (pnt[i2] < -boxExtent[i2])
724  {
725  delta = pnt[i2] + boxExtent[i2];
726  result.sqrDistance += delta * delta;
727  pnt[i2] = -boxExtent[i2];
728  }
729  else if (pnt[i2] > boxExtent[i2])
730  {
731  delta = pnt[i2] - boxExtent[i2];
732  result.sqrDistance += delta * delta;
733  pnt[i2] = boxExtent[i2];
734  }
735 }
736 
737 template <typename T>
738 void
739 LineBoxDistance3<T>::Case00(int i0, int i1,
740  int i2, VectorT & pnt, VectorT const & dir,
741  VectorT const & boxExtent, Result & result)
742 {
743  T delta;
744 
745  result.lineParameter = (boxExtent[i0] - pnt[i0]) / dir[i0];
746 
747  pnt[i0] = boxExtent[i0];
748 
749  if (pnt[i1] < -boxExtent[i1])
750  {
751  delta = pnt[i1] + boxExtent[i1];
752  result.sqrDistance += delta * delta;
753  pnt[i1] = -boxExtent[i1];
754  }
755  else if (pnt[i1] > boxExtent[i1])
756  {
757  delta = pnt[i1] - boxExtent[i1];
758  result.sqrDistance += delta * delta;
759  pnt[i1] = boxExtent[i1];
760  }
761 
762  if (pnt[i2] < -boxExtent[i2])
763  {
764  delta = pnt[i2] + boxExtent[i2];
765  result.sqrDistance += delta * delta;
766  pnt[i2] = -boxExtent[i2];
767  }
768  else if (pnt[i2] > boxExtent[i2])
769  {
770  delta = pnt[i2] - boxExtent[i2];
771  result.sqrDistance += delta * delta;
772  pnt[i2] = boxExtent[i2];
773  }
774 }
775 
776 template <typename T>
777 void
778 LineBoxDistance3<T>::Case000(VectorT & pnt, VectorT const & boxExtent, Result & result)
779 {
780  T delta;
781 
782  if (pnt[0] < -boxExtent[0])
783  {
784  delta = pnt[0] + boxExtent[0];
785  result.sqrDistance += delta * delta;
786  pnt[0] = -boxExtent[0];
787  }
788  else if (pnt[0] > boxExtent[0])
789  {
790  delta = pnt[0] - boxExtent[0];
791  result.sqrDistance += delta * delta;
792  pnt[0] = boxExtent[0];
793  }
794 
795  if (pnt[1] < -boxExtent[1])
796  {
797  delta = pnt[1] + boxExtent[1];
798  result.sqrDistance += delta * delta;
799  pnt[1] = -boxExtent[1];
800  }
801  else if (pnt[1] > boxExtent[1])
802  {
803  delta = pnt[1] - boxExtent[1];
804  result.sqrDistance += delta * delta;
805  pnt[1] = boxExtent[1];
806  }
807 
808  if (pnt[2] < -boxExtent[2])
809  {
810  delta = pnt[2] + boxExtent[2];
811  result.sqrDistance += delta * delta;
812  pnt[2] = -boxExtent[2];
813  }
814  else if (pnt[2] > boxExtent[2])
815  {
816  delta = pnt[2] - boxExtent[2];
817  result.sqrDistance += delta * delta;
818  pnt[2] = boxExtent[2];
819  }
820 }
821 
822 } // namespace Internal
823 
824 } // namespace Thea
825 
826 #endif
VectorT const & getPoint() const
Get a point on the line.
Definition: LineN.hpp:73
std::ptrdiff_t intx
A signed integer suitable for indexing a structure held in memory.
Definition: Platform.hpp:161
VectorT const & getOrigin() const
Get the origin of the ray.
Definition: RayN.hpp:50
A straight line in N-dimensional space, where N is any positive (non-zero) integer and T is a field...
A ray in N-dimensional space, having an originating point and a direction vector (not necessarily uni...
Definition: RayN.hpp:27
Root namespace for the Thea library.
[Internal] Base class for axis-aligned boxes in N-dimensional space, where N is any positive (non-zer...
VectorT const & getFaceNormal(intx i) const
Get the outward unit normal vector for face number i of the box, where i is between 0 and 5...
VectorT const & getDirection() const
Get the unnormalized direction vector of the segment from the first endpoint to the second...
std::array< uintx, 4 > const & getFace(intx i) const
Get the vertex indices of face number i of the box, where i is between 0 and 5.
AxisAlignedBoxN(VectorT const &lo_, VectorT const &hi_)
Constructor.
AxisAlignedBoxN()
Default constructor, creates a null box.
An axis-aligned box in N-dimensional space, where N is any positive (non-zero) integer and T is a fie...
std::pair< uintx, uintx > const & getEdge(intx i) const
Get the endpoint indices of edge number i of the box, where i is between 0 and 11.
void merge(VectorT const &p)
Grow to include a point.
AxisAlignedBoxN< 3, Real > AxisAlignedBox3
The default axis-aligned box class in 3-dimensional real space.
A straight line in N-dimensional space, where N is any positive (non-zero) integer and T is a field...
Definition: LineN.hpp:25
AxisAlignedBoxN transformAndBound(TransformT const &tr) const
Transform the box and return a new axis-aligned box which tightly encloses the result.
VectorT const & getDirection() const
Get the unit direction vector of the line.
Definition: LineN.hpp:76
VectorT const & getDirection() const
Get the direction of the ray.
Definition: RayN.hpp:56
Vector< N, T > VectorT
N-dimensional vector.
VectorT getEndpoint(int i) const
Get an endpoint of the line segment: 0 returns the first endpoint and 1 returns the second...
T length() const
Get the length of the line segment.
AxisAlignedBoxN(VectorT const &v)
Constructor.
void debugAssertM(CondT const &test, MessageT const &msg)
Check if a test condition is true, and immediately abort the program with an error code if not...
Definition: Common.hpp:52