curve.js 9.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345
  1. import { create as v2Create, distSquare as v2DistSquare } from './vector.js';
  2. var mathPow = Math.pow;
  3. var mathSqrt = Math.sqrt;
  4. var EPSILON = 1e-8;
  5. var EPSILON_NUMERIC = 1e-4;
  6. var THREE_SQRT = mathSqrt(3);
  7. var ONE_THIRD = 1 / 3;
  8. var _v0 = v2Create();
  9. var _v1 = v2Create();
  10. var _v2 = v2Create();
  11. function isAroundZero(val) {
  12. return val > -EPSILON && val < EPSILON;
  13. }
  14. function isNotAroundZero(val) {
  15. return val > EPSILON || val < -EPSILON;
  16. }
  17. export function cubicAt(p0, p1, p2, p3, t) {
  18. var onet = 1 - t;
  19. return onet * onet * (onet * p0 + 3 * t * p1)
  20. + t * t * (t * p3 + 3 * onet * p2);
  21. }
  22. export function cubicDerivativeAt(p0, p1, p2, p3, t) {
  23. var onet = 1 - t;
  24. return 3 * (((p1 - p0) * onet + 2 * (p2 - p1) * t) * onet
  25. + (p3 - p2) * t * t);
  26. }
  27. export function cubicRootAt(p0, p1, p2, p3, val, roots) {
  28. var a = p3 + 3 * (p1 - p2) - p0;
  29. var b = 3 * (p2 - p1 * 2 + p0);
  30. var c = 3 * (p1 - p0);
  31. var d = p0 - val;
  32. var A = b * b - 3 * a * c;
  33. var B = b * c - 9 * a * d;
  34. var C = c * c - 3 * b * d;
  35. var n = 0;
  36. if (isAroundZero(A) && isAroundZero(B)) {
  37. if (isAroundZero(b)) {
  38. roots[0] = 0;
  39. }
  40. else {
  41. var t1 = -c / b;
  42. if (t1 >= 0 && t1 <= 1) {
  43. roots[n++] = t1;
  44. }
  45. }
  46. }
  47. else {
  48. var disc = B * B - 4 * A * C;
  49. if (isAroundZero(disc)) {
  50. var K = B / A;
  51. var t1 = -b / a + K;
  52. var t2 = -K / 2;
  53. if (t1 >= 0 && t1 <= 1) {
  54. roots[n++] = t1;
  55. }
  56. if (t2 >= 0 && t2 <= 1) {
  57. roots[n++] = t2;
  58. }
  59. }
  60. else if (disc > 0) {
  61. var discSqrt = mathSqrt(disc);
  62. var Y1 = A * b + 1.5 * a * (-B + discSqrt);
  63. var Y2 = A * b + 1.5 * a * (-B - discSqrt);
  64. if (Y1 < 0) {
  65. Y1 = -mathPow(-Y1, ONE_THIRD);
  66. }
  67. else {
  68. Y1 = mathPow(Y1, ONE_THIRD);
  69. }
  70. if (Y2 < 0) {
  71. Y2 = -mathPow(-Y2, ONE_THIRD);
  72. }
  73. else {
  74. Y2 = mathPow(Y2, ONE_THIRD);
  75. }
  76. var t1 = (-b - (Y1 + Y2)) / (3 * a);
  77. if (t1 >= 0 && t1 <= 1) {
  78. roots[n++] = t1;
  79. }
  80. }
  81. else {
  82. var T = (2 * A * b - 3 * a * B) / (2 * mathSqrt(A * A * A));
  83. var theta = Math.acos(T) / 3;
  84. var ASqrt = mathSqrt(A);
  85. var tmp = Math.cos(theta);
  86. var t1 = (-b - 2 * ASqrt * tmp) / (3 * a);
  87. var t2 = (-b + ASqrt * (tmp + THREE_SQRT * Math.sin(theta))) / (3 * a);
  88. var t3 = (-b + ASqrt * (tmp - THREE_SQRT * Math.sin(theta))) / (3 * a);
  89. if (t1 >= 0 && t1 <= 1) {
  90. roots[n++] = t1;
  91. }
  92. if (t2 >= 0 && t2 <= 1) {
  93. roots[n++] = t2;
  94. }
  95. if (t3 >= 0 && t3 <= 1) {
  96. roots[n++] = t3;
  97. }
  98. }
  99. }
  100. return n;
  101. }
  102. export function cubicExtrema(p0, p1, p2, p3, extrema) {
  103. var b = 6 * p2 - 12 * p1 + 6 * p0;
  104. var a = 9 * p1 + 3 * p3 - 3 * p0 - 9 * p2;
  105. var c = 3 * p1 - 3 * p0;
  106. var n = 0;
  107. if (isAroundZero(a)) {
  108. if (isNotAroundZero(b)) {
  109. var t1 = -c / b;
  110. if (t1 >= 0 && t1 <= 1) {
  111. extrema[n++] = t1;
  112. }
  113. }
  114. }
  115. else {
  116. var disc = b * b - 4 * a * c;
  117. if (isAroundZero(disc)) {
  118. extrema[0] = -b / (2 * a);
  119. }
  120. else if (disc > 0) {
  121. var discSqrt = mathSqrt(disc);
  122. var t1 = (-b + discSqrt) / (2 * a);
  123. var t2 = (-b - discSqrt) / (2 * a);
  124. if (t1 >= 0 && t1 <= 1) {
  125. extrema[n++] = t1;
  126. }
  127. if (t2 >= 0 && t2 <= 1) {
  128. extrema[n++] = t2;
  129. }
  130. }
  131. }
  132. return n;
  133. }
  134. export function cubicSubdivide(p0, p1, p2, p3, t, out) {
  135. var p01 = (p1 - p0) * t + p0;
  136. var p12 = (p2 - p1) * t + p1;
  137. var p23 = (p3 - p2) * t + p2;
  138. var p012 = (p12 - p01) * t + p01;
  139. var p123 = (p23 - p12) * t + p12;
  140. var p0123 = (p123 - p012) * t + p012;
  141. out[0] = p0;
  142. out[1] = p01;
  143. out[2] = p012;
  144. out[3] = p0123;
  145. out[4] = p0123;
  146. out[5] = p123;
  147. out[6] = p23;
  148. out[7] = p3;
  149. }
  150. export function cubicProjectPoint(x0, y0, x1, y1, x2, y2, x3, y3, x, y, out) {
  151. var t;
  152. var interval = 0.005;
  153. var d = Infinity;
  154. var prev;
  155. var next;
  156. var d1;
  157. var d2;
  158. _v0[0] = x;
  159. _v0[1] = y;
  160. for (var _t = 0; _t < 1; _t += 0.05) {
  161. _v1[0] = cubicAt(x0, x1, x2, x3, _t);
  162. _v1[1] = cubicAt(y0, y1, y2, y3, _t);
  163. d1 = v2DistSquare(_v0, _v1);
  164. if (d1 < d) {
  165. t = _t;
  166. d = d1;
  167. }
  168. }
  169. d = Infinity;
  170. for (var i = 0; i < 32; i++) {
  171. if (interval < EPSILON_NUMERIC) {
  172. break;
  173. }
  174. prev = t - interval;
  175. next = t + interval;
  176. _v1[0] = cubicAt(x0, x1, x2, x3, prev);
  177. _v1[1] = cubicAt(y0, y1, y2, y3, prev);
  178. d1 = v2DistSquare(_v1, _v0);
  179. if (prev >= 0 && d1 < d) {
  180. t = prev;
  181. d = d1;
  182. }
  183. else {
  184. _v2[0] = cubicAt(x0, x1, x2, x3, next);
  185. _v2[1] = cubicAt(y0, y1, y2, y3, next);
  186. d2 = v2DistSquare(_v2, _v0);
  187. if (next <= 1 && d2 < d) {
  188. t = next;
  189. d = d2;
  190. }
  191. else {
  192. interval *= 0.5;
  193. }
  194. }
  195. }
  196. if (out) {
  197. out[0] = cubicAt(x0, x1, x2, x3, t);
  198. out[1] = cubicAt(y0, y1, y2, y3, t);
  199. }
  200. return mathSqrt(d);
  201. }
  202. export function cubicLength(x0, y0, x1, y1, x2, y2, x3, y3, iteration) {
  203. var px = x0;
  204. var py = y0;
  205. var d = 0;
  206. var step = 1 / iteration;
  207. for (var i = 1; i <= iteration; i++) {
  208. var t = i * step;
  209. var x = cubicAt(x0, x1, x2, x3, t);
  210. var y = cubicAt(y0, y1, y2, y3, t);
  211. var dx = x - px;
  212. var dy = y - py;
  213. d += Math.sqrt(dx * dx + dy * dy);
  214. px = x;
  215. py = y;
  216. }
  217. return d;
  218. }
  219. export function quadraticAt(p0, p1, p2, t) {
  220. var onet = 1 - t;
  221. return onet * (onet * p0 + 2 * t * p1) + t * t * p2;
  222. }
  223. export function quadraticDerivativeAt(p0, p1, p2, t) {
  224. return 2 * ((1 - t) * (p1 - p0) + t * (p2 - p1));
  225. }
  226. export function quadraticRootAt(p0, p1, p2, val, roots) {
  227. var a = p0 - 2 * p1 + p2;
  228. var b = 2 * (p1 - p0);
  229. var c = p0 - val;
  230. var n = 0;
  231. if (isAroundZero(a)) {
  232. if (isNotAroundZero(b)) {
  233. var t1 = -c / b;
  234. if (t1 >= 0 && t1 <= 1) {
  235. roots[n++] = t1;
  236. }
  237. }
  238. }
  239. else {
  240. var disc = b * b - 4 * a * c;
  241. if (isAroundZero(disc)) {
  242. var t1 = -b / (2 * a);
  243. if (t1 >= 0 && t1 <= 1) {
  244. roots[n++] = t1;
  245. }
  246. }
  247. else if (disc > 0) {
  248. var discSqrt = mathSqrt(disc);
  249. var t1 = (-b + discSqrt) / (2 * a);
  250. var t2 = (-b - discSqrt) / (2 * a);
  251. if (t1 >= 0 && t1 <= 1) {
  252. roots[n++] = t1;
  253. }
  254. if (t2 >= 0 && t2 <= 1) {
  255. roots[n++] = t2;
  256. }
  257. }
  258. }
  259. return n;
  260. }
  261. export function quadraticExtremum(p0, p1, p2) {
  262. var divider = p0 + p2 - 2 * p1;
  263. if (divider === 0) {
  264. return 0.5;
  265. }
  266. else {
  267. return (p0 - p1) / divider;
  268. }
  269. }
  270. export function quadraticSubdivide(p0, p1, p2, t, out) {
  271. var p01 = (p1 - p0) * t + p0;
  272. var p12 = (p2 - p1) * t + p1;
  273. var p012 = (p12 - p01) * t + p01;
  274. out[0] = p0;
  275. out[1] = p01;
  276. out[2] = p012;
  277. out[3] = p012;
  278. out[4] = p12;
  279. out[5] = p2;
  280. }
  281. export function quadraticProjectPoint(x0, y0, x1, y1, x2, y2, x, y, out) {
  282. var t;
  283. var interval = 0.005;
  284. var d = Infinity;
  285. _v0[0] = x;
  286. _v0[1] = y;
  287. for (var _t = 0; _t < 1; _t += 0.05) {
  288. _v1[0] = quadraticAt(x0, x1, x2, _t);
  289. _v1[1] = quadraticAt(y0, y1, y2, _t);
  290. var d1 = v2DistSquare(_v0, _v1);
  291. if (d1 < d) {
  292. t = _t;
  293. d = d1;
  294. }
  295. }
  296. d = Infinity;
  297. for (var i = 0; i < 32; i++) {
  298. if (interval < EPSILON_NUMERIC) {
  299. break;
  300. }
  301. var prev = t - interval;
  302. var next = t + interval;
  303. _v1[0] = quadraticAt(x0, x1, x2, prev);
  304. _v1[1] = quadraticAt(y0, y1, y2, prev);
  305. var d1 = v2DistSquare(_v1, _v0);
  306. if (prev >= 0 && d1 < d) {
  307. t = prev;
  308. d = d1;
  309. }
  310. else {
  311. _v2[0] = quadraticAt(x0, x1, x2, next);
  312. _v2[1] = quadraticAt(y0, y1, y2, next);
  313. var d2 = v2DistSquare(_v2, _v0);
  314. if (next <= 1 && d2 < d) {
  315. t = next;
  316. d = d2;
  317. }
  318. else {
  319. interval *= 0.5;
  320. }
  321. }
  322. }
  323. if (out) {
  324. out[0] = quadraticAt(x0, x1, x2, t);
  325. out[1] = quadraticAt(y0, y1, y2, t);
  326. }
  327. return mathSqrt(d);
  328. }
  329. export function quadraticLength(x0, y0, x1, y1, x2, y2, iteration) {
  330. var px = x0;
  331. var py = y0;
  332. var d = 0;
  333. var step = 1 / iteration;
  334. for (var i = 1; i <= iteration; i++) {
  335. var t = i * step;
  336. var x = quadraticAt(x0, x1, x2, t);
  337. var y = quadraticAt(y0, y1, y2, t);
  338. var dx = x - px;
  339. var dy = y - py;
  340. d += Math.sqrt(dx * dx + dy * dy);
  341. px = x;
  342. py = y;
  343. }
  344. return d;
  345. }