IntersectionTests.js 29 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047
  1. import Cartesian3 from "./Cartesian3.js";
  2. import Cartographic from "./Cartographic.js";
  3. import defaultValue from "./defaultValue.js";
  4. import defined from "./defined.js";
  5. import DeveloperError from "./DeveloperError.js";
  6. import Interval from "./Interval.js";
  7. import CesiumMath from "./Math.js";
  8. import Matrix3 from "./Matrix3.js";
  9. import QuadraticRealPolynomial from "./QuadraticRealPolynomial.js";
  10. import QuarticRealPolynomial from "./QuarticRealPolynomial.js";
  11. import Ray from "./Ray.js";
  12. /**
  13. * Functions for computing the intersection between geometries such as rays, planes, triangles, and ellipsoids.
  14. *
  15. * @namespace IntersectionTests
  16. */
  17. const IntersectionTests = {};
  18. /**
  19. * Computes the intersection of a ray and a plane.
  20. *
  21. * @param {Ray} ray The ray.
  22. * @param {Plane} plane The plane.
  23. * @param {Cartesian3} [result] The object onto which to store the result.
  24. * @returns {Cartesian3} The intersection point or undefined if there is no intersections.
  25. */
  26. IntersectionTests.rayPlane = function (ray, plane, result) {
  27. //>>includeStart('debug', pragmas.debug);
  28. if (!defined(ray)) {
  29. throw new DeveloperError("ray is required.");
  30. }
  31. if (!defined(plane)) {
  32. throw new DeveloperError("plane is required.");
  33. }
  34. //>>includeEnd('debug');
  35. if (!defined(result)) {
  36. result = new Cartesian3();
  37. }
  38. const origin = ray.origin;
  39. const direction = ray.direction;
  40. const normal = plane.normal;
  41. const denominator = Cartesian3.dot(normal, direction);
  42. if (Math.abs(denominator) < CesiumMath.EPSILON15) {
  43. // Ray is parallel to plane. The ray may be in the polygon's plane.
  44. return undefined;
  45. }
  46. const t = (-plane.distance - Cartesian3.dot(normal, origin)) / denominator;
  47. if (t < 0) {
  48. return undefined;
  49. }
  50. result = Cartesian3.multiplyByScalar(direction, t, result);
  51. return Cartesian3.add(origin, result, result);
  52. };
  53. const scratchEdge0 = new Cartesian3();
  54. const scratchEdge1 = new Cartesian3();
  55. const scratchPVec = new Cartesian3();
  56. const scratchTVec = new Cartesian3();
  57. const scratchQVec = new Cartesian3();
  58. /**
  59. * Computes the intersection of a ray and a triangle as a parametric distance along the input ray. The result is negative when the triangle is behind the ray.
  60. *
  61. * Implements {@link https://cadxfem.org/inf/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf|
  62. * Fast Minimum Storage Ray/Triangle Intersection} by Tomas Moller and Ben Trumbore.
  63. *
  64. * @memberof IntersectionTests
  65. *
  66. * @param {Ray} ray The ray.
  67. * @param {Cartesian3} p0 The first vertex of the triangle.
  68. * @param {Cartesian3} p1 The second vertex of the triangle.
  69. * @param {Cartesian3} p2 The third vertex of the triangle.
  70. * @param {Boolean} [cullBackFaces=false] If <code>true</code>, will only compute an intersection with the front face of the triangle
  71. * and return undefined for intersections with the back face.
  72. * @returns {Number} The intersection as a parametric distance along the ray, or undefined if there is no intersection.
  73. */
  74. IntersectionTests.rayTriangleParametric = function (
  75. ray,
  76. p0,
  77. p1,
  78. p2,
  79. cullBackFaces
  80. ) {
  81. //>>includeStart('debug', pragmas.debug);
  82. if (!defined(ray)) {
  83. throw new DeveloperError("ray is required.");
  84. }
  85. if (!defined(p0)) {
  86. throw new DeveloperError("p0 is required.");
  87. }
  88. if (!defined(p1)) {
  89. throw new DeveloperError("p1 is required.");
  90. }
  91. if (!defined(p2)) {
  92. throw new DeveloperError("p2 is required.");
  93. }
  94. //>>includeEnd('debug');
  95. cullBackFaces = defaultValue(cullBackFaces, false);
  96. const origin = ray.origin;
  97. const direction = ray.direction;
  98. const edge0 = Cartesian3.subtract(p1, p0, scratchEdge0);
  99. const edge1 = Cartesian3.subtract(p2, p0, scratchEdge1);
  100. const p = Cartesian3.cross(direction, edge1, scratchPVec);
  101. const det = Cartesian3.dot(edge0, p);
  102. let tvec;
  103. let q;
  104. let u;
  105. let v;
  106. let t;
  107. if (cullBackFaces) {
  108. if (det < CesiumMath.EPSILON6) {
  109. return undefined;
  110. }
  111. tvec = Cartesian3.subtract(origin, p0, scratchTVec);
  112. u = Cartesian3.dot(tvec, p);
  113. if (u < 0.0 || u > det) {
  114. return undefined;
  115. }
  116. q = Cartesian3.cross(tvec, edge0, scratchQVec);
  117. v = Cartesian3.dot(direction, q);
  118. if (v < 0.0 || u + v > det) {
  119. return undefined;
  120. }
  121. t = Cartesian3.dot(edge1, q) / det;
  122. } else {
  123. if (Math.abs(det) < CesiumMath.EPSILON6) {
  124. return undefined;
  125. }
  126. const invDet = 1.0 / det;
  127. tvec = Cartesian3.subtract(origin, p0, scratchTVec);
  128. u = Cartesian3.dot(tvec, p) * invDet;
  129. if (u < 0.0 || u > 1.0) {
  130. return undefined;
  131. }
  132. q = Cartesian3.cross(tvec, edge0, scratchQVec);
  133. v = Cartesian3.dot(direction, q) * invDet;
  134. if (v < 0.0 || u + v > 1.0) {
  135. return undefined;
  136. }
  137. t = Cartesian3.dot(edge1, q) * invDet;
  138. }
  139. return t;
  140. };
  141. /**
  142. * Computes the intersection of a ray and a triangle as a Cartesian3 coordinate.
  143. *
  144. * Implements {@link https://cadxfem.org/inf/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf|
  145. * Fast Minimum Storage Ray/Triangle Intersection} by Tomas Moller and Ben Trumbore.
  146. *
  147. * @memberof IntersectionTests
  148. *
  149. * @param {Ray} ray The ray.
  150. * @param {Cartesian3} p0 The first vertex of the triangle.
  151. * @param {Cartesian3} p1 The second vertex of the triangle.
  152. * @param {Cartesian3} p2 The third vertex of the triangle.
  153. * @param {Boolean} [cullBackFaces=false] If <code>true</code>, will only compute an intersection with the front face of the triangle
  154. * and return undefined for intersections with the back face.
  155. * @param {Cartesian3} [result] The <code>Cartesian3</code> onto which to store the result.
  156. * @returns {Cartesian3} The intersection point or undefined if there is no intersections.
  157. */
  158. IntersectionTests.rayTriangle = function (
  159. ray,
  160. p0,
  161. p1,
  162. p2,
  163. cullBackFaces,
  164. result
  165. ) {
  166. const t = IntersectionTests.rayTriangleParametric(
  167. ray,
  168. p0,
  169. p1,
  170. p2,
  171. cullBackFaces
  172. );
  173. if (!defined(t) || t < 0.0) {
  174. return undefined;
  175. }
  176. if (!defined(result)) {
  177. result = new Cartesian3();
  178. }
  179. Cartesian3.multiplyByScalar(ray.direction, t, result);
  180. return Cartesian3.add(ray.origin, result, result);
  181. };
  182. const scratchLineSegmentTriangleRay = new Ray();
  183. /**
  184. * Computes the intersection of a line segment and a triangle.
  185. * @memberof IntersectionTests
  186. *
  187. * @param {Cartesian3} v0 The an end point of the line segment.
  188. * @param {Cartesian3} v1 The other end point of the line segment.
  189. * @param {Cartesian3} p0 The first vertex of the triangle.
  190. * @param {Cartesian3} p1 The second vertex of the triangle.
  191. * @param {Cartesian3} p2 The third vertex of the triangle.
  192. * @param {Boolean} [cullBackFaces=false] If <code>true</code>, will only compute an intersection with the front face of the triangle
  193. * and return undefined for intersections with the back face.
  194. * @param {Cartesian3} [result] The <code>Cartesian3</code> onto which to store the result.
  195. * @returns {Cartesian3} The intersection point or undefined if there is no intersections.
  196. */
  197. IntersectionTests.lineSegmentTriangle = function (
  198. v0,
  199. v1,
  200. p0,
  201. p1,
  202. p2,
  203. cullBackFaces,
  204. result
  205. ) {
  206. //>>includeStart('debug', pragmas.debug);
  207. if (!defined(v0)) {
  208. throw new DeveloperError("v0 is required.");
  209. }
  210. if (!defined(v1)) {
  211. throw new DeveloperError("v1 is required.");
  212. }
  213. if (!defined(p0)) {
  214. throw new DeveloperError("p0 is required.");
  215. }
  216. if (!defined(p1)) {
  217. throw new DeveloperError("p1 is required.");
  218. }
  219. if (!defined(p2)) {
  220. throw new DeveloperError("p2 is required.");
  221. }
  222. //>>includeEnd('debug');
  223. const ray = scratchLineSegmentTriangleRay;
  224. Cartesian3.clone(v0, ray.origin);
  225. Cartesian3.subtract(v1, v0, ray.direction);
  226. Cartesian3.normalize(ray.direction, ray.direction);
  227. const t = IntersectionTests.rayTriangleParametric(
  228. ray,
  229. p0,
  230. p1,
  231. p2,
  232. cullBackFaces
  233. );
  234. if (!defined(t) || t < 0.0 || t > Cartesian3.distance(v0, v1)) {
  235. return undefined;
  236. }
  237. if (!defined(result)) {
  238. result = new Cartesian3();
  239. }
  240. Cartesian3.multiplyByScalar(ray.direction, t, result);
  241. return Cartesian3.add(ray.origin, result, result);
  242. };
  243. function solveQuadratic(a, b, c, result) {
  244. const det = b * b - 4.0 * a * c;
  245. if (det < 0.0) {
  246. return undefined;
  247. } else if (det > 0.0) {
  248. const denom = 1.0 / (2.0 * a);
  249. const disc = Math.sqrt(det);
  250. const root0 = (-b + disc) * denom;
  251. const root1 = (-b - disc) * denom;
  252. if (root0 < root1) {
  253. result.root0 = root0;
  254. result.root1 = root1;
  255. } else {
  256. result.root0 = root1;
  257. result.root1 = root0;
  258. }
  259. return result;
  260. }
  261. const root = -b / (2.0 * a);
  262. if (root === 0.0) {
  263. return undefined;
  264. }
  265. result.root0 = result.root1 = root;
  266. return result;
  267. }
  268. const raySphereRoots = {
  269. root0: 0.0,
  270. root1: 0.0,
  271. };
  272. function raySphere(ray, sphere, result) {
  273. if (!defined(result)) {
  274. result = new Interval();
  275. }
  276. const origin = ray.origin;
  277. const direction = ray.direction;
  278. const center = sphere.center;
  279. const radiusSquared = sphere.radius * sphere.radius;
  280. const diff = Cartesian3.subtract(origin, center, scratchPVec);
  281. const a = Cartesian3.dot(direction, direction);
  282. const b = 2.0 * Cartesian3.dot(direction, diff);
  283. const c = Cartesian3.magnitudeSquared(diff) - radiusSquared;
  284. const roots = solveQuadratic(a, b, c, raySphereRoots);
  285. if (!defined(roots)) {
  286. return undefined;
  287. }
  288. result.start = roots.root0;
  289. result.stop = roots.root1;
  290. return result;
  291. }
  292. /**
  293. * Computes the intersection points of a ray with a sphere.
  294. * @memberof IntersectionTests
  295. *
  296. * @param {Ray} ray The ray.
  297. * @param {BoundingSphere} sphere The sphere.
  298. * @param {Interval} [result] The result onto which to store the result.
  299. * @returns {Interval} The interval containing scalar points along the ray or undefined if there are no intersections.
  300. */
  301. IntersectionTests.raySphere = function (ray, sphere, result) {
  302. //>>includeStart('debug', pragmas.debug);
  303. if (!defined(ray)) {
  304. throw new DeveloperError("ray is required.");
  305. }
  306. if (!defined(sphere)) {
  307. throw new DeveloperError("sphere is required.");
  308. }
  309. //>>includeEnd('debug');
  310. result = raySphere(ray, sphere, result);
  311. if (!defined(result) || result.stop < 0.0) {
  312. return undefined;
  313. }
  314. result.start = Math.max(result.start, 0.0);
  315. return result;
  316. };
  317. const scratchLineSegmentRay = new Ray();
  318. /**
  319. * Computes the intersection points of a line segment with a sphere.
  320. * @memberof IntersectionTests
  321. *
  322. * @param {Cartesian3} p0 An end point of the line segment.
  323. * @param {Cartesian3} p1 The other end point of the line segment.
  324. * @param {BoundingSphere} sphere The sphere.
  325. * @param {Interval} [result] The result onto which to store the result.
  326. * @returns {Interval} The interval containing scalar points along the ray or undefined if there are no intersections.
  327. */
  328. IntersectionTests.lineSegmentSphere = function (p0, p1, sphere, result) {
  329. //>>includeStart('debug', pragmas.debug);
  330. if (!defined(p0)) {
  331. throw new DeveloperError("p0 is required.");
  332. }
  333. if (!defined(p1)) {
  334. throw new DeveloperError("p1 is required.");
  335. }
  336. if (!defined(sphere)) {
  337. throw new DeveloperError("sphere is required.");
  338. }
  339. //>>includeEnd('debug');
  340. const ray = scratchLineSegmentRay;
  341. Cartesian3.clone(p0, ray.origin);
  342. const direction = Cartesian3.subtract(p1, p0, ray.direction);
  343. const maxT = Cartesian3.magnitude(direction);
  344. Cartesian3.normalize(direction, direction);
  345. result = raySphere(ray, sphere, result);
  346. if (!defined(result) || result.stop < 0.0 || result.start > maxT) {
  347. return undefined;
  348. }
  349. result.start = Math.max(result.start, 0.0);
  350. result.stop = Math.min(result.stop, maxT);
  351. return result;
  352. };
  353. const scratchQ = new Cartesian3();
  354. const scratchW = new Cartesian3();
  355. /**
  356. * Computes the intersection points of a ray with an ellipsoid.
  357. *
  358. * @param {Ray} ray The ray.
  359. * @param {Ellipsoid} ellipsoid The ellipsoid.
  360. * @returns {Interval} The interval containing scalar points along the ray or undefined if there are no intersections.
  361. */
  362. IntersectionTests.rayEllipsoid = function (ray, ellipsoid) {
  363. //>>includeStart('debug', pragmas.debug);
  364. if (!defined(ray)) {
  365. throw new DeveloperError("ray is required.");
  366. }
  367. if (!defined(ellipsoid)) {
  368. throw new DeveloperError("ellipsoid is required.");
  369. }
  370. //>>includeEnd('debug');
  371. const inverseRadii = ellipsoid.oneOverRadii;
  372. const q = Cartesian3.multiplyComponents(inverseRadii, ray.origin, scratchQ);
  373. const w = Cartesian3.multiplyComponents(
  374. inverseRadii,
  375. ray.direction,
  376. scratchW
  377. );
  378. const q2 = Cartesian3.magnitudeSquared(q);
  379. const qw = Cartesian3.dot(q, w);
  380. let difference, w2, product, discriminant, temp;
  381. if (q2 > 1.0) {
  382. // Outside ellipsoid.
  383. if (qw >= 0.0) {
  384. // Looking outward or tangent (0 intersections).
  385. return undefined;
  386. }
  387. // qw < 0.0.
  388. const qw2 = qw * qw;
  389. difference = q2 - 1.0; // Positively valued.
  390. w2 = Cartesian3.magnitudeSquared(w);
  391. product = w2 * difference;
  392. if (qw2 < product) {
  393. // Imaginary roots (0 intersections).
  394. return undefined;
  395. } else if (qw2 > product) {
  396. // Distinct roots (2 intersections).
  397. discriminant = qw * qw - product;
  398. temp = -qw + Math.sqrt(discriminant); // Avoid cancellation.
  399. const root0 = temp / w2;
  400. const root1 = difference / temp;
  401. if (root0 < root1) {
  402. return new Interval(root0, root1);
  403. }
  404. return {
  405. start: root1,
  406. stop: root0,
  407. };
  408. }
  409. // qw2 == product. Repeated roots (2 intersections).
  410. const root = Math.sqrt(difference / w2);
  411. return new Interval(root, root);
  412. } else if (q2 < 1.0) {
  413. // Inside ellipsoid (2 intersections).
  414. difference = q2 - 1.0; // Negatively valued.
  415. w2 = Cartesian3.magnitudeSquared(w);
  416. product = w2 * difference; // Negatively valued.
  417. discriminant = qw * qw - product;
  418. temp = -qw + Math.sqrt(discriminant); // Positively valued.
  419. return new Interval(0.0, temp / w2);
  420. }
  421. // q2 == 1.0. On ellipsoid.
  422. if (qw < 0.0) {
  423. // Looking inward.
  424. w2 = Cartesian3.magnitudeSquared(w);
  425. return new Interval(0.0, -qw / w2);
  426. }
  427. // qw >= 0.0. Looking outward or tangent.
  428. return undefined;
  429. };
  430. function addWithCancellationCheck(left, right, tolerance) {
  431. const difference = left + right;
  432. if (
  433. CesiumMath.sign(left) !== CesiumMath.sign(right) &&
  434. Math.abs(difference / Math.max(Math.abs(left), Math.abs(right))) < tolerance
  435. ) {
  436. return 0.0;
  437. }
  438. return difference;
  439. }
  440. function quadraticVectorExpression(A, b, c, x, w) {
  441. const xSquared = x * x;
  442. const wSquared = w * w;
  443. const l2 = (A[Matrix3.COLUMN1ROW1] - A[Matrix3.COLUMN2ROW2]) * wSquared;
  444. const l1 =
  445. w *
  446. (x *
  447. addWithCancellationCheck(
  448. A[Matrix3.COLUMN1ROW0],
  449. A[Matrix3.COLUMN0ROW1],
  450. CesiumMath.EPSILON15
  451. ) +
  452. b.y);
  453. const l0 =
  454. A[Matrix3.COLUMN0ROW0] * xSquared +
  455. A[Matrix3.COLUMN2ROW2] * wSquared +
  456. x * b.x +
  457. c;
  458. const r1 =
  459. wSquared *
  460. addWithCancellationCheck(
  461. A[Matrix3.COLUMN2ROW1],
  462. A[Matrix3.COLUMN1ROW2],
  463. CesiumMath.EPSILON15
  464. );
  465. const r0 =
  466. w *
  467. (x *
  468. addWithCancellationCheck(A[Matrix3.COLUMN2ROW0], A[Matrix3.COLUMN0ROW2]) +
  469. b.z);
  470. let cosines;
  471. const solutions = [];
  472. if (r0 === 0.0 && r1 === 0.0) {
  473. cosines = QuadraticRealPolynomial.computeRealRoots(l2, l1, l0);
  474. if (cosines.length === 0) {
  475. return solutions;
  476. }
  477. const cosine0 = cosines[0];
  478. const sine0 = Math.sqrt(Math.max(1.0 - cosine0 * cosine0, 0.0));
  479. solutions.push(new Cartesian3(x, w * cosine0, w * -sine0));
  480. solutions.push(new Cartesian3(x, w * cosine0, w * sine0));
  481. if (cosines.length === 2) {
  482. const cosine1 = cosines[1];
  483. const sine1 = Math.sqrt(Math.max(1.0 - cosine1 * cosine1, 0.0));
  484. solutions.push(new Cartesian3(x, w * cosine1, w * -sine1));
  485. solutions.push(new Cartesian3(x, w * cosine1, w * sine1));
  486. }
  487. return solutions;
  488. }
  489. const r0Squared = r0 * r0;
  490. const r1Squared = r1 * r1;
  491. const l2Squared = l2 * l2;
  492. const r0r1 = r0 * r1;
  493. const c4 = l2Squared + r1Squared;
  494. const c3 = 2.0 * (l1 * l2 + r0r1);
  495. const c2 = 2.0 * l0 * l2 + l1 * l1 - r1Squared + r0Squared;
  496. const c1 = 2.0 * (l0 * l1 - r0r1);
  497. const c0 = l0 * l0 - r0Squared;
  498. if (c4 === 0.0 && c3 === 0.0 && c2 === 0.0 && c1 === 0.0) {
  499. return solutions;
  500. }
  501. cosines = QuarticRealPolynomial.computeRealRoots(c4, c3, c2, c1, c0);
  502. const length = cosines.length;
  503. if (length === 0) {
  504. return solutions;
  505. }
  506. for (let i = 0; i < length; ++i) {
  507. const cosine = cosines[i];
  508. const cosineSquared = cosine * cosine;
  509. const sineSquared = Math.max(1.0 - cosineSquared, 0.0);
  510. const sine = Math.sqrt(sineSquared);
  511. //const left = l2 * cosineSquared + l1 * cosine + l0;
  512. let left;
  513. if (CesiumMath.sign(l2) === CesiumMath.sign(l0)) {
  514. left = addWithCancellationCheck(
  515. l2 * cosineSquared + l0,
  516. l1 * cosine,
  517. CesiumMath.EPSILON12
  518. );
  519. } else if (CesiumMath.sign(l0) === CesiumMath.sign(l1 * cosine)) {
  520. left = addWithCancellationCheck(
  521. l2 * cosineSquared,
  522. l1 * cosine + l0,
  523. CesiumMath.EPSILON12
  524. );
  525. } else {
  526. left = addWithCancellationCheck(
  527. l2 * cosineSquared + l1 * cosine,
  528. l0,
  529. CesiumMath.EPSILON12
  530. );
  531. }
  532. const right = addWithCancellationCheck(
  533. r1 * cosine,
  534. r0,
  535. CesiumMath.EPSILON15
  536. );
  537. const product = left * right;
  538. if (product < 0.0) {
  539. solutions.push(new Cartesian3(x, w * cosine, w * sine));
  540. } else if (product > 0.0) {
  541. solutions.push(new Cartesian3(x, w * cosine, w * -sine));
  542. } else if (sine !== 0.0) {
  543. solutions.push(new Cartesian3(x, w * cosine, w * -sine));
  544. solutions.push(new Cartesian3(x, w * cosine, w * sine));
  545. ++i;
  546. } else {
  547. solutions.push(new Cartesian3(x, w * cosine, w * sine));
  548. }
  549. }
  550. return solutions;
  551. }
  552. const firstAxisScratch = new Cartesian3();
  553. const secondAxisScratch = new Cartesian3();
  554. const thirdAxisScratch = new Cartesian3();
  555. const referenceScratch = new Cartesian3();
  556. const bCart = new Cartesian3();
  557. const bScratch = new Matrix3();
  558. const btScratch = new Matrix3();
  559. const diScratch = new Matrix3();
  560. const dScratch = new Matrix3();
  561. const cScratch = new Matrix3();
  562. const tempMatrix = new Matrix3();
  563. const aScratch = new Matrix3();
  564. const sScratch = new Cartesian3();
  565. const closestScratch = new Cartesian3();
  566. const surfPointScratch = new Cartographic();
  567. /**
  568. * Provides the point along the ray which is nearest to the ellipsoid.
  569. *
  570. * @param {Ray} ray The ray.
  571. * @param {Ellipsoid} ellipsoid The ellipsoid.
  572. * @returns {Cartesian3} The nearest planetodetic point on the ray.
  573. */
  574. IntersectionTests.grazingAltitudeLocation = function (ray, ellipsoid) {
  575. //>>includeStart('debug', pragmas.debug);
  576. if (!defined(ray)) {
  577. throw new DeveloperError("ray is required.");
  578. }
  579. if (!defined(ellipsoid)) {
  580. throw new DeveloperError("ellipsoid is required.");
  581. }
  582. //>>includeEnd('debug');
  583. const position = ray.origin;
  584. const direction = ray.direction;
  585. if (!Cartesian3.equals(position, Cartesian3.ZERO)) {
  586. const normal = ellipsoid.geodeticSurfaceNormal(position, firstAxisScratch);
  587. if (Cartesian3.dot(direction, normal) >= 0.0) {
  588. // The location provided is the closest point in altitude
  589. return position;
  590. }
  591. }
  592. const intersects = defined(this.rayEllipsoid(ray, ellipsoid));
  593. // Compute the scaled direction vector.
  594. const f = ellipsoid.transformPositionToScaledSpace(
  595. direction,
  596. firstAxisScratch
  597. );
  598. // Constructs a basis from the unit scaled direction vector. Construct its rotation and transpose.
  599. const firstAxis = Cartesian3.normalize(f, f);
  600. const reference = Cartesian3.mostOrthogonalAxis(f, referenceScratch);
  601. const secondAxis = Cartesian3.normalize(
  602. Cartesian3.cross(reference, firstAxis, secondAxisScratch),
  603. secondAxisScratch
  604. );
  605. const thirdAxis = Cartesian3.normalize(
  606. Cartesian3.cross(firstAxis, secondAxis, thirdAxisScratch),
  607. thirdAxisScratch
  608. );
  609. const B = bScratch;
  610. B[0] = firstAxis.x;
  611. B[1] = firstAxis.y;
  612. B[2] = firstAxis.z;
  613. B[3] = secondAxis.x;
  614. B[4] = secondAxis.y;
  615. B[5] = secondAxis.z;
  616. B[6] = thirdAxis.x;
  617. B[7] = thirdAxis.y;
  618. B[8] = thirdAxis.z;
  619. const B_T = Matrix3.transpose(B, btScratch);
  620. // Get the scaling matrix and its inverse.
  621. const D_I = Matrix3.fromScale(ellipsoid.radii, diScratch);
  622. const D = Matrix3.fromScale(ellipsoid.oneOverRadii, dScratch);
  623. const C = cScratch;
  624. C[0] = 0.0;
  625. C[1] = -direction.z;
  626. C[2] = direction.y;
  627. C[3] = direction.z;
  628. C[4] = 0.0;
  629. C[5] = -direction.x;
  630. C[6] = -direction.y;
  631. C[7] = direction.x;
  632. C[8] = 0.0;
  633. const temp = Matrix3.multiply(
  634. Matrix3.multiply(B_T, D, tempMatrix),
  635. C,
  636. tempMatrix
  637. );
  638. const A = Matrix3.multiply(
  639. Matrix3.multiply(temp, D_I, aScratch),
  640. B,
  641. aScratch
  642. );
  643. const b = Matrix3.multiplyByVector(temp, position, bCart);
  644. // Solve for the solutions to the expression in standard form:
  645. const solutions = quadraticVectorExpression(
  646. A,
  647. Cartesian3.negate(b, firstAxisScratch),
  648. 0.0,
  649. 0.0,
  650. 1.0
  651. );
  652. let s;
  653. let altitude;
  654. const length = solutions.length;
  655. if (length > 0) {
  656. let closest = Cartesian3.clone(Cartesian3.ZERO, closestScratch);
  657. let maximumValue = Number.NEGATIVE_INFINITY;
  658. for (let i = 0; i < length; ++i) {
  659. s = Matrix3.multiplyByVector(
  660. D_I,
  661. Matrix3.multiplyByVector(B, solutions[i], sScratch),
  662. sScratch
  663. );
  664. const v = Cartesian3.normalize(
  665. Cartesian3.subtract(s, position, referenceScratch),
  666. referenceScratch
  667. );
  668. const dotProduct = Cartesian3.dot(v, direction);
  669. if (dotProduct > maximumValue) {
  670. maximumValue = dotProduct;
  671. closest = Cartesian3.clone(s, closest);
  672. }
  673. }
  674. const surfacePoint = ellipsoid.cartesianToCartographic(
  675. closest,
  676. surfPointScratch
  677. );
  678. maximumValue = CesiumMath.clamp(maximumValue, 0.0, 1.0);
  679. altitude =
  680. Cartesian3.magnitude(
  681. Cartesian3.subtract(closest, position, referenceScratch)
  682. ) * Math.sqrt(1.0 - maximumValue * maximumValue);
  683. altitude = intersects ? -altitude : altitude;
  684. surfacePoint.height = altitude;
  685. return ellipsoid.cartographicToCartesian(surfacePoint, new Cartesian3());
  686. }
  687. return undefined;
  688. };
  689. const lineSegmentPlaneDifference = new Cartesian3();
  690. /**
  691. * Computes the intersection of a line segment and a plane.
  692. *
  693. * @param {Cartesian3} endPoint0 An end point of the line segment.
  694. * @param {Cartesian3} endPoint1 The other end point of the line segment.
  695. * @param {Plane} plane The plane.
  696. * @param {Cartesian3} [result] The object onto which to store the result.
  697. * @returns {Cartesian3} The intersection point or undefined if there is no intersection.
  698. *
  699. * @example
  700. * const origin = Cesium.Cartesian3.fromDegrees(-75.59777, 40.03883);
  701. * const normal = ellipsoid.geodeticSurfaceNormal(origin);
  702. * const plane = Cesium.Plane.fromPointNormal(origin, normal);
  703. *
  704. * const p0 = new Cesium.Cartesian3(...);
  705. * const p1 = new Cesium.Cartesian3(...);
  706. *
  707. * // find the intersection of the line segment from p0 to p1 and the tangent plane at origin.
  708. * const intersection = Cesium.IntersectionTests.lineSegmentPlane(p0, p1, plane);
  709. */
  710. IntersectionTests.lineSegmentPlane = function (
  711. endPoint0,
  712. endPoint1,
  713. plane,
  714. result
  715. ) {
  716. //>>includeStart('debug', pragmas.debug);
  717. if (!defined(endPoint0)) {
  718. throw new DeveloperError("endPoint0 is required.");
  719. }
  720. if (!defined(endPoint1)) {
  721. throw new DeveloperError("endPoint1 is required.");
  722. }
  723. if (!defined(plane)) {
  724. throw new DeveloperError("plane is required.");
  725. }
  726. //>>includeEnd('debug');
  727. if (!defined(result)) {
  728. result = new Cartesian3();
  729. }
  730. const difference = Cartesian3.subtract(
  731. endPoint1,
  732. endPoint0,
  733. lineSegmentPlaneDifference
  734. );
  735. const normal = plane.normal;
  736. const nDotDiff = Cartesian3.dot(normal, difference);
  737. // check if the segment and plane are parallel
  738. if (Math.abs(nDotDiff) < CesiumMath.EPSILON6) {
  739. return undefined;
  740. }
  741. const nDotP0 = Cartesian3.dot(normal, endPoint0);
  742. const t = -(plane.distance + nDotP0) / nDotDiff;
  743. // intersection only if t is in [0, 1]
  744. if (t < 0.0 || t > 1.0) {
  745. return undefined;
  746. }
  747. // intersection is endPoint0 + t * (endPoint1 - endPoint0)
  748. Cartesian3.multiplyByScalar(difference, t, result);
  749. Cartesian3.add(endPoint0, result, result);
  750. return result;
  751. };
  752. /**
  753. * Computes the intersection of a triangle and a plane
  754. *
  755. * @param {Cartesian3} p0 First point of the triangle
  756. * @param {Cartesian3} p1 Second point of the triangle
  757. * @param {Cartesian3} p2 Third point of the triangle
  758. * @param {Plane} plane Intersection plane
  759. * @returns {Object} An object with properties <code>positions</code> and <code>indices</code>, which are arrays that represent three triangles that do not cross the plane. (Undefined if no intersection exists)
  760. *
  761. * @example
  762. * const origin = Cesium.Cartesian3.fromDegrees(-75.59777, 40.03883);
  763. * const normal = ellipsoid.geodeticSurfaceNormal(origin);
  764. * const plane = Cesium.Plane.fromPointNormal(origin, normal);
  765. *
  766. * const p0 = new Cesium.Cartesian3(...);
  767. * const p1 = new Cesium.Cartesian3(...);
  768. * const p2 = new Cesium.Cartesian3(...);
  769. *
  770. * // convert the triangle composed of points (p0, p1, p2) to three triangles that don't cross the plane
  771. * const triangles = Cesium.IntersectionTests.trianglePlaneIntersection(p0, p1, p2, plane);
  772. */
  773. IntersectionTests.trianglePlaneIntersection = function (p0, p1, p2, plane) {
  774. //>>includeStart('debug', pragmas.debug);
  775. if (!defined(p0) || !defined(p1) || !defined(p2) || !defined(plane)) {
  776. throw new DeveloperError("p0, p1, p2, and plane are required.");
  777. }
  778. //>>includeEnd('debug');
  779. const planeNormal = plane.normal;
  780. const planeD = plane.distance;
  781. const p0Behind = Cartesian3.dot(planeNormal, p0) + planeD < 0.0;
  782. const p1Behind = Cartesian3.dot(planeNormal, p1) + planeD < 0.0;
  783. const p2Behind = Cartesian3.dot(planeNormal, p2) + planeD < 0.0;
  784. // Given these dots products, the calls to lineSegmentPlaneIntersection
  785. // always have defined results.
  786. let numBehind = 0;
  787. numBehind += p0Behind ? 1 : 0;
  788. numBehind += p1Behind ? 1 : 0;
  789. numBehind += p2Behind ? 1 : 0;
  790. let u1, u2;
  791. if (numBehind === 1 || numBehind === 2) {
  792. u1 = new Cartesian3();
  793. u2 = new Cartesian3();
  794. }
  795. if (numBehind === 1) {
  796. if (p0Behind) {
  797. IntersectionTests.lineSegmentPlane(p0, p1, plane, u1);
  798. IntersectionTests.lineSegmentPlane(p0, p2, plane, u2);
  799. return {
  800. positions: [p0, p1, p2, u1, u2],
  801. indices: [
  802. // Behind
  803. 0,
  804. 3,
  805. 4,
  806. // In front
  807. 1,
  808. 2,
  809. 4,
  810. 1,
  811. 4,
  812. 3,
  813. ],
  814. };
  815. } else if (p1Behind) {
  816. IntersectionTests.lineSegmentPlane(p1, p2, plane, u1);
  817. IntersectionTests.lineSegmentPlane(p1, p0, plane, u2);
  818. return {
  819. positions: [p0, p1, p2, u1, u2],
  820. indices: [
  821. // Behind
  822. 1,
  823. 3,
  824. 4,
  825. // In front
  826. 2,
  827. 0,
  828. 4,
  829. 2,
  830. 4,
  831. 3,
  832. ],
  833. };
  834. } else if (p2Behind) {
  835. IntersectionTests.lineSegmentPlane(p2, p0, plane, u1);
  836. IntersectionTests.lineSegmentPlane(p2, p1, plane, u2);
  837. return {
  838. positions: [p0, p1, p2, u1, u2],
  839. indices: [
  840. // Behind
  841. 2,
  842. 3,
  843. 4,
  844. // In front
  845. 0,
  846. 1,
  847. 4,
  848. 0,
  849. 4,
  850. 3,
  851. ],
  852. };
  853. }
  854. } else if (numBehind === 2) {
  855. if (!p0Behind) {
  856. IntersectionTests.lineSegmentPlane(p1, p0, plane, u1);
  857. IntersectionTests.lineSegmentPlane(p2, p0, plane, u2);
  858. return {
  859. positions: [p0, p1, p2, u1, u2],
  860. indices: [
  861. // Behind
  862. 1,
  863. 2,
  864. 4,
  865. 1,
  866. 4,
  867. 3,
  868. // In front
  869. 0,
  870. 3,
  871. 4,
  872. ],
  873. };
  874. } else if (!p1Behind) {
  875. IntersectionTests.lineSegmentPlane(p2, p1, plane, u1);
  876. IntersectionTests.lineSegmentPlane(p0, p1, plane, u2);
  877. return {
  878. positions: [p0, p1, p2, u1, u2],
  879. indices: [
  880. // Behind
  881. 2,
  882. 0,
  883. 4,
  884. 2,
  885. 4,
  886. 3,
  887. // In front
  888. 1,
  889. 3,
  890. 4,
  891. ],
  892. };
  893. } else if (!p2Behind) {
  894. IntersectionTests.lineSegmentPlane(p0, p2, plane, u1);
  895. IntersectionTests.lineSegmentPlane(p1, p2, plane, u2);
  896. return {
  897. positions: [p0, p1, p2, u1, u2],
  898. indices: [
  899. // Behind
  900. 0,
  901. 1,
  902. 4,
  903. 0,
  904. 4,
  905. 3,
  906. // In front
  907. 2,
  908. 3,
  909. 4,
  910. ],
  911. };
  912. }
  913. }
  914. // if numBehind is 3, the triangle is completely behind the plane;
  915. // otherwise, it is completely in front (numBehind is 0).
  916. return undefined;
  917. };
  918. export default IntersectionTests;