DistLine3AlignedBox3.h 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514
  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/DistPointAlignedBox.h>
  9. #include <Mathematics/Line.h>
  10. #include <Mathematics/Vector3.h>
  11. namespace WwiseGTE
  12. {
  13. template <typename Real>
  14. class DCPQuery<Real, Line3<Real>, AlignedBox3<Real>>
  15. {
  16. public:
  17. struct Result
  18. {
  19. Real distance, sqrDistance;
  20. Real lineParameter;
  21. Vector3<Real> closestPoint[2];
  22. };
  23. Result operator()(Line3<Real> const& line, AlignedBox3<Real> const& box)
  24. {
  25. // Translate the line and box so that the box has center at the
  26. // origin.
  27. Vector3<Real> boxCenter, boxExtent;
  28. box.GetCenteredForm(boxCenter, boxExtent);
  29. Vector3<Real> point = line.origin - boxCenter;
  30. Vector3<Real> direction = line.direction;
  31. Result result;
  32. DoQuery(point, direction, boxExtent, result);
  33. // Compute the closest point on the line.
  34. result.closestPoint[0] = line.origin + result.lineParameter * line.direction;
  35. // Compute the closest point on the box.
  36. result.closestPoint[1] = boxCenter + point;
  37. return result;
  38. }
  39. protected:
  40. // Compute the distance and closest point between a line and an
  41. // axis-aligned box whose center is the origin. On input, 'point' is
  42. // the line origin and 'direction' is the line direction. On output,
  43. // 'point' is the point on the box closest to the line. The
  44. // 'direction' is non-const to allow transforming the problem into
  45. // the first octant.
  46. void DoQuery(Vector3<Real>& point, Vector3<Real>& direction,
  47. Vector3<Real> const& boxExtent, Result& result)
  48. {
  49. result.sqrDistance = (Real)0;
  50. result.lineParameter = (Real)0;
  51. // Apply reflections so that direction vector has nonnegative
  52. // components.
  53. bool reflect[3];
  54. for (int i = 0; i < 3; ++i)
  55. {
  56. if (direction[i] < (Real)0)
  57. {
  58. point[i] = -point[i];
  59. direction[i] = -direction[i];
  60. reflect[i] = true;
  61. }
  62. else
  63. {
  64. reflect[i] = false;
  65. }
  66. }
  67. if (direction[0] > (Real)0)
  68. {
  69. if (direction[1] > (Real)0)
  70. {
  71. if (direction[2] > (Real)0) // (+,+,+)
  72. {
  73. CaseNoZeros(point, direction, boxExtent, result);
  74. }
  75. else // (+,+,0)
  76. {
  77. Case0(0, 1, 2, point, direction, boxExtent, result);
  78. }
  79. }
  80. else
  81. {
  82. if (direction[2] > (Real)0) // (+,0,+)
  83. {
  84. Case0(0, 2, 1, point, direction, boxExtent, result);
  85. }
  86. else // (+,0,0)
  87. {
  88. Case00(0, 1, 2, point, direction, boxExtent, result);
  89. }
  90. }
  91. }
  92. else
  93. {
  94. if (direction[1] > (Real)0)
  95. {
  96. if (direction[2] > (Real)0) // (0,+,+)
  97. {
  98. Case0(1, 2, 0, point, direction, boxExtent, result);
  99. }
  100. else // (0,+,0)
  101. {
  102. Case00(1, 0, 2, point, direction, boxExtent, result);
  103. }
  104. }
  105. else
  106. {
  107. if (direction[2] > (Real)0) // (0,0,+)
  108. {
  109. Case00(2, 0, 1, point, direction, boxExtent, result);
  110. }
  111. else // (0,0,0)
  112. {
  113. Case000(point, boxExtent, result);
  114. }
  115. }
  116. }
  117. // Undo the reflections applied previously.
  118. for (int i = 0; i < 3; ++i)
  119. {
  120. if (reflect[i])
  121. {
  122. point[i] = -point[i];
  123. }
  124. }
  125. result.distance = std::sqrt(result.sqrDistance);
  126. }
  127. private:
  128. void Face(int i0, int i1, int i2, Vector3<Real>& pnt,
  129. Vector3<Real> const& dir, Vector3<Real> const& PmE,
  130. Vector3<Real> const& boxExtent, Result& result)
  131. {
  132. Vector3<Real> PpE;
  133. Real lenSqr, inv, tmp, param, t, delta;
  134. PpE[i1] = pnt[i1] + boxExtent[i1];
  135. PpE[i2] = pnt[i2] + boxExtent[i2];
  136. if (dir[i0] * PpE[i1] >= dir[i1] * PmE[i0])
  137. {
  138. if (dir[i0] * PpE[i2] >= dir[i2] * PmE[i0])
  139. {
  140. // v[i1] >= -e[i1], v[i2] >= -e[i2] (distance = 0)
  141. pnt[i0] = boxExtent[i0];
  142. inv = ((Real)1) / dir[i0];
  143. pnt[i1] -= dir[i1] * PmE[i0] * inv;
  144. pnt[i2] -= dir[i2] * PmE[i0] * inv;
  145. result.lineParameter = -PmE[i0] * inv;
  146. }
  147. else
  148. {
  149. // v[i1] >= -e[i1], v[i2] < -e[i2]
  150. lenSqr = dir[i0] * dir[i0] + dir[i2] * dir[i2];
  151. tmp = lenSqr * PpE[i1] - dir[i1] * (dir[i0] * PmE[i0] +
  152. dir[i2] * PpE[i2]);
  153. if (tmp <= ((Real)2) * lenSqr * boxExtent[i1])
  154. {
  155. t = tmp / lenSqr;
  156. lenSqr += dir[i1] * dir[i1];
  157. tmp = PpE[i1] - t;
  158. delta = dir[i0] * PmE[i0] + dir[i1] * tmp + dir[i2] * PpE[i2];
  159. param = -delta / lenSqr;
  160. result.sqrDistance += PmE[i0] * PmE[i0] + tmp * tmp +
  161. PpE[i2] * PpE[i2] + delta * param;
  162. result.lineParameter = param;
  163. pnt[i0] = boxExtent[i0];
  164. pnt[i1] = t - boxExtent[i1];
  165. pnt[i2] = -boxExtent[i2];
  166. }
  167. else
  168. {
  169. lenSqr += dir[i1] * dir[i1];
  170. delta = dir[i0] * PmE[i0] + dir[i1] * PmE[i1] + dir[i2] * PpE[i2];
  171. param = -delta / lenSqr;
  172. result.sqrDistance += PmE[i0] * PmE[i0] + PmE[i1] * PmE[i1]
  173. + PpE[i2] * PpE[i2] + delta * param;
  174. result.lineParameter = param;
  175. pnt[i0] = boxExtent[i0];
  176. pnt[i1] = boxExtent[i1];
  177. pnt[i2] = -boxExtent[i2];
  178. }
  179. }
  180. }
  181. else
  182. {
  183. if (dir[i0] * PpE[i2] >= dir[i2] * PmE[i0])
  184. {
  185. // v[i1] < -e[i1], v[i2] >= -e[i2]
  186. lenSqr = dir[i0] * dir[i0] + dir[i1] * dir[i1];
  187. tmp = lenSqr * PpE[i2] - dir[i2] * (dir[i0] * PmE[i0] +
  188. dir[i1] * PpE[i1]);
  189. if (tmp <= ((Real)2) * lenSqr * boxExtent[i2])
  190. {
  191. t = tmp / lenSqr;
  192. lenSqr += dir[i2] * dir[i2];
  193. tmp = PpE[i2] - t;
  194. delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * tmp;
  195. param = -delta / lenSqr;
  196. result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] +
  197. tmp * tmp + delta * param;
  198. result.lineParameter = param;
  199. pnt[i0] = boxExtent[i0];
  200. pnt[i1] = -boxExtent[i1];
  201. pnt[i2] = t - boxExtent[i2];
  202. }
  203. else
  204. {
  205. lenSqr += dir[i2] * dir[i2];
  206. delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * PmE[i2];
  207. param = -delta / lenSqr;
  208. result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] +
  209. PmE[i2] * PmE[i2] + delta * param;
  210. result.lineParameter = param;
  211. pnt[i0] = boxExtent[i0];
  212. pnt[i1] = -boxExtent[i1];
  213. pnt[i2] = boxExtent[i2];
  214. }
  215. }
  216. else
  217. {
  218. // v[i1] < -e[i1], v[i2] < -e[i2]
  219. lenSqr = dir[i0] * dir[i0] + dir[i2] * dir[i2];
  220. tmp = lenSqr * PpE[i1] - dir[i1] * (dir[i0] * PmE[i0] +
  221. dir[i2] * PpE[i2]);
  222. if (tmp >= (Real)0)
  223. {
  224. // v[i1]-edge is closest
  225. if (tmp <= ((Real)2) * lenSqr * boxExtent[i1])
  226. {
  227. t = tmp / lenSqr;
  228. lenSqr += dir[i1] * dir[i1];
  229. tmp = PpE[i1] - t;
  230. delta = dir[i0] * PmE[i0] + dir[i1] * tmp + dir[i2] * PpE[i2];
  231. param = -delta / lenSqr;
  232. result.sqrDistance += PmE[i0] * PmE[i0] + tmp * tmp +
  233. PpE[i2] * PpE[i2] + delta * param;
  234. result.lineParameter = param;
  235. pnt[i0] = boxExtent[i0];
  236. pnt[i1] = t - boxExtent[i1];
  237. pnt[i2] = -boxExtent[i2];
  238. }
  239. else
  240. {
  241. lenSqr += dir[i1] * dir[i1];
  242. delta = dir[i0] * PmE[i0] + dir[i1] * PmE[i1]
  243. + dir[i2] * PpE[i2];
  244. param = -delta / lenSqr;
  245. result.sqrDistance += PmE[i0] * PmE[i0] + PmE[i1] * PmE[i1]
  246. + PpE[i2] * PpE[i2] + delta * param;
  247. result.lineParameter = param;
  248. pnt[i0] = boxExtent[i0];
  249. pnt[i1] = boxExtent[i1];
  250. pnt[i2] = -boxExtent[i2];
  251. }
  252. return;
  253. }
  254. lenSqr = dir[i0] * dir[i0] + dir[i1] * dir[i1];
  255. tmp = lenSqr * PpE[i2] - dir[i2] * (dir[i0] * PmE[i0] +
  256. dir[i1] * PpE[i1]);
  257. if (tmp >= (Real)0)
  258. {
  259. // v[i2]-edge is closest
  260. if (tmp <= ((Real)2) * lenSqr * boxExtent[i2])
  261. {
  262. t = tmp / lenSqr;
  263. lenSqr += dir[i2] * dir[i2];
  264. tmp = PpE[i2] - t;
  265. delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * tmp;
  266. param = -delta / lenSqr;
  267. result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] +
  268. tmp * tmp + delta * param;
  269. result.lineParameter = param;
  270. pnt[i0] = boxExtent[i0];
  271. pnt[i1] = -boxExtent[i1];
  272. pnt[i2] = t - boxExtent[i2];
  273. }
  274. else
  275. {
  276. lenSqr += dir[i2] * dir[i2];
  277. delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1]
  278. + dir[i2] * PmE[i2];
  279. param = -delta / lenSqr;
  280. result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1]
  281. + PmE[i2] * PmE[i2] + delta * param;
  282. result.lineParameter = param;
  283. pnt[i0] = boxExtent[i0];
  284. pnt[i1] = -boxExtent[i1];
  285. pnt[i2] = boxExtent[i2];
  286. }
  287. return;
  288. }
  289. // (v[i1],v[i2])-corner is closest
  290. lenSqr += dir[i2] * dir[i2];
  291. delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * PpE[i2];
  292. param = -delta / lenSqr;
  293. result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1]
  294. + PpE[i2] * PpE[i2] + delta * param;
  295. result.lineParameter = param;
  296. pnt[i0] = boxExtent[i0];
  297. pnt[i1] = -boxExtent[i1];
  298. pnt[i2] = -boxExtent[i2];
  299. }
  300. }
  301. }
  302. void CaseNoZeros(Vector3<Real>& pnt, Vector3<Real> const& dir,
  303. Vector3<Real> const& boxExtent, Result& result)
  304. {
  305. Vector3<Real> PmE = pnt - boxExtent;
  306. Real prodDxPy = dir[0] * PmE[1];
  307. Real prodDyPx = dir[1] * PmE[0];
  308. Real prodDzPx, prodDxPz, prodDzPy, prodDyPz;
  309. if (prodDyPx >= prodDxPy)
  310. {
  311. prodDzPx = dir[2] * PmE[0];
  312. prodDxPz = dir[0] * PmE[2];
  313. if (prodDzPx >= prodDxPz)
  314. {
  315. // line intersects x = e0
  316. Face(0, 1, 2, pnt, dir, PmE, boxExtent, result);
  317. }
  318. else
  319. {
  320. // line intersects z = e2
  321. Face(2, 0, 1, pnt, dir, PmE, boxExtent, result);
  322. }
  323. }
  324. else
  325. {
  326. prodDzPy = dir[2] * PmE[1];
  327. prodDyPz = dir[1] * PmE[2];
  328. if (prodDzPy >= prodDyPz)
  329. {
  330. // line intersects y = e1
  331. Face(1, 2, 0, pnt, dir, PmE, boxExtent, result);
  332. }
  333. else
  334. {
  335. // line intersects z = e2
  336. Face(2, 0, 1, pnt, dir, PmE, boxExtent, result);
  337. }
  338. }
  339. }
  340. void Case0(int i0, int i1, int i2, Vector3<Real>& pnt,
  341. Vector3<Real> const& dir, Vector3<Real> const& boxExtent, Result& result)
  342. {
  343. Real PmE0 = pnt[i0] - boxExtent[i0];
  344. Real PmE1 = pnt[i1] - boxExtent[i1];
  345. Real prod0 = dir[i1] * PmE0;
  346. Real prod1 = dir[i0] * PmE1;
  347. Real delta, invLSqr, inv;
  348. if (prod0 >= prod1)
  349. {
  350. // line intersects P[i0] = e[i0]
  351. pnt[i0] = boxExtent[i0];
  352. Real PpE1 = pnt[i1] + boxExtent[i1];
  353. delta = prod0 - dir[i0] * PpE1;
  354. if (delta >= (Real)0)
  355. {
  356. invLSqr = ((Real)1) / (dir[i0] * dir[i0] + dir[i1] * dir[i1]);
  357. result.sqrDistance += delta * delta * invLSqr;
  358. pnt[i1] = -boxExtent[i1];
  359. result.lineParameter = -(dir[i0] * PmE0 + dir[i1] * PpE1) * invLSqr;
  360. }
  361. else
  362. {
  363. inv = ((Real)1) / dir[i0];
  364. pnt[i1] -= prod0 * inv;
  365. result.lineParameter = -PmE0 * inv;
  366. }
  367. }
  368. else
  369. {
  370. // line intersects P[i1] = e[i1]
  371. pnt[i1] = boxExtent[i1];
  372. Real PpE0 = pnt[i0] + boxExtent[i0];
  373. delta = prod1 - dir[i1] * PpE0;
  374. if (delta >= (Real)0)
  375. {
  376. invLSqr = ((Real)1) / (dir[i0] * dir[i0] + dir[i1] * dir[i1]);
  377. result.sqrDistance += delta * delta * invLSqr;
  378. pnt[i0] = -boxExtent[i0];
  379. result.lineParameter = -(dir[i0] * PpE0 + dir[i1] * PmE1) * invLSqr;
  380. }
  381. else
  382. {
  383. inv = ((Real)1) / dir[i1];
  384. pnt[i0] -= prod1 * inv;
  385. result.lineParameter = -PmE1 * inv;
  386. }
  387. }
  388. if (pnt[i2] < -boxExtent[i2])
  389. {
  390. delta = pnt[i2] + boxExtent[i2];
  391. result.sqrDistance += delta * delta;
  392. pnt[i2] = -boxExtent[i2];
  393. }
  394. else if (pnt[i2] > boxExtent[i2])
  395. {
  396. delta = pnt[i2] - boxExtent[i2];
  397. result.sqrDistance += delta * delta;
  398. pnt[i2] = boxExtent[i2];
  399. }
  400. }
  401. void Case00(int i0, int i1, int i2, Vector3<Real>& pnt,
  402. Vector3<Real> const& dir, Vector3<Real> const& boxExtent, Result& result)
  403. {
  404. Real delta;
  405. result.lineParameter = (boxExtent[i0] - pnt[i0]) / dir[i0];
  406. pnt[i0] = boxExtent[i0];
  407. if (pnt[i1] < -boxExtent[i1])
  408. {
  409. delta = pnt[i1] + boxExtent[i1];
  410. result.sqrDistance += delta * delta;
  411. pnt[i1] = -boxExtent[i1];
  412. }
  413. else if (pnt[i1] > boxExtent[i1])
  414. {
  415. delta = pnt[i1] - boxExtent[i1];
  416. result.sqrDistance += delta * delta;
  417. pnt[i1] = boxExtent[i1];
  418. }
  419. if (pnt[i2] < -boxExtent[i2])
  420. {
  421. delta = pnt[i2] + boxExtent[i2];
  422. result.sqrDistance += delta * delta;
  423. pnt[i2] = -boxExtent[i2];
  424. }
  425. else if (pnt[i2] > boxExtent[i2])
  426. {
  427. delta = pnt[i2] - boxExtent[i2];
  428. result.sqrDistance += delta * delta;
  429. pnt[i2] = boxExtent[i2];
  430. }
  431. }
  432. void Case000(Vector3<Real>& pnt, Vector3<Real> const& boxExtent, Result& result)
  433. {
  434. Real delta;
  435. if (pnt[0] < -boxExtent[0])
  436. {
  437. delta = pnt[0] + boxExtent[0];
  438. result.sqrDistance += delta * delta;
  439. pnt[0] = -boxExtent[0];
  440. }
  441. else if (pnt[0] > boxExtent[0])
  442. {
  443. delta = pnt[0] - boxExtent[0];
  444. result.sqrDistance += delta * delta;
  445. pnt[0] = boxExtent[0];
  446. }
  447. if (pnt[1] < -boxExtent[1])
  448. {
  449. delta = pnt[1] + boxExtent[1];
  450. result.sqrDistance += delta * delta;
  451. pnt[1] = -boxExtent[1];
  452. }
  453. else if (pnt[1] > boxExtent[1])
  454. {
  455. delta = pnt[1] - boxExtent[1];
  456. result.sqrDistance += delta * delta;
  457. pnt[1] = boxExtent[1];
  458. }
  459. if (pnt[2] < -boxExtent[2])
  460. {
  461. delta = pnt[2] + boxExtent[2];
  462. result.sqrDistance += delta * delta;
  463. pnt[2] = -boxExtent[2];
  464. }
  465. else if (pnt[2] > boxExtent[2])
  466. {
  467. delta = pnt[2] - boxExtent[2];
  468. result.sqrDistance += delta * delta;
  469. pnt[2] = boxExtent[2];
  470. }
  471. }
  472. };
  473. }