IntersectionTests-c2360ffa.js 66 KB

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