IntersectionTests-ca40c01c.js 67 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652
  1. /**
  2. * Cesium - https://github.com/CesiumGS/cesium
  3. *
  4. * Copyright 2011-2020 Cesium Contributors
  5. *
  6. * Licensed under the Apache License, Version 2.0 (the "License");
  7. * you may not use this file except in compliance with the License.
  8. * You may obtain a copy of the License at
  9. *
  10. * http://www.apache.org/licenses/LICENSE-2.0
  11. *
  12. * Unless required by applicable law or agreed to in writing, software
  13. * distributed under the License is distributed on an "AS IS" BASIS,
  14. * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  15. * See the License for the specific language governing permissions and
  16. * limitations under the License.
  17. *
  18. * Columbus View (Pat. Pend.)
  19. *
  20. * Portions licensed separately.
  21. * See https://github.com/CesiumGS/cesium/blob/master/LICENSE.md for full licensing details.
  22. */
  23. define(['exports', './when-8d13db60', './Check-70bec281', './Math-61ede240', './Cartographic-fe4be337', './BoundingSphere-8f8a682c'], function (exports, when, Check, _Math, Cartographic, BoundingSphere) { 'use strict';
  24. /**
  25. * Defines functions for 2nd order polynomial functions of one variable with only real coefficients.
  26. *
  27. * @exports QuadraticRealPolynomial
  28. */
  29. var QuadraticRealPolynomial = {};
  30. /**
  31. * Provides the discriminant of the quadratic equation from the supplied coefficients.
  32. *
  33. * @param {Number} a The coefficient of the 2nd order monomial.
  34. * @param {Number} b The coefficient of the 1st order monomial.
  35. * @param {Number} c The coefficient of the 0th order monomial.
  36. * @returns {Number} The value of the discriminant.
  37. */
  38. QuadraticRealPolynomial.computeDiscriminant = function(a, b, c) {
  39. //>>includeStart('debug', pragmas.debug);
  40. if (typeof a !== 'number') {
  41. throw new Check.DeveloperError('a is a required number.');
  42. }
  43. if (typeof b !== 'number') {
  44. throw new Check.DeveloperError('b is a required number.');
  45. }
  46. if (typeof c !== 'number') {
  47. throw new Check.DeveloperError('c is a required number.');
  48. }
  49. //>>includeEnd('debug');
  50. var discriminant = b * b - 4.0 * a * c;
  51. return discriminant;
  52. };
  53. function addWithCancellationCheck(left, right, tolerance) {
  54. var difference = left + right;
  55. if ((_Math.CesiumMath.sign(left) !== _Math.CesiumMath.sign(right)) &&
  56. Math.abs(difference / Math.max(Math.abs(left), Math.abs(right))) < tolerance) {
  57. return 0.0;
  58. }
  59. return difference;
  60. }
  61. /**
  62. * Provides the real valued roots of the quadratic polynomial with the provided coefficients.
  63. *
  64. * @param {Number} a The coefficient of the 2nd order monomial.
  65. * @param {Number} b The coefficient of the 1st order monomial.
  66. * @param {Number} c The coefficient of the 0th order monomial.
  67. * @returns {Number[]} The real valued roots.
  68. */
  69. QuadraticRealPolynomial.computeRealRoots = function(a, b, c) {
  70. //>>includeStart('debug', pragmas.debug);
  71. if (typeof a !== 'number') {
  72. throw new Check.DeveloperError('a is a required number.');
  73. }
  74. if (typeof b !== 'number') {
  75. throw new Check.DeveloperError('b is a required number.');
  76. }
  77. if (typeof c !== 'number') {
  78. throw new Check.DeveloperError('c is a required number.');
  79. }
  80. //>>includeEnd('debug');
  81. var ratio;
  82. if (a === 0.0) {
  83. if (b === 0.0) {
  84. // Constant function: c = 0.
  85. return [];
  86. }
  87. // Linear function: b * x + c = 0.
  88. return [-c / b];
  89. } else if (b === 0.0) {
  90. if (c === 0.0) {
  91. // 2nd order monomial: a * x^2 = 0.
  92. return [0.0, 0.0];
  93. }
  94. var cMagnitude = Math.abs(c);
  95. var aMagnitude = Math.abs(a);
  96. if ((cMagnitude < aMagnitude) && (cMagnitude / aMagnitude < _Math.CesiumMath.EPSILON14)) { // c ~= 0.0.
  97. // 2nd order monomial: a * x^2 = 0.
  98. return [0.0, 0.0];
  99. } else if ((cMagnitude > aMagnitude) && (aMagnitude / cMagnitude < _Math.CesiumMath.EPSILON14)) { // a ~= 0.0.
  100. // Constant function: c = 0.
  101. return [];
  102. }
  103. // a * x^2 + c = 0
  104. ratio = -c / a;
  105. if (ratio < 0.0) {
  106. // Both roots are complex.
  107. return [];
  108. }
  109. // Both roots are real.
  110. var root = Math.sqrt(ratio);
  111. return [-root, root];
  112. } else if (c === 0.0) {
  113. // a * x^2 + b * x = 0
  114. ratio = -b / a;
  115. if (ratio < 0.0) {
  116. return [ratio, 0.0];
  117. }
  118. return [0.0, ratio];
  119. }
  120. // a * x^2 + b * x + c = 0
  121. var b2 = b * b;
  122. var four_ac = 4.0 * a * c;
  123. var radicand = addWithCancellationCheck(b2, -four_ac, _Math.CesiumMath.EPSILON14);
  124. if (radicand < 0.0) {
  125. // Both roots are complex.
  126. return [];
  127. }
  128. var q = -0.5 * addWithCancellationCheck(b, _Math.CesiumMath.sign(b) * Math.sqrt(radicand), _Math.CesiumMath.EPSILON14);
  129. if (b > 0.0) {
  130. return [q / a, c / q];
  131. }
  132. return [c / q, q / a];
  133. };
  134. /**
  135. * Defines functions for 3rd order polynomial functions of one variable with only real coefficients.
  136. *
  137. * @exports CubicRealPolynomial
  138. */
  139. var CubicRealPolynomial = {};
  140. /**
  141. * Provides the discriminant of the cubic equation from the supplied coefficients.
  142. *
  143. * @param {Number} a The coefficient of the 3rd order monomial.
  144. * @param {Number} b The coefficient of the 2nd order monomial.
  145. * @param {Number} c The coefficient of the 1st order monomial.
  146. * @param {Number} d The coefficient of the 0th order monomial.
  147. * @returns {Number} The value of the discriminant.
  148. */
  149. CubicRealPolynomial.computeDiscriminant = function(a, b, c, d) {
  150. //>>includeStart('debug', pragmas.debug);
  151. if (typeof a !== 'number') {
  152. throw new Check.DeveloperError('a is a required number.');
  153. }
  154. if (typeof b !== 'number') {
  155. throw new Check.DeveloperError('b is a required number.');
  156. }
  157. if (typeof c !== 'number') {
  158. throw new Check.DeveloperError('c is a required number.');
  159. }
  160. if (typeof d !== 'number') {
  161. throw new Check.DeveloperError('d is a required number.');
  162. }
  163. //>>includeEnd('debug');
  164. var a2 = a * a;
  165. var b2 = b * b;
  166. var c2 = c * c;
  167. var d2 = d * d;
  168. var discriminant = 18.0 * a * b * c * d + b2 * c2 - 27.0 * a2 * d2 - 4.0 * (a * c2 * c + b2 * b * d);
  169. return discriminant;
  170. };
  171. function computeRealRoots(a, b, c, d) {
  172. var A = a;
  173. var B = b / 3.0;
  174. var C = c / 3.0;
  175. var D = d;
  176. var AC = A * C;
  177. var BD = B * D;
  178. var B2 = B * B;
  179. var C2 = C * C;
  180. var delta1 = A * C - B2;
  181. var delta2 = A * D - B * C;
  182. var delta3 = B * D - C2;
  183. var discriminant = 4.0 * delta1 * delta3 - delta2 * delta2;
  184. var temp;
  185. var temp1;
  186. if (discriminant < 0.0) {
  187. var ABar;
  188. var CBar;
  189. var DBar;
  190. if (B2 * BD >= AC * C2) {
  191. ABar = A;
  192. CBar = delta1;
  193. DBar = -2.0 * B * delta1 + A * delta2;
  194. } else {
  195. ABar = D;
  196. CBar = delta3;
  197. DBar = -D * delta2 + 2.0 * C * delta3;
  198. }
  199. var s = (DBar < 0.0) ? -1.0 : 1.0; // This is not Math.Sign()!
  200. var temp0 = -s * Math.abs(ABar) * Math.sqrt(-discriminant);
  201. temp1 = -DBar + temp0;
  202. var x = temp1 / 2.0;
  203. var p = x < 0.0 ? -Math.pow(-x, 1.0 / 3.0) : Math.pow(x, 1.0 / 3.0);
  204. var q = (temp1 === temp0) ? -p : -CBar / p;
  205. temp = (CBar <= 0.0) ? p + q : -DBar / (p * p + q * q + CBar);
  206. if (B2 * BD >= AC * C2) {
  207. return [(temp - B) / A];
  208. }
  209. return [-D / (temp + C)];
  210. }
  211. var CBarA = delta1;
  212. var DBarA = -2.0 * B * delta1 + A * delta2;
  213. var CBarD = delta3;
  214. var DBarD = -D * delta2 + 2.0 * C * delta3;
  215. var squareRootOfDiscriminant = Math.sqrt(discriminant);
  216. var halfSquareRootOf3 = Math.sqrt(3.0) / 2.0;
  217. var theta = Math.abs(Math.atan2(A * squareRootOfDiscriminant, -DBarA) / 3.0);
  218. temp = 2.0 * Math.sqrt(-CBarA);
  219. var cosine = Math.cos(theta);
  220. temp1 = temp * cosine;
  221. var temp3 = temp * (-cosine / 2.0 - halfSquareRootOf3 * Math.sin(theta));
  222. var numeratorLarge = (temp1 + temp3 > 2.0 * B) ? temp1 - B : temp3 - B;
  223. var denominatorLarge = A;
  224. var root1 = numeratorLarge / denominatorLarge;
  225. theta = Math.abs(Math.atan2(D * squareRootOfDiscriminant, -DBarD) / 3.0);
  226. temp = 2.0 * Math.sqrt(-CBarD);
  227. cosine = Math.cos(theta);
  228. temp1 = temp * cosine;
  229. temp3 = temp * (-cosine / 2.0 - halfSquareRootOf3 * Math.sin(theta));
  230. var numeratorSmall = -D;
  231. var denominatorSmall = (temp1 + temp3 < 2.0 * C) ? temp1 + C : temp3 + C;
  232. var root3 = numeratorSmall / denominatorSmall;
  233. var E = denominatorLarge * denominatorSmall;
  234. var F = -numeratorLarge * denominatorSmall - denominatorLarge * numeratorSmall;
  235. var G = numeratorLarge * numeratorSmall;
  236. var root2 = (C * F - B * G) / (-B * F + C * E);
  237. if (root1 <= root2) {
  238. if (root1 <= root3) {
  239. if (root2 <= root3) {
  240. return [root1, root2, root3];
  241. }
  242. return [root1, root3, root2];
  243. }
  244. return [root3, root1, root2];
  245. }
  246. if (root1 <= root3) {
  247. return [root2, root1, root3];
  248. }
  249. if (root2 <= root3) {
  250. return [root2, root3, root1];
  251. }
  252. return [root3, root2, root1];
  253. }
  254. /**
  255. * Provides the real valued roots of the cubic polynomial with the provided coefficients.
  256. *
  257. * @param {Number} a The coefficient of the 3rd order monomial.
  258. * @param {Number} b The coefficient of the 2nd order monomial.
  259. * @param {Number} c The coefficient of the 1st order monomial.
  260. * @param {Number} d The coefficient of the 0th order monomial.
  261. * @returns {Number[]} The real valued roots.
  262. */
  263. CubicRealPolynomial.computeRealRoots = function(a, b, c, d) {
  264. //>>includeStart('debug', pragmas.debug);
  265. if (typeof a !== 'number') {
  266. throw new Check.DeveloperError('a is a required number.');
  267. }
  268. if (typeof b !== 'number') {
  269. throw new Check.DeveloperError('b is a required number.');
  270. }
  271. if (typeof c !== 'number') {
  272. throw new Check.DeveloperError('c is a required number.');
  273. }
  274. if (typeof d !== 'number') {
  275. throw new Check.DeveloperError('d is a required number.');
  276. }
  277. //>>includeEnd('debug');
  278. var roots;
  279. var ratio;
  280. if (a === 0.0) {
  281. // Quadratic function: b * x^2 + c * x + d = 0.
  282. return QuadraticRealPolynomial.computeRealRoots(b, c, d);
  283. } else if (b === 0.0) {
  284. if (c === 0.0) {
  285. if (d === 0.0) {
  286. // 3rd order monomial: a * x^3 = 0.
  287. return [0.0, 0.0, 0.0];
  288. }
  289. // a * x^3 + d = 0
  290. ratio = -d / a;
  291. var root = (ratio < 0.0) ? -Math.pow(-ratio, 1.0 / 3.0) : Math.pow(ratio, 1.0 / 3.0);
  292. return [root, root, root];
  293. } else if (d === 0.0) {
  294. // x * (a * x^2 + c) = 0.
  295. roots = QuadraticRealPolynomial.computeRealRoots(a, 0, c);
  296. // Return the roots in ascending order.
  297. if (roots.Length === 0) {
  298. return [0.0];
  299. }
  300. return [roots[0], 0.0, roots[1]];
  301. }
  302. // Deflated cubic polynomial: a * x^3 + c * x + d= 0.
  303. return computeRealRoots(a, 0, c, d);
  304. } else if (c === 0.0) {
  305. if (d === 0.0) {
  306. // x^2 * (a * x + b) = 0.
  307. ratio = -b / a;
  308. if (ratio < 0.0) {
  309. return [ratio, 0.0, 0.0];
  310. }
  311. return [0.0, 0.0, ratio];
  312. }
  313. // a * x^3 + b * x^2 + d = 0.
  314. return computeRealRoots(a, b, 0, d);
  315. } else if (d === 0.0) {
  316. // x * (a * x^2 + b * x + c) = 0
  317. roots = QuadraticRealPolynomial.computeRealRoots(a, b, c);
  318. // Return the roots in ascending order.
  319. if (roots.length === 0) {
  320. return [0.0];
  321. } else if (roots[1] <= 0.0) {
  322. return [roots[0], roots[1], 0.0];
  323. } else if (roots[0] >= 0.0) {
  324. return [0.0, roots[0], roots[1]];
  325. }
  326. return [roots[0], 0.0, roots[1]];
  327. }
  328. return computeRealRoots(a, b, c, d);
  329. };
  330. /**
  331. * Defines functions for 4th order polynomial functions of one variable with only real coefficients.
  332. *
  333. * @exports QuarticRealPolynomial
  334. */
  335. var QuarticRealPolynomial = {};
  336. /**
  337. * Provides the discriminant of the quartic equation from the supplied coefficients.
  338. *
  339. * @param {Number} a The coefficient of the 4th order monomial.
  340. * @param {Number} b The coefficient of the 3rd order monomial.
  341. * @param {Number} c The coefficient of the 2nd order monomial.
  342. * @param {Number} d The coefficient of the 1st order monomial.
  343. * @param {Number} e The coefficient of the 0th order monomial.
  344. * @returns {Number} The value of the discriminant.
  345. */
  346. QuarticRealPolynomial.computeDiscriminant = function(a, b, c, d, e) {
  347. //>>includeStart('debug', pragmas.debug);
  348. if (typeof a !== 'number') {
  349. throw new Check.DeveloperError('a is a required number.');
  350. }
  351. if (typeof b !== 'number') {
  352. throw new Check.DeveloperError('b is a required number.');
  353. }
  354. if (typeof c !== 'number') {
  355. throw new Check.DeveloperError('c is a required number.');
  356. }
  357. if (typeof d !== 'number') {
  358. throw new Check.DeveloperError('d is a required number.');
  359. }
  360. if (typeof e !== 'number') {
  361. throw new Check.DeveloperError('e is a required number.');
  362. }
  363. //>>includeEnd('debug');
  364. var a2 = a * a;
  365. var a3 = a2 * a;
  366. var b2 = b * b;
  367. var b3 = b2 * b;
  368. var c2 = c * c;
  369. var c3 = c2 * c;
  370. var d2 = d * d;
  371. var d3 = d2 * d;
  372. var e2 = e * e;
  373. var e3 = e2 * e;
  374. var discriminant = (b2 * c2 * d2 - 4.0 * b3 * d3 - 4.0 * a * c3 * d2 + 18 * a * b * c * d3 - 27.0 * a2 * d2 * d2 + 256.0 * a3 * e3) +
  375. e * (18.0 * b3 * c * d - 4.0 * b2 * c3 + 16.0 * a * c2 * c2 - 80.0 * a * b * c2 * d - 6.0 * a * b2 * d2 + 144.0 * a2 * c * d2) +
  376. e2 * (144.0 * a * b2 * c - 27.0 * b2 * b2 - 128.0 * a2 * c2 - 192.0 * a2 * b * d);
  377. return discriminant;
  378. };
  379. function original(a3, a2, a1, a0) {
  380. var a3Squared = a3 * a3;
  381. var p = a2 - 3.0 * a3Squared / 8.0;
  382. var q = a1 - a2 * a3 / 2.0 + a3Squared * a3 / 8.0;
  383. var r = a0 - a1 * a3 / 4.0 + a2 * a3Squared / 16.0 - 3.0 * a3Squared * a3Squared / 256.0;
  384. // Find the roots of the cubic equations: h^6 + 2 p h^4 + (p^2 - 4 r) h^2 - q^2 = 0.
  385. var cubicRoots = CubicRealPolynomial.computeRealRoots(1.0, 2.0 * p, p * p - 4.0 * r, -q * q);
  386. if (cubicRoots.length > 0) {
  387. var temp = -a3 / 4.0;
  388. // Use the largest positive root.
  389. var hSquared = cubicRoots[cubicRoots.length - 1];
  390. if (Math.abs(hSquared) < _Math.CesiumMath.EPSILON14) {
  391. // y^4 + p y^2 + r = 0.
  392. var roots = QuadraticRealPolynomial.computeRealRoots(1.0, p, r);
  393. if (roots.length === 2) {
  394. var root0 = roots[0];
  395. var root1 = roots[1];
  396. var y;
  397. if (root0 >= 0.0 && root1 >= 0.0) {
  398. var y0 = Math.sqrt(root0);
  399. var y1 = Math.sqrt(root1);
  400. return [temp - y1, temp - y0, temp + y0, temp + y1];
  401. } else if (root0 >= 0.0 && root1 < 0.0) {
  402. y = Math.sqrt(root0);
  403. return [temp - y, temp + y];
  404. } else if (root0 < 0.0 && root1 >= 0.0) {
  405. y = Math.sqrt(root1);
  406. return [temp - y, temp + y];
  407. }
  408. }
  409. return [];
  410. } else if (hSquared > 0.0) {
  411. var h = Math.sqrt(hSquared);
  412. var m = (p + hSquared - q / h) / 2.0;
  413. var n = (p + hSquared + q / h) / 2.0;
  414. // Now solve the two quadratic factors: (y^2 + h y + m)(y^2 - h y + n);
  415. var roots1 = QuadraticRealPolynomial.computeRealRoots(1.0, h, m);
  416. var roots2 = QuadraticRealPolynomial.computeRealRoots(1.0, -h, n);
  417. if (roots1.length !== 0) {
  418. roots1[0] += temp;
  419. roots1[1] += temp;
  420. if (roots2.length !== 0) {
  421. roots2[0] += temp;
  422. roots2[1] += temp;
  423. if (roots1[1] <= roots2[0]) {
  424. return [roots1[0], roots1[1], roots2[0], roots2[1]];
  425. } else if (roots2[1] <= roots1[0]) {
  426. return [roots2[0], roots2[1], roots1[0], roots1[1]];
  427. } else if (roots1[0] >= roots2[0] && roots1[1] <= roots2[1]) {
  428. return [roots2[0], roots1[0], roots1[1], roots2[1]];
  429. } else if (roots2[0] >= roots1[0] && roots2[1] <= roots1[1]) {
  430. return [roots1[0], roots2[0], roots2[1], roots1[1]];
  431. } else if (roots1[0] > roots2[0] && roots1[0] < roots2[1]) {
  432. return [roots2[0], roots1[0], roots2[1], roots1[1]];
  433. }
  434. return [roots1[0], roots2[0], roots1[1], roots2[1]];
  435. }
  436. return roots1;
  437. }
  438. if (roots2.length !== 0) {
  439. roots2[0] += temp;
  440. roots2[1] += temp;
  441. return roots2;
  442. }
  443. return [];
  444. }
  445. }
  446. return [];
  447. }
  448. function neumark(a3, a2, a1, a0) {
  449. var a1Squared = a1 * a1;
  450. var a2Squared = a2 * a2;
  451. var a3Squared = a3 * a3;
  452. var p = -2.0 * a2;
  453. var q = a1 * a3 + a2Squared - 4.0 * a0;
  454. var r = a3Squared * a0 - a1 * a2 * a3 + a1Squared;
  455. var cubicRoots = CubicRealPolynomial.computeRealRoots(1.0, p, q, r);
  456. if (cubicRoots.length > 0) {
  457. // Use the most positive root
  458. var y = cubicRoots[0];
  459. var temp = (a2 - y);
  460. var tempSquared = temp * temp;
  461. var g1 = a3 / 2.0;
  462. var h1 = temp / 2.0;
  463. var m = tempSquared - 4.0 * a0;
  464. var mError = tempSquared + 4.0 * Math.abs(a0);
  465. var n = a3Squared - 4.0 * y;
  466. var nError = a3Squared + 4.0 * Math.abs(y);
  467. var g2;
  468. var h2;
  469. if (y < 0.0 || (m * nError < n * mError)) {
  470. var squareRootOfN = Math.sqrt(n);
  471. g2 = squareRootOfN / 2.0;
  472. h2 = squareRootOfN === 0.0 ? 0.0 : (a3 * h1 - a1) / squareRootOfN;
  473. } else {
  474. var squareRootOfM = Math.sqrt(m);
  475. g2 = squareRootOfM === 0.0 ? 0.0 : (a3 * h1 - a1) / squareRootOfM;
  476. h2 = squareRootOfM / 2.0;
  477. }
  478. var G;
  479. var g;
  480. if (g1 === 0.0 && g2 === 0.0) {
  481. G = 0.0;
  482. g = 0.0;
  483. } else if (_Math.CesiumMath.sign(g1) === _Math.CesiumMath.sign(g2)) {
  484. G = g1 + g2;
  485. g = y / G;
  486. } else {
  487. g = g1 - g2;
  488. G = y / g;
  489. }
  490. var H;
  491. var h;
  492. if (h1 === 0.0 && h2 === 0.0) {
  493. H = 0.0;
  494. h = 0.0;
  495. } else if (_Math.CesiumMath.sign(h1) === _Math.CesiumMath.sign(h2)) {
  496. H = h1 + h2;
  497. h = a0 / H;
  498. } else {
  499. h = h1 - h2;
  500. H = a0 / h;
  501. }
  502. // Now solve the two quadratic factors: (y^2 + G y + H)(y^2 + g y + h);
  503. var roots1 = QuadraticRealPolynomial.computeRealRoots(1.0, G, H);
  504. var roots2 = QuadraticRealPolynomial.computeRealRoots(1.0, g, h);
  505. if (roots1.length !== 0) {
  506. if (roots2.length !== 0) {
  507. if (roots1[1] <= roots2[0]) {
  508. return [roots1[0], roots1[1], roots2[0], roots2[1]];
  509. } else if (roots2[1] <= roots1[0]) {
  510. return [roots2[0], roots2[1], roots1[0], roots1[1]];
  511. } else if (roots1[0] >= roots2[0] && roots1[1] <= roots2[1]) {
  512. return [roots2[0], roots1[0], roots1[1], roots2[1]];
  513. } else if (roots2[0] >= roots1[0] && roots2[1] <= roots1[1]) {
  514. return [roots1[0], roots2[0], roots2[1], roots1[1]];
  515. } else if (roots1[0] > roots2[0] && roots1[0] < roots2[1]) {
  516. return [roots2[0], roots1[0], roots2[1], roots1[1]];
  517. }
  518. return [roots1[0], roots2[0], roots1[1], roots2[1]];
  519. }
  520. return roots1;
  521. }
  522. if (roots2.length !== 0) {
  523. return roots2;
  524. }
  525. }
  526. return [];
  527. }
  528. /**
  529. * Provides the real valued roots of the quartic polynomial with the provided coefficients.
  530. *
  531. * @param {Number} a The coefficient of the 4th order monomial.
  532. * @param {Number} b The coefficient of the 3rd order monomial.
  533. * @param {Number} c The coefficient of the 2nd order monomial.
  534. * @param {Number} d The coefficient of the 1st order monomial.
  535. * @param {Number} e The coefficient of the 0th order monomial.
  536. * @returns {Number[]} The real valued roots.
  537. */
  538. QuarticRealPolynomial.computeRealRoots = function(a, b, c, d, e) {
  539. //>>includeStart('debug', pragmas.debug);
  540. if (typeof a !== 'number') {
  541. throw new Check.DeveloperError('a is a required number.');
  542. }
  543. if (typeof b !== 'number') {
  544. throw new Check.DeveloperError('b is a required number.');
  545. }
  546. if (typeof c !== 'number') {
  547. throw new Check.DeveloperError('c is a required number.');
  548. }
  549. if (typeof d !== 'number') {
  550. throw new Check.DeveloperError('d is a required number.');
  551. }
  552. if (typeof e !== 'number') {
  553. throw new Check.DeveloperError('e is a required number.');
  554. }
  555. //>>includeEnd('debug');
  556. if (Math.abs(a) < _Math.CesiumMath.EPSILON15) {
  557. return CubicRealPolynomial.computeRealRoots(b, c, d, e);
  558. }
  559. var a3 = b / a;
  560. var a2 = c / a;
  561. var a1 = d / a;
  562. var a0 = e / a;
  563. var k = (a3 < 0.0) ? 1 : 0;
  564. k += (a2 < 0.0) ? k + 1 : k;
  565. k += (a1 < 0.0) ? k + 1 : k;
  566. k += (a0 < 0.0) ? k + 1 : k;
  567. switch (k) {
  568. case 0:
  569. return original(a3, a2, a1, a0);
  570. case 1:
  571. return neumark(a3, a2, a1, a0);
  572. case 2:
  573. return neumark(a3, a2, a1, a0);
  574. case 3:
  575. return original(a3, a2, a1, a0);
  576. case 4:
  577. return original(a3, a2, a1, a0);
  578. case 5:
  579. return neumark(a3, a2, a1, a0);
  580. case 6:
  581. return original(a3, a2, a1, a0);
  582. case 7:
  583. return original(a3, a2, a1, a0);
  584. case 8:
  585. return neumark(a3, a2, a1, a0);
  586. case 9:
  587. return original(a3, a2, a1, a0);
  588. case 10:
  589. return original(a3, a2, a1, a0);
  590. case 11:
  591. return neumark(a3, a2, a1, a0);
  592. case 12:
  593. return original(a3, a2, a1, a0);
  594. case 13:
  595. return original(a3, a2, a1, a0);
  596. case 14:
  597. return original(a3, a2, a1, a0);
  598. case 15:
  599. return original(a3, a2, a1, a0);
  600. default:
  601. return undefined;
  602. }
  603. };
  604. /**
  605. * Represents a ray that extends infinitely from the provided origin in the provided direction.
  606. * @alias Ray
  607. * @constructor
  608. *
  609. * @param {Cartesian3} [origin=Cartesian3.ZERO] The origin of the ray.
  610. * @param {Cartesian3} [direction=Cartesian3.ZERO] The direction of the ray.
  611. */
  612. function Ray(origin, direction) {
  613. direction = Cartographic.Cartesian3.clone(when.defaultValue(direction, Cartographic.Cartesian3.ZERO));
  614. if (!Cartographic.Cartesian3.equals(direction, Cartographic.Cartesian3.ZERO)) {
  615. Cartographic.Cartesian3.normalize(direction, direction);
  616. }
  617. /**
  618. * The origin of the ray.
  619. * @type {Cartesian3}
  620. * @default {@link Cartesian3.ZERO}
  621. */
  622. this.origin = Cartographic.Cartesian3.clone(when.defaultValue(origin, Cartographic.Cartesian3.ZERO));
  623. /**
  624. * The direction of the ray.
  625. * @type {Cartesian3}
  626. */
  627. this.direction = direction;
  628. }
  629. /**
  630. * Duplicates a Ray instance.
  631. *
  632. * @param {Ray} ray The ray to duplicate.
  633. * @param {Ray} [result] The object onto which to store the result.
  634. * @returns {Ray} The modified result parameter or a new Ray instance if one was not provided. (Returns undefined if ray is undefined)
  635. */
  636. Ray.clone = function(ray, result) {
  637. if (!when.defined(ray)) {
  638. return undefined;
  639. }
  640. if (!when.defined(result)) {
  641. return new Ray(ray.origin, ray.direction);
  642. }
  643. result.origin = Cartographic.Cartesian3.clone(ray.origin);
  644. result.direction = Cartographic.Cartesian3.clone(ray.direction);
  645. return result;
  646. };
  647. /**
  648. * Computes the point along the ray given by r(t) = o + t*d,
  649. * where o is the origin of the ray and d is the direction.
  650. *
  651. * @param {Ray} ray The ray.
  652. * @param {Number} t A scalar value.
  653. * @param {Cartesian3} [result] The object in which the result will be stored.
  654. * @returns {Cartesian3} The modified result parameter, or a new instance if none was provided.
  655. *
  656. * @example
  657. * //Get the first intersection point of a ray and an ellipsoid.
  658. * var intersection = Cesium.IntersectionTests.rayEllipsoid(ray, ellipsoid);
  659. * var point = Cesium.Ray.getPoint(ray, intersection.start);
  660. */
  661. Ray.getPoint = function(ray, t, result) {
  662. //>>includeStart('debug', pragmas.debug);
  663. Check.Check.typeOf.object('ray', ray);
  664. Check.Check.typeOf.number('t', t);
  665. //>>includeEnd('debug');
  666. if (!when.defined(result)) {
  667. result = new Cartographic.Cartesian3();
  668. }
  669. result = Cartographic.Cartesian3.multiplyByScalar(ray.direction, t, result);
  670. return Cartographic.Cartesian3.add(ray.origin, result, result);
  671. };
  672. /**
  673. * Functions for computing the intersection between geometries such as rays, planes, triangles, and ellipsoids.
  674. *
  675. * @exports IntersectionTests
  676. * @namespace
  677. */
  678. var IntersectionTests = {};
  679. /**
  680. * Computes the intersection of a ray and a plane.
  681. *
  682. * @param {Ray} ray The ray.
  683. * @param {Plane} plane The plane.
  684. * @param {Cartesian3} [result] The object onto which to store the result.
  685. * @returns {Cartesian3} The intersection point or undefined if there is no intersections.
  686. */
  687. IntersectionTests.rayPlane = function(ray, plane, result) {
  688. //>>includeStart('debug', pragmas.debug);
  689. if (!when.defined(ray)) {
  690. throw new Check.DeveloperError('ray is required.');
  691. }
  692. if (!when.defined(plane)) {
  693. throw new Check.DeveloperError('plane is required.');
  694. }
  695. //>>includeEnd('debug');
  696. if (!when.defined(result)) {
  697. result = new Cartographic.Cartesian3();
  698. }
  699. var origin = ray.origin;
  700. var direction = ray.direction;
  701. var normal = plane.normal;
  702. var denominator = Cartographic.Cartesian3.dot(normal, direction);
  703. if (Math.abs(denominator) < _Math.CesiumMath.EPSILON15) {
  704. // Ray is parallel to plane. The ray may be in the polygon's plane.
  705. return undefined;
  706. }
  707. var t = (-plane.distance - Cartographic.Cartesian3.dot(normal, origin)) / denominator;
  708. if (t < 0) {
  709. return undefined;
  710. }
  711. result = Cartographic.Cartesian3.multiplyByScalar(direction, t, result);
  712. return Cartographic.Cartesian3.add(origin, result, result);
  713. };
  714. var scratchEdge0 = new Cartographic.Cartesian3();
  715. var scratchEdge1 = new Cartographic.Cartesian3();
  716. var scratchPVec = new Cartographic.Cartesian3();
  717. var scratchTVec = new Cartographic.Cartesian3();
  718. var scratchQVec = new Cartographic.Cartesian3();
  719. /**
  720. * 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.
  721. *
  722. * Implements {@link https://cadxfem.org/inf/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf|
  723. * Fast Minimum Storage Ray/Triangle Intersection} by Tomas Moller and Ben Trumbore.
  724. *
  725. * @memberof IntersectionTests
  726. *
  727. * @param {Ray} ray The ray.
  728. * @param {Cartesian3} p0 The first vertex of the triangle.
  729. * @param {Cartesian3} p1 The second vertex of the triangle.
  730. * @param {Cartesian3} p2 The third vertex of the triangle.
  731. * @param {Boolean} [cullBackFaces=false] If <code>true</code>, will only compute an intersection with the front face of the triangle
  732. * and return undefined for intersections with the back face.
  733. * @returns {Number} The intersection as a parametric distance along the ray, or undefined if there is no intersection.
  734. */
  735. IntersectionTests.rayTriangleParametric = function(ray, p0, p1, p2, cullBackFaces) {
  736. //>>includeStart('debug', pragmas.debug);
  737. if (!when.defined(ray)) {
  738. throw new Check.DeveloperError('ray is required.');
  739. }
  740. if (!when.defined(p0)) {
  741. throw new Check.DeveloperError('p0 is required.');
  742. }
  743. if (!when.defined(p1)) {
  744. throw new Check.DeveloperError('p1 is required.');
  745. }
  746. if (!when.defined(p2)) {
  747. throw new Check.DeveloperError('p2 is required.');
  748. }
  749. //>>includeEnd('debug');
  750. cullBackFaces = when.defaultValue(cullBackFaces, false);
  751. var origin = ray.origin;
  752. var direction = ray.direction;
  753. var edge0 = Cartographic.Cartesian3.subtract(p1, p0, scratchEdge0);
  754. var edge1 = Cartographic.Cartesian3.subtract(p2, p0, scratchEdge1);
  755. var p = Cartographic.Cartesian3.cross(direction, edge1, scratchPVec);
  756. var det = Cartographic.Cartesian3.dot(edge0, p);
  757. var tvec;
  758. var q;
  759. var u;
  760. var v;
  761. var t;
  762. if (cullBackFaces) {
  763. if (det < _Math.CesiumMath.EPSILON6) {
  764. return undefined;
  765. }
  766. tvec = Cartographic.Cartesian3.subtract(origin, p0, scratchTVec);
  767. u = Cartographic.Cartesian3.dot(tvec, p);
  768. if (u < 0.0 || u > det) {
  769. return undefined;
  770. }
  771. q = Cartographic.Cartesian3.cross(tvec, edge0, scratchQVec);
  772. v = Cartographic.Cartesian3.dot(direction, q);
  773. if (v < 0.0 || u + v > det) {
  774. return undefined;
  775. }
  776. t = Cartographic.Cartesian3.dot(edge1, q) / det;
  777. } else {
  778. if (Math.abs(det) < _Math.CesiumMath.EPSILON6) {
  779. return undefined;
  780. }
  781. var invDet = 1.0 / det;
  782. tvec = Cartographic.Cartesian3.subtract(origin, p0, scratchTVec);
  783. u = Cartographic.Cartesian3.dot(tvec, p) * invDet;
  784. if (u < 0.0 || u > 1.0) {
  785. return undefined;
  786. }
  787. q = Cartographic.Cartesian3.cross(tvec, edge0, scratchQVec);
  788. v = Cartographic.Cartesian3.dot(direction, q) * invDet;
  789. if (v < 0.0 || u + v > 1.0) {
  790. return undefined;
  791. }
  792. t = Cartographic.Cartesian3.dot(edge1, q) * invDet;
  793. }
  794. return t;
  795. };
  796. /**
  797. * Computes the intersection of a ray and a triangle as a Cartesian3 coordinate.
  798. *
  799. * Implements {@link https://cadxfem.org/inf/Fast%20MinimumStorage%20RayTriangle%20Intersection.pdf|
  800. * Fast Minimum Storage Ray/Triangle Intersection} by Tomas Moller and Ben Trumbore.
  801. *
  802. * @memberof IntersectionTests
  803. *
  804. * @param {Ray} ray The ray.
  805. * @param {Cartesian3} p0 The first vertex of the triangle.
  806. * @param {Cartesian3} p1 The second vertex of the triangle.
  807. * @param {Cartesian3} p2 The third vertex of the triangle.
  808. * @param {Boolean} [cullBackFaces=false] If <code>true</code>, will only compute an intersection with the front face of the triangle
  809. * and return undefined for intersections with the back face.
  810. * @param {Cartesian3} [result] The <code>Cartesian3</code> onto which to store the result.
  811. * @returns {Cartesian3} The intersection point or undefined if there is no intersections.
  812. */
  813. IntersectionTests.rayTriangle = function(ray, p0, p1, p2, cullBackFaces, result) {
  814. var t = IntersectionTests.rayTriangleParametric(ray, p0, p1, p2, cullBackFaces);
  815. if (!when.defined(t) || t < 0.0) {
  816. return undefined;
  817. }
  818. if (!when.defined(result)) {
  819. result = new Cartographic.Cartesian3();
  820. }
  821. Cartographic.Cartesian3.multiplyByScalar(ray.direction, t, result);
  822. return Cartographic.Cartesian3.add(ray.origin, result, result);
  823. };
  824. var scratchLineSegmentTriangleRay = new Ray();
  825. /**
  826. * Computes the intersection of a line segment and a triangle.
  827. * @memberof IntersectionTests
  828. *
  829. * @param {Cartesian3} v0 The an end point of the line segment.
  830. * @param {Cartesian3} v1 The other end point of the line segment.
  831. * @param {Cartesian3} p0 The first vertex of the triangle.
  832. * @param {Cartesian3} p1 The second vertex of the triangle.
  833. * @param {Cartesian3} p2 The third vertex of the triangle.
  834. * @param {Boolean} [cullBackFaces=false] If <code>true</code>, will only compute an intersection with the front face of the triangle
  835. * and return undefined for intersections with the back face.
  836. * @param {Cartesian3} [result] The <code>Cartesian3</code> onto which to store the result.
  837. * @returns {Cartesian3} The intersection point or undefined if there is no intersections.
  838. */
  839. IntersectionTests.lineSegmentTriangle = function(v0, v1, p0, p1, p2, cullBackFaces, result) {
  840. //>>includeStart('debug', pragmas.debug);
  841. if (!when.defined(v0)) {
  842. throw new Check.DeveloperError('v0 is required.');
  843. }
  844. if (!when.defined(v1)) {
  845. throw new Check.DeveloperError('v1 is required.');
  846. }
  847. if (!when.defined(p0)) {
  848. throw new Check.DeveloperError('p0 is required.');
  849. }
  850. if (!when.defined(p1)) {
  851. throw new Check.DeveloperError('p1 is required.');
  852. }
  853. if (!when.defined(p2)) {
  854. throw new Check.DeveloperError('p2 is required.');
  855. }
  856. //>>includeEnd('debug');
  857. var ray = scratchLineSegmentTriangleRay;
  858. Cartographic.Cartesian3.clone(v0, ray.origin);
  859. Cartographic.Cartesian3.subtract(v1, v0, ray.direction);
  860. Cartographic.Cartesian3.normalize(ray.direction, ray.direction);
  861. var t = IntersectionTests.rayTriangleParametric(ray, p0, p1, p2, cullBackFaces);
  862. if (!when.defined(t) || t < 0.0 || t > Cartographic.Cartesian3.distance(v0, v1)) {
  863. return undefined;
  864. }
  865. if (!when.defined(result)) {
  866. result = new Cartographic.Cartesian3();
  867. }
  868. Cartographic.Cartesian3.multiplyByScalar(ray.direction, t, result);
  869. return Cartographic.Cartesian3.add(ray.origin, result, result);
  870. };
  871. function solveQuadratic(a, b, c, result) {
  872. var det = b * b - 4.0 * a * c;
  873. if (det < 0.0) {
  874. return undefined;
  875. } else if (det > 0.0) {
  876. var denom = 1.0 / (2.0 * a);
  877. var disc = Math.sqrt(det);
  878. var root0 = (-b + disc) * denom;
  879. var root1 = (-b - disc) * denom;
  880. if (root0 < root1) {
  881. result.root0 = root0;
  882. result.root1 = root1;
  883. } else {
  884. result.root0 = root1;
  885. result.root1 = root0;
  886. }
  887. return result;
  888. }
  889. var root = -b / (2.0 * a);
  890. if (root === 0.0) {
  891. return undefined;
  892. }
  893. result.root0 = result.root1 = root;
  894. return result;
  895. }
  896. var raySphereRoots = {
  897. root0 : 0.0,
  898. root1 : 0.0
  899. };
  900. function raySphere(ray, sphere, result) {
  901. if (!when.defined(result)) {
  902. result = new BoundingSphere.Interval();
  903. }
  904. var origin = ray.origin;
  905. var direction = ray.direction;
  906. var center = sphere.center;
  907. var radiusSquared = sphere.radius * sphere.radius;
  908. var diff = Cartographic.Cartesian3.subtract(origin, center, scratchPVec);
  909. var a = Cartographic.Cartesian3.dot(direction, direction);
  910. var b = 2.0 * Cartographic.Cartesian3.dot(direction, diff);
  911. var c = Cartographic.Cartesian3.magnitudeSquared(diff) - radiusSquared;
  912. var roots = solveQuadratic(a, b, c, raySphereRoots);
  913. if (!when.defined(roots)) {
  914. return undefined;
  915. }
  916. result.start = roots.root0;
  917. result.stop = roots.root1;
  918. return result;
  919. }
  920. /**
  921. * Computes the intersection points of a ray with a sphere.
  922. * @memberof IntersectionTests
  923. *
  924. * @param {Ray} ray The ray.
  925. * @param {BoundingSphere} sphere The sphere.
  926. * @param {Interval} [result] The result onto which to store the result.
  927. * @returns {Interval} The interval containing scalar points along the ray or undefined if there are no intersections.
  928. */
  929. IntersectionTests.raySphere = function(ray, sphere, result) {
  930. //>>includeStart('debug', pragmas.debug);
  931. if (!when.defined(ray)) {
  932. throw new Check.DeveloperError('ray is required.');
  933. }
  934. if (!when.defined(sphere)) {
  935. throw new Check.DeveloperError('sphere is required.');
  936. }
  937. //>>includeEnd('debug');
  938. result = raySphere(ray, sphere, result);
  939. if (!when.defined(result) || result.stop < 0.0) {
  940. return undefined;
  941. }
  942. result.start = Math.max(result.start, 0.0);
  943. return result;
  944. };
  945. var scratchLineSegmentRay = new Ray();
  946. /**
  947. * Computes the intersection points of a line segment with a sphere.
  948. * @memberof IntersectionTests
  949. *
  950. * @param {Cartesian3} p0 An end point of the line segment.
  951. * @param {Cartesian3} p1 The other end point of the line segment.
  952. * @param {BoundingSphere} sphere The sphere.
  953. * @param {Interval} [result] The result onto which to store the result.
  954. * @returns {Interval} The interval containing scalar points along the ray or undefined if there are no intersections.
  955. */
  956. IntersectionTests.lineSegmentSphere = function(p0, p1, sphere, result) {
  957. //>>includeStart('debug', pragmas.debug);
  958. if (!when.defined(p0)) {
  959. throw new Check.DeveloperError('p0 is required.');
  960. }
  961. if (!when.defined(p1)) {
  962. throw new Check.DeveloperError('p1 is required.');
  963. }
  964. if (!when.defined(sphere)) {
  965. throw new Check.DeveloperError('sphere is required.');
  966. }
  967. //>>includeEnd('debug');
  968. var ray = scratchLineSegmentRay;
  969. Cartographic.Cartesian3.clone(p0, ray.origin);
  970. var direction = Cartographic.Cartesian3.subtract(p1, p0, ray.direction);
  971. var maxT = Cartographic.Cartesian3.magnitude(direction);
  972. Cartographic.Cartesian3.normalize(direction, direction);
  973. result = raySphere(ray, sphere, result);
  974. if (!when.defined(result) || result.stop < 0.0 || result.start > maxT) {
  975. return undefined;
  976. }
  977. result.start = Math.max(result.start, 0.0);
  978. result.stop = Math.min(result.stop, maxT);
  979. return result;
  980. };
  981. var scratchQ = new Cartographic.Cartesian3();
  982. var scratchW = new Cartographic.Cartesian3();
  983. /**
  984. * Computes the intersection points of a ray with an ellipsoid.
  985. *
  986. * @param {Ray} ray The ray.
  987. * @param {Ellipsoid} ellipsoid The ellipsoid.
  988. * @returns {Interval} The interval containing scalar points along the ray or undefined if there are no intersections.
  989. */
  990. IntersectionTests.rayEllipsoid = function(ray, ellipsoid) {
  991. //>>includeStart('debug', pragmas.debug);
  992. if (!when.defined(ray)) {
  993. throw new Check.DeveloperError('ray is required.');
  994. }
  995. if (!when.defined(ellipsoid)) {
  996. throw new Check.DeveloperError('ellipsoid is required.');
  997. }
  998. //>>includeEnd('debug');
  999. var inverseRadii = ellipsoid.oneOverRadii;
  1000. var q = Cartographic.Cartesian3.multiplyComponents(inverseRadii, ray.origin, scratchQ);
  1001. var w = Cartographic.Cartesian3.multiplyComponents(inverseRadii, ray.direction, scratchW);
  1002. var q2 = Cartographic.Cartesian3.magnitudeSquared(q);
  1003. var qw = Cartographic.Cartesian3.dot(q, w);
  1004. var difference, w2, product, discriminant, temp;
  1005. if (q2 > 1.0) {
  1006. // Outside ellipsoid.
  1007. if (qw >= 0.0) {
  1008. // Looking outward or tangent (0 intersections).
  1009. return undefined;
  1010. }
  1011. // qw < 0.0.
  1012. var qw2 = qw * qw;
  1013. difference = q2 - 1.0; // Positively valued.
  1014. w2 = Cartographic.Cartesian3.magnitudeSquared(w);
  1015. product = w2 * difference;
  1016. if (qw2 < product) {
  1017. // Imaginary roots (0 intersections).
  1018. return undefined;
  1019. } else if (qw2 > product) {
  1020. // Distinct roots (2 intersections).
  1021. discriminant = qw * qw - product;
  1022. temp = -qw + Math.sqrt(discriminant); // Avoid cancellation.
  1023. var root0 = temp / w2;
  1024. var root1 = difference / temp;
  1025. if (root0 < root1) {
  1026. return new BoundingSphere.Interval(root0, root1);
  1027. }
  1028. return {
  1029. start : root1,
  1030. stop : root0
  1031. };
  1032. }
  1033. // qw2 == product. Repeated roots (2 intersections).
  1034. var root = Math.sqrt(difference / w2);
  1035. return new BoundingSphere.Interval(root, root);
  1036. } else if (q2 < 1.0) {
  1037. // Inside ellipsoid (2 intersections).
  1038. difference = q2 - 1.0; // Negatively valued.
  1039. w2 = Cartographic.Cartesian3.magnitudeSquared(w);
  1040. product = w2 * difference; // Negatively valued.
  1041. discriminant = qw * qw - product;
  1042. temp = -qw + Math.sqrt(discriminant); // Positively valued.
  1043. return new BoundingSphere.Interval(0.0, temp / w2);
  1044. }
  1045. // q2 == 1.0. On ellipsoid.
  1046. if (qw < 0.0) {
  1047. // Looking inward.
  1048. w2 = Cartographic.Cartesian3.magnitudeSquared(w);
  1049. return new BoundingSphere.Interval(0.0, -qw / w2);
  1050. }
  1051. // qw >= 0.0. Looking outward or tangent.
  1052. return undefined;
  1053. };
  1054. function addWithCancellationCheck$1(left, right, tolerance) {
  1055. var difference = left + right;
  1056. if ((_Math.CesiumMath.sign(left) !== _Math.CesiumMath.sign(right)) &&
  1057. Math.abs(difference / Math.max(Math.abs(left), Math.abs(right))) < tolerance) {
  1058. return 0.0;
  1059. }
  1060. return difference;
  1061. }
  1062. function quadraticVectorExpression(A, b, c, x, w) {
  1063. var xSquared = x * x;
  1064. var wSquared = w * w;
  1065. var l2 = (A[BoundingSphere.Matrix3.COLUMN1ROW1] - A[BoundingSphere.Matrix3.COLUMN2ROW2]) * wSquared;
  1066. var l1 = w * (x * addWithCancellationCheck$1(A[BoundingSphere.Matrix3.COLUMN1ROW0], A[BoundingSphere.Matrix3.COLUMN0ROW1], _Math.CesiumMath.EPSILON15) + b.y);
  1067. var l0 = (A[BoundingSphere.Matrix3.COLUMN0ROW0] * xSquared + A[BoundingSphere.Matrix3.COLUMN2ROW2] * wSquared) + x * b.x + c;
  1068. var r1 = wSquared * addWithCancellationCheck$1(A[BoundingSphere.Matrix3.COLUMN2ROW1], A[BoundingSphere.Matrix3.COLUMN1ROW2], _Math.CesiumMath.EPSILON15);
  1069. var r0 = w * (x * addWithCancellationCheck$1(A[BoundingSphere.Matrix3.COLUMN2ROW0], A[BoundingSphere.Matrix3.COLUMN0ROW2]) + b.z);
  1070. var cosines;
  1071. var solutions = [];
  1072. if (r0 === 0.0 && r1 === 0.0) {
  1073. cosines = QuadraticRealPolynomial.computeRealRoots(l2, l1, l0);
  1074. if (cosines.length === 0) {
  1075. return solutions;
  1076. }
  1077. var cosine0 = cosines[0];
  1078. var sine0 = Math.sqrt(Math.max(1.0 - cosine0 * cosine0, 0.0));
  1079. solutions.push(new Cartographic.Cartesian3(x, w * cosine0, w * -sine0));
  1080. solutions.push(new Cartographic.Cartesian3(x, w * cosine0, w * sine0));
  1081. if (cosines.length === 2) {
  1082. var cosine1 = cosines[1];
  1083. var sine1 = Math.sqrt(Math.max(1.0 - cosine1 * cosine1, 0.0));
  1084. solutions.push(new Cartographic.Cartesian3(x, w * cosine1, w * -sine1));
  1085. solutions.push(new Cartographic.Cartesian3(x, w * cosine1, w * sine1));
  1086. }
  1087. return solutions;
  1088. }
  1089. var r0Squared = r0 * r0;
  1090. var r1Squared = r1 * r1;
  1091. var l2Squared = l2 * l2;
  1092. var r0r1 = r0 * r1;
  1093. var c4 = l2Squared + r1Squared;
  1094. var c3 = 2.0 * (l1 * l2 + r0r1);
  1095. var c2 = 2.0 * l0 * l2 + l1 * l1 - r1Squared + r0Squared;
  1096. var c1 = 2.0 * (l0 * l1 - r0r1);
  1097. var c0 = l0 * l0 - r0Squared;
  1098. if (c4 === 0.0 && c3 === 0.0 && c2 === 0.0 && c1 === 0.0) {
  1099. return solutions;
  1100. }
  1101. cosines = QuarticRealPolynomial.computeRealRoots(c4, c3, c2, c1, c0);
  1102. var length = cosines.length;
  1103. if (length === 0) {
  1104. return solutions;
  1105. }
  1106. for ( var i = 0; i < length; ++i) {
  1107. var cosine = cosines[i];
  1108. var cosineSquared = cosine * cosine;
  1109. var sineSquared = Math.max(1.0 - cosineSquared, 0.0);
  1110. var sine = Math.sqrt(sineSquared);
  1111. //var left = l2 * cosineSquared + l1 * cosine + l0;
  1112. var left;
  1113. if (_Math.CesiumMath.sign(l2) === _Math.CesiumMath.sign(l0)) {
  1114. left = addWithCancellationCheck$1(l2 * cosineSquared + l0, l1 * cosine, _Math.CesiumMath.EPSILON12);
  1115. } else if (_Math.CesiumMath.sign(l0) === _Math.CesiumMath.sign(l1 * cosine)) {
  1116. left = addWithCancellationCheck$1(l2 * cosineSquared, l1 * cosine + l0, _Math.CesiumMath.EPSILON12);
  1117. } else {
  1118. left = addWithCancellationCheck$1(l2 * cosineSquared + l1 * cosine, l0, _Math.CesiumMath.EPSILON12);
  1119. }
  1120. var right = addWithCancellationCheck$1(r1 * cosine, r0, _Math.CesiumMath.EPSILON15);
  1121. var product = left * right;
  1122. if (product < 0.0) {
  1123. solutions.push(new Cartographic.Cartesian3(x, w * cosine, w * sine));
  1124. } else if (product > 0.0) {
  1125. solutions.push(new Cartographic.Cartesian3(x, w * cosine, w * -sine));
  1126. } else if (sine !== 0.0) {
  1127. solutions.push(new Cartographic.Cartesian3(x, w * cosine, w * -sine));
  1128. solutions.push(new Cartographic.Cartesian3(x, w * cosine, w * sine));
  1129. ++i;
  1130. } else {
  1131. solutions.push(new Cartographic.Cartesian3(x, w * cosine, w * sine));
  1132. }
  1133. }
  1134. return solutions;
  1135. }
  1136. var firstAxisScratch = new Cartographic.Cartesian3();
  1137. var secondAxisScratch = new Cartographic.Cartesian3();
  1138. var thirdAxisScratch = new Cartographic.Cartesian3();
  1139. var referenceScratch = new Cartographic.Cartesian3();
  1140. var bCart = new Cartographic.Cartesian3();
  1141. var bScratch = new BoundingSphere.Matrix3();
  1142. var btScratch = new BoundingSphere.Matrix3();
  1143. var diScratch = new BoundingSphere.Matrix3();
  1144. var dScratch = new BoundingSphere.Matrix3();
  1145. var cScratch = new BoundingSphere.Matrix3();
  1146. var tempMatrix = new BoundingSphere.Matrix3();
  1147. var aScratch = new BoundingSphere.Matrix3();
  1148. var sScratch = new Cartographic.Cartesian3();
  1149. var closestScratch = new Cartographic.Cartesian3();
  1150. var surfPointScratch = new Cartographic.Cartographic();
  1151. /**
  1152. * Provides the point along the ray which is nearest to the ellipsoid.
  1153. *
  1154. * @param {Ray} ray The ray.
  1155. * @param {Ellipsoid} ellipsoid The ellipsoid.
  1156. * @returns {Cartesian3} The nearest planetodetic point on the ray.
  1157. */
  1158. IntersectionTests.grazingAltitudeLocation = function(ray, ellipsoid) {
  1159. //>>includeStart('debug', pragmas.debug);
  1160. if (!when.defined(ray)) {
  1161. throw new Check.DeveloperError('ray is required.');
  1162. }
  1163. if (!when.defined(ellipsoid)) {
  1164. throw new Check.DeveloperError('ellipsoid is required.');
  1165. }
  1166. //>>includeEnd('debug');
  1167. var position = ray.origin;
  1168. var direction = ray.direction;
  1169. if (!Cartographic.Cartesian3.equals(position, Cartographic.Cartesian3.ZERO)) {
  1170. var normal = ellipsoid.geodeticSurfaceNormal(position, firstAxisScratch);
  1171. if (Cartographic.Cartesian3.dot(direction, normal) >= 0.0) { // The location provided is the closest point in altitude
  1172. return position;
  1173. }
  1174. }
  1175. var intersects = when.defined(this.rayEllipsoid(ray, ellipsoid));
  1176. // Compute the scaled direction vector.
  1177. var f = ellipsoid.transformPositionToScaledSpace(direction, firstAxisScratch);
  1178. // Constructs a basis from the unit scaled direction vector. Construct its rotation and transpose.
  1179. var firstAxis = Cartographic.Cartesian3.normalize(f, f);
  1180. var reference = Cartographic.Cartesian3.mostOrthogonalAxis(f, referenceScratch);
  1181. var secondAxis = Cartographic.Cartesian3.normalize(Cartographic.Cartesian3.cross(reference, firstAxis, secondAxisScratch), secondAxisScratch);
  1182. var thirdAxis = Cartographic.Cartesian3.normalize(Cartographic.Cartesian3.cross(firstAxis, secondAxis, thirdAxisScratch), thirdAxisScratch);
  1183. var B = bScratch;
  1184. B[0] = firstAxis.x;
  1185. B[1] = firstAxis.y;
  1186. B[2] = firstAxis.z;
  1187. B[3] = secondAxis.x;
  1188. B[4] = secondAxis.y;
  1189. B[5] = secondAxis.z;
  1190. B[6] = thirdAxis.x;
  1191. B[7] = thirdAxis.y;
  1192. B[8] = thirdAxis.z;
  1193. var B_T = BoundingSphere.Matrix3.transpose(B, btScratch);
  1194. // Get the scaling matrix and its inverse.
  1195. var D_I = BoundingSphere.Matrix3.fromScale(ellipsoid.radii, diScratch);
  1196. var D = BoundingSphere.Matrix3.fromScale(ellipsoid.oneOverRadii, dScratch);
  1197. var C = cScratch;
  1198. C[0] = 0.0;
  1199. C[1] = -direction.z;
  1200. C[2] = direction.y;
  1201. C[3] = direction.z;
  1202. C[4] = 0.0;
  1203. C[5] = -direction.x;
  1204. C[6] = -direction.y;
  1205. C[7] = direction.x;
  1206. C[8] = 0.0;
  1207. var temp = BoundingSphere.Matrix3.multiply(BoundingSphere.Matrix3.multiply(B_T, D, tempMatrix), C, tempMatrix);
  1208. var A = BoundingSphere.Matrix3.multiply(BoundingSphere.Matrix3.multiply(temp, D_I, aScratch), B, aScratch);
  1209. var b = BoundingSphere.Matrix3.multiplyByVector(temp, position, bCart);
  1210. // Solve for the solutions to the expression in standard form:
  1211. var solutions = quadraticVectorExpression(A, Cartographic.Cartesian3.negate(b, firstAxisScratch), 0.0, 0.0, 1.0);
  1212. var s;
  1213. var altitude;
  1214. var length = solutions.length;
  1215. if (length > 0) {
  1216. var closest = Cartographic.Cartesian3.clone(Cartographic.Cartesian3.ZERO, closestScratch);
  1217. var maximumValue = Number.NEGATIVE_INFINITY;
  1218. for ( var i = 0; i < length; ++i) {
  1219. s = BoundingSphere.Matrix3.multiplyByVector(D_I, BoundingSphere.Matrix3.multiplyByVector(B, solutions[i], sScratch), sScratch);
  1220. var v = Cartographic.Cartesian3.normalize(Cartographic.Cartesian3.subtract(s, position, referenceScratch), referenceScratch);
  1221. var dotProduct = Cartographic.Cartesian3.dot(v, direction);
  1222. if (dotProduct > maximumValue) {
  1223. maximumValue = dotProduct;
  1224. closest = Cartographic.Cartesian3.clone(s, closest);
  1225. }
  1226. }
  1227. var surfacePoint = ellipsoid.cartesianToCartographic(closest, surfPointScratch);
  1228. maximumValue = _Math.CesiumMath.clamp(maximumValue, 0.0, 1.0);
  1229. altitude = Cartographic.Cartesian3.magnitude(Cartographic.Cartesian3.subtract(closest, position, referenceScratch)) * Math.sqrt(1.0 - maximumValue * maximumValue);
  1230. altitude = intersects ? -altitude : altitude;
  1231. surfacePoint.height = altitude;
  1232. return ellipsoid.cartographicToCartesian(surfacePoint, new Cartographic.Cartesian3());
  1233. }
  1234. return undefined;
  1235. };
  1236. var lineSegmentPlaneDifference = new Cartographic.Cartesian3();
  1237. /**
  1238. * Computes the intersection of a line segment and a plane.
  1239. *
  1240. * @param {Cartesian3} endPoint0 An end point of the line segment.
  1241. * @param {Cartesian3} endPoint1 The other end point of the line segment.
  1242. * @param {Plane} plane The plane.
  1243. * @param {Cartesian3} [result] The object onto which to store the result.
  1244. * @returns {Cartesian3} The intersection point or undefined if there is no intersection.
  1245. *
  1246. * @example
  1247. * var origin = Cesium.Cartesian3.fromDegrees(-75.59777, 40.03883);
  1248. * var normal = ellipsoid.geodeticSurfaceNormal(origin);
  1249. * var plane = Cesium.Plane.fromPointNormal(origin, normal);
  1250. *
  1251. * var p0 = new Cesium.Cartesian3(...);
  1252. * var p1 = new Cesium.Cartesian3(...);
  1253. *
  1254. * // find the intersection of the line segment from p0 to p1 and the tangent plane at origin.
  1255. * var intersection = Cesium.IntersectionTests.lineSegmentPlane(p0, p1, plane);
  1256. */
  1257. IntersectionTests.lineSegmentPlane = function(endPoint0, endPoint1, plane, result) {
  1258. //>>includeStart('debug', pragmas.debug);
  1259. if (!when.defined(endPoint0)) {
  1260. throw new Check.DeveloperError('endPoint0 is required.');
  1261. }
  1262. if (!when.defined(endPoint1)) {
  1263. throw new Check.DeveloperError('endPoint1 is required.');
  1264. }
  1265. if (!when.defined(plane)) {
  1266. throw new Check.DeveloperError('plane is required.');
  1267. }
  1268. //>>includeEnd('debug');
  1269. if (!when.defined(result)) {
  1270. result = new Cartographic.Cartesian3();
  1271. }
  1272. var difference = Cartographic.Cartesian3.subtract(endPoint1, endPoint0, lineSegmentPlaneDifference);
  1273. var normal = plane.normal;
  1274. var nDotDiff = Cartographic.Cartesian3.dot(normal, difference);
  1275. // check if the segment and plane are parallel
  1276. if (Math.abs(nDotDiff) < _Math.CesiumMath.EPSILON6) {
  1277. return undefined;
  1278. }
  1279. var nDotP0 = Cartographic.Cartesian3.dot(normal, endPoint0);
  1280. var t = -(plane.distance + nDotP0) / nDotDiff;
  1281. // intersection only if t is in [0, 1]
  1282. if (t < 0.0 || t > 1.0) {
  1283. return undefined;
  1284. }
  1285. // intersection is endPoint0 + t * (endPoint1 - endPoint0)
  1286. Cartographic.Cartesian3.multiplyByScalar(difference, t, result);
  1287. Cartographic.Cartesian3.add(endPoint0, result, result);
  1288. return result;
  1289. };
  1290. /**
  1291. * Computes the intersection of a triangle and a plane
  1292. *
  1293. * @param {Cartesian3} p0 First point of the triangle
  1294. * @param {Cartesian3} p1 Second point of the triangle
  1295. * @param {Cartesian3} p2 Third point of the triangle
  1296. * @param {Plane} plane Intersection plane
  1297. * @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)
  1298. *
  1299. * @example
  1300. * var origin = Cesium.Cartesian3.fromDegrees(-75.59777, 40.03883);
  1301. * var normal = ellipsoid.geodeticSurfaceNormal(origin);
  1302. * var plane = Cesium.Plane.fromPointNormal(origin, normal);
  1303. *
  1304. * var p0 = new Cesium.Cartesian3(...);
  1305. * var p1 = new Cesium.Cartesian3(...);
  1306. * var p2 = new Cesium.Cartesian3(...);
  1307. *
  1308. * // convert the triangle composed of points (p0, p1, p2) to three triangles that don't cross the plane
  1309. * var triangles = Cesium.IntersectionTests.trianglePlaneIntersection(p0, p1, p2, plane);
  1310. */
  1311. IntersectionTests.trianglePlaneIntersection = function(p0, p1, p2, plane) {
  1312. //>>includeStart('debug', pragmas.debug);
  1313. if ((!when.defined(p0)) ||
  1314. (!when.defined(p1)) ||
  1315. (!when.defined(p2)) ||
  1316. (!when.defined(plane))) {
  1317. throw new Check.DeveloperError('p0, p1, p2, and plane are required.');
  1318. }
  1319. //>>includeEnd('debug');
  1320. var planeNormal = plane.normal;
  1321. var planeD = plane.distance;
  1322. var p0Behind = (Cartographic.Cartesian3.dot(planeNormal, p0) + planeD) < 0.0;
  1323. var p1Behind = (Cartographic.Cartesian3.dot(planeNormal, p1) + planeD) < 0.0;
  1324. var p2Behind = (Cartographic.Cartesian3.dot(planeNormal, p2) + planeD) < 0.0;
  1325. // Given these dots products, the calls to lineSegmentPlaneIntersection
  1326. // always have defined results.
  1327. var numBehind = 0;
  1328. numBehind += p0Behind ? 1 : 0;
  1329. numBehind += p1Behind ? 1 : 0;
  1330. numBehind += p2Behind ? 1 : 0;
  1331. var u1, u2;
  1332. if (numBehind === 1 || numBehind === 2) {
  1333. u1 = new Cartographic.Cartesian3();
  1334. u2 = new Cartographic.Cartesian3();
  1335. }
  1336. if (numBehind === 1) {
  1337. if (p0Behind) {
  1338. IntersectionTests.lineSegmentPlane(p0, p1, plane, u1);
  1339. IntersectionTests.lineSegmentPlane(p0, p2, plane, u2);
  1340. return {
  1341. positions : [p0, p1, p2, u1, u2 ],
  1342. indices : [
  1343. // Behind
  1344. 0, 3, 4,
  1345. // In front
  1346. 1, 2, 4,
  1347. 1, 4, 3
  1348. ]
  1349. };
  1350. } else if (p1Behind) {
  1351. IntersectionTests.lineSegmentPlane(p1, p2, plane, u1);
  1352. IntersectionTests.lineSegmentPlane(p1, p0, plane, u2);
  1353. return {
  1354. positions : [p0, p1, p2, u1, u2 ],
  1355. indices : [
  1356. // Behind
  1357. 1, 3, 4,
  1358. // In front
  1359. 2, 0, 4,
  1360. 2, 4, 3
  1361. ]
  1362. };
  1363. } else if (p2Behind) {
  1364. IntersectionTests.lineSegmentPlane(p2, p0, plane, u1);
  1365. IntersectionTests.lineSegmentPlane(p2, p1, plane, u2);
  1366. return {
  1367. positions : [p0, p1, p2, u1, u2 ],
  1368. indices : [
  1369. // Behind
  1370. 2, 3, 4,
  1371. // In front
  1372. 0, 1, 4,
  1373. 0, 4, 3
  1374. ]
  1375. };
  1376. }
  1377. } else if (numBehind === 2) {
  1378. if (!p0Behind) {
  1379. IntersectionTests.lineSegmentPlane(p1, p0, plane, u1);
  1380. IntersectionTests.lineSegmentPlane(p2, p0, plane, u2);
  1381. return {
  1382. positions : [p0, p1, p2, u1, u2 ],
  1383. indices : [
  1384. // Behind
  1385. 1, 2, 4,
  1386. 1, 4, 3,
  1387. // In front
  1388. 0, 3, 4
  1389. ]
  1390. };
  1391. } else if (!p1Behind) {
  1392. IntersectionTests.lineSegmentPlane(p2, p1, plane, u1);
  1393. IntersectionTests.lineSegmentPlane(p0, p1, plane, u2);
  1394. return {
  1395. positions : [p0, p1, p2, u1, u2 ],
  1396. indices : [
  1397. // Behind
  1398. 2, 0, 4,
  1399. 2, 4, 3,
  1400. // In front
  1401. 1, 3, 4
  1402. ]
  1403. };
  1404. } else if (!p2Behind) {
  1405. IntersectionTests.lineSegmentPlane(p0, p2, plane, u1);
  1406. IntersectionTests.lineSegmentPlane(p1, p2, plane, u2);
  1407. return {
  1408. positions : [p0, p1, p2, u1, u2 ],
  1409. indices : [
  1410. // Behind
  1411. 0, 1, 4,
  1412. 0, 4, 3,
  1413. // In front
  1414. 2, 3, 4
  1415. ]
  1416. };
  1417. }
  1418. }
  1419. // if numBehind is 3, the triangle is completely behind the plane;
  1420. // otherwise, it is completely in front (numBehind is 0).
  1421. return undefined;
  1422. };
  1423. exports.IntersectionTests = IntersectionTests;
  1424. exports.Ray = Ray;
  1425. });