DistAlignedBox3OrientedBox3.h 5.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150
  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/OrientedBox.h>
  12. #include <Mathematics/Vector3.h>
  13. // Compute the distance between an aligned box and an oriented 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, AlignedBox3<Real>, OrientedBox3<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> box0Parameter, box1Parameter;
  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. // Default maximum iterations is 144 (n = 12, 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()(AlignedBox3<Real> const& box0, OrientedBox3<Real> const& box1)
  43. {
  44. Result result;
  45. // Translate the boxes so that the aligned box becomes a canonical
  46. // box. Modify the oriented box coefficients to be nonnegative.
  47. Vector3<Real> K = box0.max - box0.min;
  48. Vector3<Real> delta = box1.center - box0.min;
  49. for (int i = 0; i < 3; ++i)
  50. {
  51. delta -= box1.extent[i] * box1.axis[i];
  52. }
  53. Vector3<Real> rotDelta;
  54. for (int i = 0; i < 3; ++i)
  55. {
  56. rotDelta[i] = Dot(box1.axis[i], delta);
  57. }
  58. Vector3<Real> twoExtent = box1.extent * (Real)2;
  59. // The LCP has 6 variables and 6 (nontrivial) inequality
  60. // constraints.
  61. std::array<Real, 12> q =
  62. {
  63. -delta[0], -delta[1], -delta[2], rotDelta[0], rotDelta[1], rotDelta[2],
  64. K[0], K[1], K[2], twoExtent[0], twoExtent[1], twoExtent[2]
  65. };
  66. std::array<std::array<Real, 12>, 12> M;
  67. {
  68. Real const z = (Real)0;
  69. Real const p = (Real)1;
  70. Real const m = (Real)-1;
  71. Vector3<Real> const& U0 = box1.axis[0];
  72. Vector3<Real> const& U1 = box1.axis[1];
  73. Vector3<Real> const& U2 = box1.axis[2];
  74. M[0] = { p, z, z, -U0[0], -U1[0], -U1[2], p, z, z, z, z, z };
  75. M[1] = { z, p, z, -U0[1], -U1[1], -U1[1], z, p, z, z, z, z };
  76. M[2] = { z, z, p, -U0[2], -U1[2], -U1[2], z, z, p, z, z, z };
  77. M[3] = { -U0[0], -U0[1], -U0[2], p, z, z, z, z, z, p, z, z };
  78. M[4] = { -U1[0], -U1[1], -U1[2], z, p, z, z, z, z, z, p, z };
  79. M[5] = { -U2[0], -U2[1], -U2[2], z, z, p, z, z, z, z, z, p };
  80. M[6] = { m, z, z, z, z, z, z, z, z, z, z, z };
  81. M[7] = { z, m, z, z, z, z, z, z, z, z, z, z };
  82. M[8] = { z, z, m, z, z, z, z, z, z, z, z, z };
  83. M[9] = { z, z, z, m, z, z, z, z, z, z, z, z };
  84. M[10] = { z, z, z, z, m, z, z, z, z, z, z, z };
  85. M[11] = { z, z, z, z, z, m, z, z, z, z, z, z };
  86. }
  87. std::array<Real, 12> w, z;
  88. if (mLCP.Solve(q, M, w, z))
  89. {
  90. result.queryIsSuccessful = true;
  91. for (int i = 0; i < 3; ++i)
  92. {
  93. result.box0Parameter[i] = z[i] + box0.min[i];
  94. result.closestPoint[0][i] = result.box0Parameter[i];
  95. }
  96. result.closestPoint[1] = box1.center;
  97. for (int i = 0, j = 3; i < 3; ++i, ++j)
  98. {
  99. result.box1Parameter[i] = z[j] - box1.extent[i];
  100. result.closestPoint[1] += result.box1Parameter[i] * box1.axis[i];
  101. }
  102. Vector3<Real> diff = result.closestPoint[1] - result.closestPoint[0];
  103. result.sqrDistance = Dot(diff, diff);
  104. result.distance = std::sqrt(result.sqrDistance);
  105. }
  106. else
  107. {
  108. // If you reach this case, the maximum number of iterations
  109. // was not specified to be large enough or there is a problem
  110. // due to floating-point rounding errors. If you believe the
  111. // latter is true, file a bug report.
  112. result.queryIsSuccessful = false;
  113. for (int i = 0; i < 3; ++i)
  114. {
  115. result.box0Parameter[i] = (Real)0;
  116. result.box1Parameter[i] = (Real)0;
  117. result.closestPoint[0][i] = (Real)0;
  118. result.closestPoint[1][i] = (Real)0;
  119. }
  120. result.distance = (Real)0;
  121. result.sqrDistance = (Real)0;
  122. }
  123. result.numLCPIterations = mLCP.GetNumIterations();
  124. return result;
  125. }
  126. private:
  127. LCPSolver<Real, 12> mLCP;
  128. };
  129. }