DistTriangle3AlignedBox3.h 5.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131
  1. // David Eberly, Geometric Tools, Redmond WA 98052
  2. // Copyright (c) 1998-2020
  3. // Distributed under the Boost Software License, Version 1.0.
  4. // https://www.boost.org/LICENSE_1_0.txt
  5. // https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
  6. // Version: 4.0.2019.08.13
  7. #pragma once
  8. #include <Mathematics/DCPQuery.h>
  9. #include <Mathematics/LCPSolver.h>
  10. #include <Mathematics/AlignedBox.h>
  11. #include <Mathematics/Triangle.h>
  12. #include <Mathematics/Vector3.h>
  13. // Compute the distance between a triangle and an aligned box in 3D. The
  14. // algorithm is based on using an LCP solver for the convex quadratic
  15. // programming problem. For details, see
  16. // https://www.geometrictools.com/Documentation/ConvexQuadraticProgramming.pdf
  17. namespace WwiseGTE
  18. {
  19. template <typename Real>
  20. class DCPQuery<Real, Triangle3<Real>, AlignedBox3<Real>>
  21. {
  22. public:
  23. struct Result
  24. {
  25. bool queryIsSuccessful;
  26. // These members are valid only when queryIsSuccessful is true;
  27. // otherwise, they are all set to zero.
  28. Real distance, sqrDistance;
  29. std::array<Real, 3> triangleParameter, boxParameter;
  30. Vector3<Real> closestPoint[2];
  31. // The number of iterations used by LCPSolver regardless of
  32. // whether the query is successful.
  33. int numLCPIterations;
  34. };
  35. // The default maximum iterations is 81 (n = 9, maxIterations = n*n).
  36. // If the solver fails to converge, try increasing the maximum number
  37. // of iterations.
  38. void SetMaxLCPIterations(int maxLCPIterations)
  39. {
  40. mLCP.SetMaxIterations(maxLCPIterations);
  41. }
  42. Result operator()(Triangle3<Real> const& triangle, AlignedBox3<Real> const& box)
  43. {
  44. Result result;
  45. // Translate the triangle and aligned box so that the aligned box
  46. // becomes a canonical box.
  47. Vector3<Real> K = box.max - box.min;
  48. Vector3<Real> V = triangle.v[0] - box.min;
  49. Vector3<Real> E0 = triangle.v[1] - triangle.v[0];
  50. Vector3<Real> E1 = triangle.v[2] - triangle.v[0];
  51. // Compute quantities to initialize q and M in the LCP.
  52. Real dotVE0 = Dot(V, E0);
  53. Real dotVE1 = Dot(V, E1);
  54. Real dotE0E0 = Dot(E0, E0);
  55. Real dotE0E1 = Dot(E0, E1);
  56. Real dotE1E1 = Dot(E1, E1);
  57. // The LCP has 5 variables and 4 (nontrivial) inequality
  58. // constraints.
  59. std::array<Real, 9> q =
  60. {
  61. -V[0], -V[1], -V[2], dotVE0, dotVE1, K[0], K[1], K[2], (Real)1
  62. };
  63. std::array<std::array<Real, 9>, 9> M;
  64. M[0] = { (Real)1, (Real)0, (Real)0, -E0[0], -E1[0], (Real)1, (Real)0, (Real)0, (Real)0 };
  65. M[1] = { (Real)0, (Real)1, (Real)0, -E0[1], -E1[1], (Real)0, (Real)1, (Real)0, (Real)0 };
  66. M[2] = { (Real)0, (Real)0, (Real)1, -E0[2], -E1[2], (Real)0, (Real)0, (Real)1, (Real)0 };
  67. M[3] = { -E0[0], -E0[1], -E0[2], dotE0E0, dotE0E1, (Real)0, (Real)0, (Real)0, (Real)1 };
  68. M[4] = { -E1[0], -E1[1], -E1[2], dotE0E1, dotE1E1, (Real)0, (Real)0, (Real)0, (Real)1 };
  69. M[5] = { (Real)-1, (Real)0, (Real)0, (Real)0, (Real)0, (Real)0, (Real)0, (Real)0, (Real)0 };
  70. M[6] = { (Real)0, (Real)-1, (Real)0, (Real)0, (Real)0, (Real)0, (Real)0, (Real)0, (Real)0 };
  71. M[7] = { (Real)0, (Real)0, (Real)-1, (Real)0, (Real)0, (Real)0, (Real)0, (Real)0, (Real)0 };
  72. M[8] = { (Real)0, (Real)0, (Real)0, (Real)-1, (Real)-1, (Real)0, (Real)0, (Real)0, (Real)0 };
  73. std::array<Real, 9> w, z;
  74. if (mLCP.Solve(q, M, w, z))
  75. {
  76. result.queryIsSuccessful = true;
  77. result.triangleParameter[0] = (Real)1 - z[3] - z[4];
  78. result.triangleParameter[1] = z[3];
  79. result.triangleParameter[2] = z[4];
  80. result.closestPoint[0] = triangle.v[0] + z[3] * E0 + z[4] * E1;
  81. for (int i = 0; i < 3; ++i)
  82. {
  83. result.boxParameter[i] = z[i] + box.min[i];
  84. result.closestPoint[1][i] = result.boxParameter[i];
  85. }
  86. Vector3<Real> diff = result.closestPoint[1] - result.closestPoint[0];
  87. result.sqrDistance = Dot(diff, diff);
  88. result.distance = std::sqrt(result.sqrDistance);
  89. }
  90. else
  91. {
  92. // If you reach this case, the maximum number of iterations
  93. // was not specified to be large enough or there is a problem
  94. // due to floating-point rounding errors. If you believe the
  95. // latter is true, file a bug report.
  96. result.queryIsSuccessful = false;
  97. for (int i = 0; i < 3; ++i)
  98. {
  99. result.triangleParameter[i] = (Real)0;
  100. result.boxParameter[i] = (Real)0;
  101. result.closestPoint[0][i] = (Real)0;
  102. result.closestPoint[1][i] = (Real)0;
  103. }
  104. result.distance = (Real)0;
  105. result.sqrDistance = (Real)0;
  106. }
  107. result.numLCPIterations = mLCP.GetNumIterations();
  108. return result;
  109. }
  110. private:
  111. LCPSolver<Real, 9> mLCP;
  112. };
  113. }