IntersectionTests.js 32 KB

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