21302dac73bb8c92b2b9f06f939d32dc654b8ebba3e110391f3361272054231584a41460cee1e90ba8ddaab258f2309f4d8ce59f24a972012ed085d7dc254d 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500
  1. /**
  2. * 曲线辅助模块
  3. */
  4. import {
  5. create as v2Create,
  6. distSquare as v2DistSquare,
  7. VectorArray
  8. } from './vector';
  9. const mathPow = Math.pow;
  10. const mathSqrt = Math.sqrt;
  11. const EPSILON = 1e-8;
  12. const EPSILON_NUMERIC = 1e-4;
  13. const THREE_SQRT = mathSqrt(3);
  14. const ONE_THIRD = 1 / 3;
  15. // 临时变量
  16. const _v0 = v2Create();
  17. const _v1 = v2Create();
  18. const _v2 = v2Create();
  19. function isAroundZero(val: number) {
  20. return val > -EPSILON && val < EPSILON;
  21. }
  22. function isNotAroundZero(val: number) {
  23. return val > EPSILON || val < -EPSILON;
  24. }
  25. /**
  26. * 计算三次贝塞尔值
  27. */
  28. export function cubicAt(p0: number, p1: number, p2: number, p3: number, t: number): number {
  29. const onet = 1 - t;
  30. return onet * onet * (onet * p0 + 3 * t * p1)
  31. + t * t * (t * p3 + 3 * onet * p2);
  32. }
  33. /**
  34. * 计算三次贝塞尔导数值
  35. */
  36. export function cubicDerivativeAt(p0: number, p1: number, p2: number, p3: number, t: number): number {
  37. const onet = 1 - t;
  38. return 3 * (
  39. ((p1 - p0) * onet + 2 * (p2 - p1) * t) * onet
  40. + (p3 - p2) * t * t
  41. );
  42. }
  43. /**
  44. * 计算三次贝塞尔方程根,使用盛金公式
  45. */
  46. export function cubicRootAt(p0: number, p1: number, p2: number, p3: number, val: number, roots: number[]): number {
  47. // Evaluate roots of cubic functions
  48. const a = p3 + 3 * (p1 - p2) - p0;
  49. const b = 3 * (p2 - p1 * 2 + p0);
  50. const c = 3 * (p1 - p0);
  51. const d = p0 - val;
  52. const A = b * b - 3 * a * c;
  53. const B = b * c - 9 * a * d;
  54. const C = c * c - 3 * b * d;
  55. let n = 0;
  56. if (isAroundZero(A) && isAroundZero(B)) {
  57. if (isAroundZero(b)) {
  58. roots[0] = 0;
  59. }
  60. else {
  61. const t1 = -c / b; //t1, t2, t3, b is not zero
  62. if (t1 >= 0 && t1 <= 1) {
  63. roots[n++] = t1;
  64. }
  65. }
  66. }
  67. else {
  68. const disc = B * B - 4 * A * C;
  69. if (isAroundZero(disc)) {
  70. const K = B / A;
  71. const t1 = -b / a + K; // t1, a is not zero
  72. const t2 = -K / 2; // t2, t3
  73. if (t1 >= 0 && t1 <= 1) {
  74. roots[n++] = t1;
  75. }
  76. if (t2 >= 0 && t2 <= 1) {
  77. roots[n++] = t2;
  78. }
  79. }
  80. else if (disc > 0) {
  81. const discSqrt = mathSqrt(disc);
  82. let Y1 = A * b + 1.5 * a * (-B + discSqrt);
  83. let Y2 = A * b + 1.5 * a * (-B - discSqrt);
  84. if (Y1 < 0) {
  85. Y1 = -mathPow(-Y1, ONE_THIRD);
  86. }
  87. else {
  88. Y1 = mathPow(Y1, ONE_THIRD);
  89. }
  90. if (Y2 < 0) {
  91. Y2 = -mathPow(-Y2, ONE_THIRD);
  92. }
  93. else {
  94. Y2 = mathPow(Y2, ONE_THIRD);
  95. }
  96. const t1 = (-b - (Y1 + Y2)) / (3 * a);
  97. if (t1 >= 0 && t1 <= 1) {
  98. roots[n++] = t1;
  99. }
  100. }
  101. else {
  102. const T = (2 * A * b - 3 * a * B) / (2 * mathSqrt(A * A * A));
  103. const theta = Math.acos(T) / 3;
  104. const ASqrt = mathSqrt(A);
  105. const tmp = Math.cos(theta);
  106. const t1 = (-b - 2 * ASqrt * tmp) / (3 * a);
  107. const t2 = (-b + ASqrt * (tmp + THREE_SQRT * Math.sin(theta))) / (3 * a);
  108. const t3 = (-b + ASqrt * (tmp - THREE_SQRT * Math.sin(theta))) / (3 * a);
  109. if (t1 >= 0 && t1 <= 1) {
  110. roots[n++] = t1;
  111. }
  112. if (t2 >= 0 && t2 <= 1) {
  113. roots[n++] = t2;
  114. }
  115. if (t3 >= 0 && t3 <= 1) {
  116. roots[n++] = t3;
  117. }
  118. }
  119. }
  120. return n;
  121. }
  122. /**
  123. * 计算三次贝塞尔方程极限值的位置
  124. * @return 有效数目
  125. */
  126. export function cubicExtrema(p0: number, p1: number, p2: number, p3: number, extrema: number[]): number {
  127. const b = 6 * p2 - 12 * p1 + 6 * p0;
  128. const a = 9 * p1 + 3 * p3 - 3 * p0 - 9 * p2;
  129. const c = 3 * p1 - 3 * p0;
  130. let n = 0;
  131. if (isAroundZero(a)) {
  132. if (isNotAroundZero(b)) {
  133. const t1 = -c / b;
  134. if (t1 >= 0 && t1 <= 1) {
  135. extrema[n++] = t1;
  136. }
  137. }
  138. }
  139. else {
  140. const disc = b * b - 4 * a * c;
  141. if (isAroundZero(disc)) {
  142. extrema[0] = -b / (2 * a);
  143. }
  144. else if (disc > 0) {
  145. const discSqrt = mathSqrt(disc);
  146. const t1 = (-b + discSqrt) / (2 * a);
  147. const t2 = (-b - discSqrt) / (2 * a);
  148. if (t1 >= 0 && t1 <= 1) {
  149. extrema[n++] = t1;
  150. }
  151. if (t2 >= 0 && t2 <= 1) {
  152. extrema[n++] = t2;
  153. }
  154. }
  155. }
  156. return n;
  157. }
  158. /**
  159. * 细分三次贝塞尔曲线
  160. */
  161. export function cubicSubdivide(p0: number, p1: number, p2: number, p3: number, t: number, out: number[]) {
  162. const p01 = (p1 - p0) * t + p0;
  163. const p12 = (p2 - p1) * t + p1;
  164. const p23 = (p3 - p2) * t + p2;
  165. const p012 = (p12 - p01) * t + p01;
  166. const p123 = (p23 - p12) * t + p12;
  167. const p0123 = (p123 - p012) * t + p012;
  168. // Seg0
  169. out[0] = p0;
  170. out[1] = p01;
  171. out[2] = p012;
  172. out[3] = p0123;
  173. // Seg1
  174. out[4] = p0123;
  175. out[5] = p123;
  176. out[6] = p23;
  177. out[7] = p3;
  178. }
  179. /**
  180. * 投射点到三次贝塞尔曲线上,返回投射距离。
  181. * 投射点有可能会有一个或者多个,这里只返回其中距离最短的一个。
  182. */
  183. export function cubicProjectPoint(
  184. x0: number, y0: number, x1: number, y1: number, x2: number, y2: number, x3: number, y3: number,
  185. x: number, y: number, out: VectorArray
  186. ): number {
  187. // http://pomax.github.io/bezierinfo/#projections
  188. let t;
  189. let interval = 0.005;
  190. let d = Infinity;
  191. let prev;
  192. let next;
  193. let d1;
  194. let d2;
  195. _v0[0] = x;
  196. _v0[1] = y;
  197. // 先粗略估计一下可能的最小距离的 t 值
  198. // PENDING
  199. for (let _t = 0; _t < 1; _t += 0.05) {
  200. _v1[0] = cubicAt(x0, x1, x2, x3, _t);
  201. _v1[1] = cubicAt(y0, y1, y2, y3, _t);
  202. d1 = v2DistSquare(_v0, _v1);
  203. if (d1 < d) {
  204. t = _t;
  205. d = d1;
  206. }
  207. }
  208. d = Infinity;
  209. // At most 32 iteration
  210. for (let i = 0; i < 32; i++) {
  211. if (interval < EPSILON_NUMERIC) {
  212. break;
  213. }
  214. prev = t - interval;
  215. next = t + interval;
  216. // t - interval
  217. _v1[0] = cubicAt(x0, x1, x2, x3, prev);
  218. _v1[1] = cubicAt(y0, y1, y2, y3, prev);
  219. d1 = v2DistSquare(_v1, _v0);
  220. if (prev >= 0 && d1 < d) {
  221. t = prev;
  222. d = d1;
  223. }
  224. else {
  225. // t + interval
  226. _v2[0] = cubicAt(x0, x1, x2, x3, next);
  227. _v2[1] = cubicAt(y0, y1, y2, y3, next);
  228. d2 = v2DistSquare(_v2, _v0);
  229. if (next <= 1 && d2 < d) {
  230. t = next;
  231. d = d2;
  232. }
  233. else {
  234. interval *= 0.5;
  235. }
  236. }
  237. }
  238. // t
  239. if (out) {
  240. out[0] = cubicAt(x0, x1, x2, x3, t);
  241. out[1] = cubicAt(y0, y1, y2, y3, t);
  242. }
  243. // console.log(interval, i);
  244. return mathSqrt(d);
  245. }
  246. /**
  247. * 计算三次贝塞尔曲线长度
  248. */
  249. export function cubicLength(
  250. x0: number, y0: number, x1: number, y1: number, x2: number, y2: number, x3: number, y3: number,
  251. iteration: number
  252. ) {
  253. let px = x0;
  254. let py = y0;
  255. let d = 0;
  256. const step = 1 / iteration;
  257. for (let i = 1; i <= iteration; i++) {
  258. let t = i * step;
  259. const x = cubicAt(x0, x1, x2, x3, t);
  260. const y = cubicAt(y0, y1, y2, y3, t);
  261. const dx = x - px;
  262. const dy = y - py;
  263. d += Math.sqrt(dx * dx + dy * dy);
  264. px = x;
  265. py = y;
  266. }
  267. return d;
  268. }
  269. /**
  270. * 计算二次方贝塞尔值
  271. */
  272. export function quadraticAt(p0: number, p1: number, p2: number, t: number): number {
  273. const onet = 1 - t;
  274. return onet * (onet * p0 + 2 * t * p1) + t * t * p2;
  275. }
  276. /**
  277. * 计算二次方贝塞尔导数值
  278. */
  279. export function quadraticDerivativeAt(p0: number, p1: number, p2: number, t: number): number {
  280. return 2 * ((1 - t) * (p1 - p0) + t * (p2 - p1));
  281. }
  282. /**
  283. * 计算二次方贝塞尔方程根
  284. * @return 有效根数目
  285. */
  286. export function quadraticRootAt(p0: number, p1: number, p2: number, val: number, roots: number[]): number {
  287. const a = p0 - 2 * p1 + p2;
  288. const b = 2 * (p1 - p0);
  289. const c = p0 - val;
  290. let n = 0;
  291. if (isAroundZero(a)) {
  292. if (isNotAroundZero(b)) {
  293. const t1 = -c / b;
  294. if (t1 >= 0 && t1 <= 1) {
  295. roots[n++] = t1;
  296. }
  297. }
  298. }
  299. else {
  300. const disc = b * b - 4 * a * c;
  301. if (isAroundZero(disc)) {
  302. const t1 = -b / (2 * a);
  303. if (t1 >= 0 && t1 <= 1) {
  304. roots[n++] = t1;
  305. }
  306. }
  307. else if (disc > 0) {
  308. const discSqrt = mathSqrt(disc);
  309. const t1 = (-b + discSqrt) / (2 * a);
  310. const t2 = (-b - discSqrt) / (2 * a);
  311. if (t1 >= 0 && t1 <= 1) {
  312. roots[n++] = t1;
  313. }
  314. if (t2 >= 0 && t2 <= 1) {
  315. roots[n++] = t2;
  316. }
  317. }
  318. }
  319. return n;
  320. }
  321. /**
  322. * 计算二次贝塞尔方程极限值
  323. */
  324. export function quadraticExtremum(p0: number, p1: number, p2: number): number {
  325. const divider = p0 + p2 - 2 * p1;
  326. if (divider === 0) {
  327. // p1 is center of p0 and p2
  328. return 0.5;
  329. }
  330. else {
  331. return (p0 - p1) / divider;
  332. }
  333. }
  334. /**
  335. * 细分二次贝塞尔曲线
  336. */
  337. export function quadraticSubdivide(p0: number, p1: number, p2: number, t: number, out: number[]) {
  338. const p01 = (p1 - p0) * t + p0;
  339. const p12 = (p2 - p1) * t + p1;
  340. const p012 = (p12 - p01) * t + p01;
  341. // Seg0
  342. out[0] = p0;
  343. out[1] = p01;
  344. out[2] = p012;
  345. // Seg1
  346. out[3] = p012;
  347. out[4] = p12;
  348. out[5] = p2;
  349. }
  350. /**
  351. * 投射点到二次贝塞尔曲线上,返回投射距离。
  352. * 投射点有可能会有一个或者多个,这里只返回其中距离最短的一个。
  353. * @param {number} x0
  354. * @param {number} y0
  355. * @param {number} x1
  356. * @param {number} y1
  357. * @param {number} x2
  358. * @param {number} y2
  359. * @param {number} x
  360. * @param {number} y
  361. * @param {Array.<number>} out 投射点
  362. * @return {number}
  363. */
  364. export function quadraticProjectPoint(
  365. x0: number, y0: number, x1: number, y1: number, x2: number, y2: number,
  366. x: number, y: number, out: VectorArray
  367. ): number {
  368. // http://pomax.github.io/bezierinfo/#projections
  369. let t: number;
  370. let interval = 0.005;
  371. let d = Infinity;
  372. _v0[0] = x;
  373. _v0[1] = y;
  374. // 先粗略估计一下可能的最小距离的 t 值
  375. // PENDING
  376. for (let _t = 0; _t < 1; _t += 0.05) {
  377. _v1[0] = quadraticAt(x0, x1, x2, _t);
  378. _v1[1] = quadraticAt(y0, y1, y2, _t);
  379. const d1 = v2DistSquare(_v0, _v1);
  380. if (d1 < d) {
  381. t = _t;
  382. d = d1;
  383. }
  384. }
  385. d = Infinity;
  386. // At most 32 iteration
  387. for (let i = 0; i < 32; i++) {
  388. if (interval < EPSILON_NUMERIC) {
  389. break;
  390. }
  391. const prev = t - interval;
  392. const next = t + interval;
  393. // t - interval
  394. _v1[0] = quadraticAt(x0, x1, x2, prev);
  395. _v1[1] = quadraticAt(y0, y1, y2, prev);
  396. const d1 = v2DistSquare(_v1, _v0);
  397. if (prev >= 0 && d1 < d) {
  398. t = prev;
  399. d = d1;
  400. }
  401. else {
  402. // t + interval
  403. _v2[0] = quadraticAt(x0, x1, x2, next);
  404. _v2[1] = quadraticAt(y0, y1, y2, next);
  405. const d2 = v2DistSquare(_v2, _v0);
  406. if (next <= 1 && d2 < d) {
  407. t = next;
  408. d = d2;
  409. }
  410. else {
  411. interval *= 0.5;
  412. }
  413. }
  414. }
  415. // t
  416. if (out) {
  417. out[0] = quadraticAt(x0, x1, x2, t);
  418. out[1] = quadraticAt(y0, y1, y2, t);
  419. }
  420. // console.log(interval, i);
  421. return mathSqrt(d);
  422. }
  423. /**
  424. * 计算二次贝塞尔曲线长度
  425. */
  426. export function quadraticLength(
  427. x0: number, y0: number, x1: number, y1: number, x2: number, y2: number,
  428. iteration: number
  429. ) {
  430. let px = x0;
  431. let py = y0;
  432. let d = 0;
  433. const step = 1 / iteration;
  434. for (let i = 1; i <= iteration; i++) {
  435. let t = i * step;
  436. const x = quadraticAt(x0, x1, x2, t);
  437. const y = quadraticAt(y0, y1, y2, t);
  438. const dx = x - px;
  439. const dy = y - py;
  440. d += Math.sqrt(dx * dx + dy * dy);
  441. px = x;
  442. py = y;
  443. }
  444. return d;
  445. }