EllipseGeometryLibrary.js 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365
  1. import Cartesian3 from "./Cartesian3.js";
  2. import CesiumMath from "./Math.js";
  3. import Matrix3 from "./Matrix3.js";
  4. import Quaternion from "./Quaternion.js";
  5. const EllipseGeometryLibrary = {};
  6. const rotAxis = new Cartesian3();
  7. const tempVec = new Cartesian3();
  8. const unitQuat = new Quaternion();
  9. const rotMtx = new Matrix3();
  10. function pointOnEllipsoid(
  11. theta,
  12. rotation,
  13. northVec,
  14. eastVec,
  15. aSqr,
  16. ab,
  17. bSqr,
  18. mag,
  19. unitPos,
  20. result
  21. ) {
  22. const azimuth = theta + rotation;
  23. Cartesian3.multiplyByScalar(eastVec, Math.cos(azimuth), rotAxis);
  24. Cartesian3.multiplyByScalar(northVec, Math.sin(azimuth), tempVec);
  25. Cartesian3.add(rotAxis, tempVec, rotAxis);
  26. let cosThetaSquared = Math.cos(theta);
  27. cosThetaSquared = cosThetaSquared * cosThetaSquared;
  28. let sinThetaSquared = Math.sin(theta);
  29. sinThetaSquared = sinThetaSquared * sinThetaSquared;
  30. const radius =
  31. ab / Math.sqrt(bSqr * cosThetaSquared + aSqr * sinThetaSquared);
  32. const angle = radius / mag;
  33. // Create the quaternion to rotate the position vector to the boundary of the ellipse.
  34. Quaternion.fromAxisAngle(rotAxis, angle, unitQuat);
  35. Matrix3.fromQuaternion(unitQuat, rotMtx);
  36. Matrix3.multiplyByVector(rotMtx, unitPos, result);
  37. Cartesian3.normalize(result, result);
  38. Cartesian3.multiplyByScalar(result, mag, result);
  39. return result;
  40. }
  41. const scratchCartesian1 = new Cartesian3();
  42. const scratchCartesian2 = new Cartesian3();
  43. const scratchCartesian3 = new Cartesian3();
  44. const scratchNormal = new Cartesian3();
  45. /**
  46. * Returns the positions raised to the given heights
  47. * @private
  48. */
  49. EllipseGeometryLibrary.raisePositionsToHeight = function (
  50. positions,
  51. options,
  52. extrude
  53. ) {
  54. const ellipsoid = options.ellipsoid;
  55. const height = options.height;
  56. const extrudedHeight = options.extrudedHeight;
  57. const size = extrude ? (positions.length / 3) * 2 : positions.length / 3;
  58. const finalPositions = new Float64Array(size * 3);
  59. const length = positions.length;
  60. const bottomOffset = extrude ? length : 0;
  61. for (let i = 0; i < length; i += 3) {
  62. const i1 = i + 1;
  63. const i2 = i + 2;
  64. const position = Cartesian3.fromArray(positions, i, scratchCartesian1);
  65. ellipsoid.scaleToGeodeticSurface(position, position);
  66. const extrudedPosition = Cartesian3.clone(position, scratchCartesian2);
  67. const normal = ellipsoid.geodeticSurfaceNormal(position, scratchNormal);
  68. const scaledNormal = Cartesian3.multiplyByScalar(
  69. normal,
  70. height,
  71. scratchCartesian3
  72. );
  73. Cartesian3.add(position, scaledNormal, position);
  74. if (extrude) {
  75. Cartesian3.multiplyByScalar(normal, extrudedHeight, scaledNormal);
  76. Cartesian3.add(extrudedPosition, scaledNormal, extrudedPosition);
  77. finalPositions[i + bottomOffset] = extrudedPosition.x;
  78. finalPositions[i1 + bottomOffset] = extrudedPosition.y;
  79. finalPositions[i2 + bottomOffset] = extrudedPosition.z;
  80. }
  81. finalPositions[i] = position.x;
  82. finalPositions[i1] = position.y;
  83. finalPositions[i2] = position.z;
  84. }
  85. return finalPositions;
  86. };
  87. const unitPosScratch = new Cartesian3();
  88. const eastVecScratch = new Cartesian3();
  89. const northVecScratch = new Cartesian3();
  90. /**
  91. * Returns an array of positions that make up the ellipse.
  92. * @private
  93. */
  94. EllipseGeometryLibrary.computeEllipsePositions = function (
  95. options,
  96. addFillPositions,
  97. addEdgePositions
  98. ) {
  99. const semiMinorAxis = options.semiMinorAxis;
  100. const semiMajorAxis = options.semiMajorAxis;
  101. const rotation = options.rotation;
  102. const center = options.center;
  103. // Computing the arc-length of the ellipse is too expensive to be practical. Estimating it using the
  104. // arc length of the sphere is too inaccurate and creates sharp edges when either the semi-major or
  105. // semi-minor axis is much bigger than the other. Instead, scale the angle delta to make
  106. // the distance along the ellipse boundary more closely match the granularity.
  107. const granularity = options.granularity * 8.0;
  108. const aSqr = semiMinorAxis * semiMinorAxis;
  109. const bSqr = semiMajorAxis * semiMajorAxis;
  110. const ab = semiMajorAxis * semiMinorAxis;
  111. const mag = Cartesian3.magnitude(center);
  112. const unitPos = Cartesian3.normalize(center, unitPosScratch);
  113. let eastVec = Cartesian3.cross(Cartesian3.UNIT_Z, center, eastVecScratch);
  114. eastVec = Cartesian3.normalize(eastVec, eastVec);
  115. const northVec = Cartesian3.cross(unitPos, eastVec, northVecScratch);
  116. // The number of points in the first quadrant
  117. let numPts = 1 + Math.ceil(CesiumMath.PI_OVER_TWO / granularity);
  118. const deltaTheta = CesiumMath.PI_OVER_TWO / (numPts - 1);
  119. let theta = CesiumMath.PI_OVER_TWO - numPts * deltaTheta;
  120. if (theta < 0.0) {
  121. numPts -= Math.ceil(Math.abs(theta) / deltaTheta);
  122. }
  123. // If the number of points were three, the ellipse
  124. // would be tessellated like below:
  125. //
  126. // *---*
  127. // / | \ | \
  128. // *---*---*---*
  129. // / | \ | \ | \ | \
  130. // / .*---*---*---*. \
  131. // * ` | \ | \ | \ | `*
  132. // \`.*---*---*---*.`/
  133. // \ | \ | \ | \ | /
  134. // *---*---*---*
  135. // \ | \ | /
  136. // *---*
  137. // The first and last column have one position and fan to connect to the adjacent column.
  138. // Each other vertical column contains an even number of positions.
  139. const size = 2 * (numPts * (numPts + 2));
  140. const positions = addFillPositions ? new Array(size * 3) : undefined;
  141. let positionIndex = 0;
  142. let position = scratchCartesian1;
  143. let reflectedPosition = scratchCartesian2;
  144. const outerPositionsLength = numPts * 4 * 3;
  145. let outerRightIndex = outerPositionsLength - 1;
  146. let outerLeftIndex = 0;
  147. const outerPositions = addEdgePositions
  148. ? new Array(outerPositionsLength)
  149. : undefined;
  150. let i;
  151. let j;
  152. let numInterior;
  153. let t;
  154. let interiorPosition;
  155. // Compute points in the 'eastern' half of the ellipse
  156. theta = CesiumMath.PI_OVER_TWO;
  157. position = pointOnEllipsoid(
  158. theta,
  159. rotation,
  160. northVec,
  161. eastVec,
  162. aSqr,
  163. ab,
  164. bSqr,
  165. mag,
  166. unitPos,
  167. position
  168. );
  169. if (addFillPositions) {
  170. positions[positionIndex++] = position.x;
  171. positions[positionIndex++] = position.y;
  172. positions[positionIndex++] = position.z;
  173. }
  174. if (addEdgePositions) {
  175. outerPositions[outerRightIndex--] = position.z;
  176. outerPositions[outerRightIndex--] = position.y;
  177. outerPositions[outerRightIndex--] = position.x;
  178. }
  179. theta = CesiumMath.PI_OVER_TWO - deltaTheta;
  180. for (i = 1; i < numPts + 1; ++i) {
  181. position = pointOnEllipsoid(
  182. theta,
  183. rotation,
  184. northVec,
  185. eastVec,
  186. aSqr,
  187. ab,
  188. bSqr,
  189. mag,
  190. unitPos,
  191. position
  192. );
  193. reflectedPosition = pointOnEllipsoid(
  194. Math.PI - theta,
  195. rotation,
  196. northVec,
  197. eastVec,
  198. aSqr,
  199. ab,
  200. bSqr,
  201. mag,
  202. unitPos,
  203. reflectedPosition
  204. );
  205. if (addFillPositions) {
  206. positions[positionIndex++] = position.x;
  207. positions[positionIndex++] = position.y;
  208. positions[positionIndex++] = position.z;
  209. numInterior = 2 * i + 2;
  210. for (j = 1; j < numInterior - 1; ++j) {
  211. t = j / (numInterior - 1);
  212. interiorPosition = Cartesian3.lerp(
  213. position,
  214. reflectedPosition,
  215. t,
  216. scratchCartesian3
  217. );
  218. positions[positionIndex++] = interiorPosition.x;
  219. positions[positionIndex++] = interiorPosition.y;
  220. positions[positionIndex++] = interiorPosition.z;
  221. }
  222. positions[positionIndex++] = reflectedPosition.x;
  223. positions[positionIndex++] = reflectedPosition.y;
  224. positions[positionIndex++] = reflectedPosition.z;
  225. }
  226. if (addEdgePositions) {
  227. outerPositions[outerRightIndex--] = position.z;
  228. outerPositions[outerRightIndex--] = position.y;
  229. outerPositions[outerRightIndex--] = position.x;
  230. outerPositions[outerLeftIndex++] = reflectedPosition.x;
  231. outerPositions[outerLeftIndex++] = reflectedPosition.y;
  232. outerPositions[outerLeftIndex++] = reflectedPosition.z;
  233. }
  234. theta = CesiumMath.PI_OVER_TWO - (i + 1) * deltaTheta;
  235. }
  236. // Compute points in the 'western' half of the ellipse
  237. for (i = numPts; i > 1; --i) {
  238. theta = CesiumMath.PI_OVER_TWO - (i - 1) * deltaTheta;
  239. position = pointOnEllipsoid(
  240. -theta,
  241. rotation,
  242. northVec,
  243. eastVec,
  244. aSqr,
  245. ab,
  246. bSqr,
  247. mag,
  248. unitPos,
  249. position
  250. );
  251. reflectedPosition = pointOnEllipsoid(
  252. theta + Math.PI,
  253. rotation,
  254. northVec,
  255. eastVec,
  256. aSqr,
  257. ab,
  258. bSqr,
  259. mag,
  260. unitPos,
  261. reflectedPosition
  262. );
  263. if (addFillPositions) {
  264. positions[positionIndex++] = position.x;
  265. positions[positionIndex++] = position.y;
  266. positions[positionIndex++] = position.z;
  267. numInterior = 2 * (i - 1) + 2;
  268. for (j = 1; j < numInterior - 1; ++j) {
  269. t = j / (numInterior - 1);
  270. interiorPosition = Cartesian3.lerp(
  271. position,
  272. reflectedPosition,
  273. t,
  274. scratchCartesian3
  275. );
  276. positions[positionIndex++] = interiorPosition.x;
  277. positions[positionIndex++] = interiorPosition.y;
  278. positions[positionIndex++] = interiorPosition.z;
  279. }
  280. positions[positionIndex++] = reflectedPosition.x;
  281. positions[positionIndex++] = reflectedPosition.y;
  282. positions[positionIndex++] = reflectedPosition.z;
  283. }
  284. if (addEdgePositions) {
  285. outerPositions[outerRightIndex--] = position.z;
  286. outerPositions[outerRightIndex--] = position.y;
  287. outerPositions[outerRightIndex--] = position.x;
  288. outerPositions[outerLeftIndex++] = reflectedPosition.x;
  289. outerPositions[outerLeftIndex++] = reflectedPosition.y;
  290. outerPositions[outerLeftIndex++] = reflectedPosition.z;
  291. }
  292. }
  293. theta = CesiumMath.PI_OVER_TWO;
  294. position = pointOnEllipsoid(
  295. -theta,
  296. rotation,
  297. northVec,
  298. eastVec,
  299. aSqr,
  300. ab,
  301. bSqr,
  302. mag,
  303. unitPos,
  304. position
  305. );
  306. const r = {};
  307. if (addFillPositions) {
  308. positions[positionIndex++] = position.x;
  309. positions[positionIndex++] = position.y;
  310. positions[positionIndex++] = position.z;
  311. r.positions = positions;
  312. r.numPts = numPts;
  313. }
  314. if (addEdgePositions) {
  315. outerPositions[outerRightIndex--] = position.z;
  316. outerPositions[outerRightIndex--] = position.y;
  317. outerPositions[outerRightIndex--] = position.x;
  318. r.outerPositions = outerPositions;
  319. }
  320. return r;
  321. };
  322. export default EllipseGeometryLibrary;