HermitePolynomialApproximation.js 9.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342
  1. import defaultValue from "./defaultValue.js";
  2. import defined from "./defined.js";
  3. import DeveloperError from "./DeveloperError.js";
  4. import CesiumMath from "./Math.js";
  5. const factorial = CesiumMath.factorial;
  6. function calculateCoefficientTerm(
  7. x,
  8. zIndices,
  9. xTable,
  10. derivOrder,
  11. termOrder,
  12. reservedIndices
  13. ) {
  14. let result = 0;
  15. let reserved;
  16. let i;
  17. let j;
  18. if (derivOrder > 0) {
  19. for (i = 0; i < termOrder; i++) {
  20. reserved = false;
  21. for (j = 0; j < reservedIndices.length && !reserved; j++) {
  22. if (i === reservedIndices[j]) {
  23. reserved = true;
  24. }
  25. }
  26. if (!reserved) {
  27. reservedIndices.push(i);
  28. result += calculateCoefficientTerm(
  29. x,
  30. zIndices,
  31. xTable,
  32. derivOrder - 1,
  33. termOrder,
  34. reservedIndices
  35. );
  36. reservedIndices.splice(reservedIndices.length - 1, 1);
  37. }
  38. }
  39. return result;
  40. }
  41. result = 1;
  42. for (i = 0; i < termOrder; i++) {
  43. reserved = false;
  44. for (j = 0; j < reservedIndices.length && !reserved; j++) {
  45. if (i === reservedIndices[j]) {
  46. reserved = true;
  47. }
  48. }
  49. if (!reserved) {
  50. result *= x - xTable[zIndices[i]];
  51. }
  52. }
  53. return result;
  54. }
  55. /**
  56. * An {@link InterpolationAlgorithm} for performing Hermite interpolation.
  57. *
  58. * @namespace HermitePolynomialApproximation
  59. */
  60. const HermitePolynomialApproximation = {
  61. type: "Hermite",
  62. };
  63. /**
  64. * Given the desired degree, returns the number of data points required for interpolation.
  65. *
  66. * @param {number} degree The desired degree of interpolation.
  67. * @param {number} [inputOrder=0] The order of the inputs (0 means just the data, 1 means the data and its derivative, etc).
  68. * @returns {number} The number of required data points needed for the desired degree of interpolation.
  69. * @exception {DeveloperError} degree must be 0 or greater.
  70. * @exception {DeveloperError} inputOrder must be 0 or greater.
  71. */
  72. HermitePolynomialApproximation.getRequiredDataPoints = function (
  73. degree,
  74. inputOrder
  75. ) {
  76. inputOrder = defaultValue(inputOrder, 0);
  77. //>>includeStart('debug', pragmas.debug);
  78. if (!defined(degree)) {
  79. throw new DeveloperError("degree is required.");
  80. }
  81. if (degree < 0) {
  82. throw new DeveloperError("degree must be 0 or greater.");
  83. }
  84. if (inputOrder < 0) {
  85. throw new DeveloperError("inputOrder must be 0 or greater.");
  86. }
  87. //>>includeEnd('debug');
  88. return Math.max(Math.floor((degree + 1) / (inputOrder + 1)), 2);
  89. };
  90. /**
  91. * Interpolates values using Hermite Polynomial Approximation.
  92. *
  93. * @param {number} x The independent variable for which the dependent variables will be interpolated.
  94. * @param {number[]} xTable The array of independent variables to use to interpolate. The values
  95. * in this array must be in increasing order and the same value must not occur twice in the array.
  96. * @param {number[]} yTable The array of dependent variables to use to interpolate. For a set of three
  97. * dependent values (p,q,w) at time 1 and time 2 this should be as follows: {p1, q1, w1, p2, q2, w2}.
  98. * @param {number} yStride The number of dependent variable values in yTable corresponding to
  99. * each independent variable value in xTable.
  100. * @param {number[]} [result] An existing array into which to store the result.
  101. * @returns {number[]} The array of interpolated values, or the result parameter if one was provided.
  102. */
  103. HermitePolynomialApproximation.interpolateOrderZero = function (
  104. x,
  105. xTable,
  106. yTable,
  107. yStride,
  108. result
  109. ) {
  110. if (!defined(result)) {
  111. result = new Array(yStride);
  112. }
  113. let i;
  114. let j;
  115. let d;
  116. let s;
  117. let len;
  118. let index;
  119. const length = xTable.length;
  120. const coefficients = new Array(yStride);
  121. for (i = 0; i < yStride; i++) {
  122. result[i] = 0;
  123. const l = new Array(length);
  124. coefficients[i] = l;
  125. for (j = 0; j < length; j++) {
  126. l[j] = [];
  127. }
  128. }
  129. const zIndicesLength = length,
  130. zIndices = new Array(zIndicesLength);
  131. for (i = 0; i < zIndicesLength; i++) {
  132. zIndices[i] = i;
  133. }
  134. let highestNonZeroCoef = length - 1;
  135. for (s = 0; s < yStride; s++) {
  136. for (j = 0; j < zIndicesLength; j++) {
  137. index = zIndices[j] * yStride + s;
  138. coefficients[s][0].push(yTable[index]);
  139. }
  140. for (i = 1; i < zIndicesLength; i++) {
  141. let nonZeroCoefficients = false;
  142. for (j = 0; j < zIndicesLength - i; j++) {
  143. const zj = xTable[zIndices[j]];
  144. const zn = xTable[zIndices[j + i]];
  145. let numerator;
  146. if (zn - zj <= 0) {
  147. index = zIndices[j] * yStride + yStride * i + s;
  148. numerator = yTable[index];
  149. coefficients[s][i].push(numerator / factorial(i));
  150. } else {
  151. numerator = coefficients[s][i - 1][j + 1] - coefficients[s][i - 1][j];
  152. coefficients[s][i].push(numerator / (zn - zj));
  153. }
  154. nonZeroCoefficients = nonZeroCoefficients || numerator !== 0;
  155. }
  156. if (!nonZeroCoefficients) {
  157. highestNonZeroCoef = i - 1;
  158. }
  159. }
  160. }
  161. for (d = 0, len = 0; d <= len; d++) {
  162. for (i = d; i <= highestNonZeroCoef; i++) {
  163. const tempTerm = calculateCoefficientTerm(x, zIndices, xTable, d, i, []);
  164. for (s = 0; s < yStride; s++) {
  165. const coeff = coefficients[s][i][0];
  166. result[s + d * yStride] += coeff * tempTerm;
  167. }
  168. }
  169. }
  170. return result;
  171. };
  172. const arrayScratch = [];
  173. /**
  174. * Interpolates values using Hermite Polynomial Approximation.
  175. *
  176. * @param {number} x The independent variable for which the dependent variables will be interpolated.
  177. * @param {number[]} xTable The array of independent variables to use to interpolate. The values
  178. * in this array must be in increasing order and the same value must not occur twice in the array.
  179. * @param {number[]} yTable The array of dependent variables to use to interpolate. For a set of three
  180. * dependent values (p,q,w) at time 1 and time 2 this should be as follows: {p1, q1, w1, p2, q2, w2}.
  181. * @param {number} yStride The number of dependent variable values in yTable corresponding to
  182. * each independent variable value in xTable.
  183. * @param {number} inputOrder The number of derivatives supplied for input.
  184. * @param {number} outputOrder The number of derivatives desired for output.
  185. * @param {number[]} [result] An existing array into which to store the result.
  186. *
  187. * @returns {number[]} The array of interpolated values, or the result parameter if one was provided.
  188. */
  189. HermitePolynomialApproximation.interpolate = function (
  190. x,
  191. xTable,
  192. yTable,
  193. yStride,
  194. inputOrder,
  195. outputOrder,
  196. result
  197. ) {
  198. const resultLength = yStride * (outputOrder + 1);
  199. if (!defined(result)) {
  200. result = new Array(resultLength);
  201. }
  202. for (let r = 0; r < resultLength; r++) {
  203. result[r] = 0;
  204. }
  205. const length = xTable.length;
  206. // The zIndices array holds copies of the addresses of the xTable values
  207. // in the range we're looking at. Even though this just holds information already
  208. // available in xTable this is a much more convenient format.
  209. const zIndices = new Array(length * (inputOrder + 1));
  210. let i;
  211. for (i = 0; i < length; i++) {
  212. for (let j = 0; j < inputOrder + 1; j++) {
  213. zIndices[i * (inputOrder + 1) + j] = i;
  214. }
  215. }
  216. const zIndiceslength = zIndices.length;
  217. const coefficients = arrayScratch;
  218. const highestNonZeroCoef = fillCoefficientList(
  219. coefficients,
  220. zIndices,
  221. xTable,
  222. yTable,
  223. yStride,
  224. inputOrder
  225. );
  226. const reservedIndices = [];
  227. const tmp = (zIndiceslength * (zIndiceslength + 1)) / 2;
  228. const loopStop = Math.min(highestNonZeroCoef, outputOrder);
  229. for (let d = 0; d <= loopStop; d++) {
  230. for (i = d; i <= highestNonZeroCoef; i++) {
  231. reservedIndices.length = 0;
  232. const tempTerm = calculateCoefficientTerm(
  233. x,
  234. zIndices,
  235. xTable,
  236. d,
  237. i,
  238. reservedIndices
  239. );
  240. const dimTwo = Math.floor((i * (1 - i)) / 2) + zIndiceslength * i;
  241. for (let s = 0; s < yStride; s++) {
  242. const dimOne = Math.floor(s * tmp);
  243. const coef = coefficients[dimOne + dimTwo];
  244. result[s + d * yStride] += coef * tempTerm;
  245. }
  246. }
  247. }
  248. return result;
  249. };
  250. function fillCoefficientList(
  251. coefficients,
  252. zIndices,
  253. xTable,
  254. yTable,
  255. yStride,
  256. inputOrder
  257. ) {
  258. let j;
  259. let index;
  260. let highestNonZero = -1;
  261. const zIndiceslength = zIndices.length;
  262. const tmp = (zIndiceslength * (zIndiceslength + 1)) / 2;
  263. for (let s = 0; s < yStride; s++) {
  264. const dimOne = Math.floor(s * tmp);
  265. for (j = 0; j < zIndiceslength; j++) {
  266. index = zIndices[j] * yStride * (inputOrder + 1) + s;
  267. coefficients[dimOne + j] = yTable[index];
  268. }
  269. for (let i = 1; i < zIndiceslength; i++) {
  270. let coefIndex = 0;
  271. const dimTwo = Math.floor((i * (1 - i)) / 2) + zIndiceslength * i;
  272. let nonZeroCoefficients = false;
  273. for (j = 0; j < zIndiceslength - i; j++) {
  274. const zj = xTable[zIndices[j]];
  275. const zn = xTable[zIndices[j + i]];
  276. let numerator;
  277. let coefficient;
  278. if (zn - zj <= 0) {
  279. index = zIndices[j] * yStride * (inputOrder + 1) + yStride * i + s;
  280. numerator = yTable[index];
  281. coefficient = numerator / CesiumMath.factorial(i);
  282. coefficients[dimOne + dimTwo + coefIndex] = coefficient;
  283. coefIndex++;
  284. } else {
  285. const dimTwoMinusOne =
  286. Math.floor(((i - 1) * (2 - i)) / 2) + zIndiceslength * (i - 1);
  287. numerator =
  288. coefficients[dimOne + dimTwoMinusOne + j + 1] -
  289. coefficients[dimOne + dimTwoMinusOne + j];
  290. coefficient = numerator / (zn - zj);
  291. coefficients[dimOne + dimTwo + coefIndex] = coefficient;
  292. coefIndex++;
  293. }
  294. nonZeroCoefficients = nonZeroCoefficients || numerator !== 0.0;
  295. }
  296. if (nonZeroCoefficients) {
  297. highestNonZero = Math.max(highestNonZero, i);
  298. }
  299. }
  300. }
  301. return highestNonZero;
  302. }
  303. export default HermitePolynomialApproximation;