1 /* 2 Copyright 2008-2023 3 Matthias Ehmann, 4 Michael Gerhaeuser, 5 Carsten Miller, 6 Bianca Valentin, 7 Alfred Wassermann, 8 Peter Wilfahrt 9 10 This file is part of JSXGraph. 11 12 JSXGraph is free software dual licensed under the GNU LGPL or MIT License. 13 14 You can redistribute it and/or modify it under the terms of the 15 16 * GNU Lesser General Public License as published by 17 the Free Software Foundation, either version 3 of the License, or 18 (at your option) any later version 19 OR 20 * MIT License: https://github.com/jsxgraph/jsxgraph/blob/master/LICENSE.MIT 21 22 JSXGraph is distributed in the hope that it will be useful, 23 but WITHOUT ANY WARRANTY; without even the implied warranty of 24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 25 GNU Lesser General Public License for more details. 26 27 You should have received a copy of the GNU Lesser General Public License and 28 the MIT License along with JSXGraph. If not, see <https://www.gnu.org/licenses/> 29 and <https://opensource.org/licenses/MIT/>. 30 31 32 Metapost/Hobby curves, see e.g. https://bosker.wordpress.com/2013/11/13/beyond-bezier-curves/ 33 34 * Ported to Python for the project PyX. Copyright (C) 2011 Michael Schindler <m-schindler@users.sourceforge.net> 35 * Ported to javascript from the PyX implementation (https://pyx-project.org/) by Vlad-X. 36 * Adapted to JSXGraph and some code changes by Alfred Wassermann 2020. 37 38 This program is distributed in the hope that it will be useful, 39 but WITHOUT ANY WARRANTY; without even the implied warranty of 40 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 41 GNU General Public License for more details. 42 43 You should have received a copy of the GNU General Public License 44 along with this program; if not, write to the Free Software 45 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 46 47 Internal functions of MetaPost 48 This file re-implements some of the functionality of MetaPost 49 (https://tug.org/metapost.html). MetaPost was developed by John D. Hobby and 50 others. The code of Metapost is in the public domain, which we understand as 51 an implicit permission to reuse the code here (see the comment at 52 https://www.gnu.org/licenses/license-list.html) 53 54 This file is based on the MetaPost version distributed by TeXLive: 55 svn://tug.org/texlive/trunk/Build/source/texk/web2c/mplibdir revision 22737 # 56 (2011-05-31) 57 */ 58 59 /*global JXG: true, define: true*/ 60 /*jslint nomen: true, plusplus: true*/ 61 62 /** 63 * @fileoverview In this file the namespace Math.Metapost is defined which holds algorithms translated from Metapost 64 * by D.E. Knuth and J.D. Hobby. 65 */ 66 67 import Type from "../utils/type"; 68 import Mat from "./math"; 69 70 /** 71 * The JXG.Math.Metapost namespace holds algorithms translated from Metapost 72 * by D.E. Knuth and J.D. Hobby. 73 * 74 * @name JXG.Math.Metapost 75 * @exports Mat.Metapost as JXG.Math.Metapost 76 * @namespace 77 */ 78 Mat.Metapost = { 79 MP_ENDPOINT: 0, 80 MP_EXPLICIT: 1, 81 MP_GIVEN: 2, 82 MP_CURL: 3, 83 MP_OPEN: 4, 84 MP_END_CYCLE: 5, 85 86 UNITY: 1.0, 87 // two: 2, 88 // fraction_half: 0.5, 89 FRACTION_ONE: 1.0, 90 FRACTION_THREE: 3.0, 91 ONE_EIGHTY_DEG: Math.PI, 92 THREE_SIXTY_DEG: 2 * Math.PI, 93 // EPSILON: 1e-5, 94 EPS_SQ: 1e-5 * 1e-5, 95 96 /** 97 * @private 98 */ 99 make_choices: function (knots) { 100 var dely, h, k, delx, n, q, p, s, cosine, t, sine, delta_x, delta_y, delta, psi, 101 endless = true; 102 103 p = knots[0]; 104 do { 105 if (!p) { 106 break; 107 } 108 q = p.next; 109 110 // Join two identical knots by setting the control points to the same 111 // coordinates. 112 // MP 291 113 if ( 114 p.rtype > this.MP_EXPLICIT && 115 (p.x - q.x) * (p.x - q.x) + (p.y - q.y) * (p.y - q.y) < this.EPS_SQ 116 ) { 117 p.rtype = this.MP_EXPLICIT; 118 if (p.ltype === this.MP_OPEN) { 119 p.ltype = this.MP_CURL; 120 p.set_left_curl(this.UNITY); 121 } 122 123 q.ltype = this.MP_EXPLICIT; 124 if (q.rtype === this.MP_OPEN) { 125 q.rtype = this.MP_CURL; 126 q.set_right_curl(this.UNITY); 127 } 128 129 p.rx = p.x; 130 q.lx = p.x; 131 p.ry = p.y; 132 q.ly = p.y; 133 } 134 p = q; 135 } while (p !== knots[0]); 136 137 // Find the first breakpoint, h, on the path 138 // MP 292 139 h = knots[0]; 140 while (endless) { 141 if (h.ltype !== this.MP_OPEN || h.rtype !== this.MP_OPEN) { 142 break; 143 } 144 h = h.next; 145 if (h === knots[0]) { 146 h.ltype = this.MP_END_CYCLE; 147 break; 148 } 149 } 150 151 p = h; 152 while (endless) { 153 if (!p) { 154 break; 155 } 156 157 // Fill in the control points between p and the next breakpoint, 158 // then advance p to that breakpoint 159 // MP 299 160 q = p.next; 161 if (p.rtype >= this.MP_GIVEN) { 162 while (q.ltype === this.MP_OPEN && q.rtype === this.MP_OPEN) { 163 q = q.next; 164 } 165 166 // Calculate the turning angles psi_ k and the distances d_{k,k+1}; 167 // set n to the length of the path 168 // MP 302 169 k = 0; 170 s = p; 171 n = knots.length; 172 173 delta_x = []; 174 delta_y = []; 175 delta = []; 176 psi = [null]; 177 178 // tuple([]) = tuple([[], [], [], [null]]); 179 while (endless) { 180 t = s.next; 181 // None; 182 delta_x.push(t.x - s.x); 183 delta_y.push(t.y - s.y); 184 delta.push(this.mp_pyth_add(delta_x[k], delta_y[k])); 185 if (k > 0) { 186 sine = delta_y[k - 1] / delta[k - 1]; 187 cosine = delta_x[k - 1] / delta[k - 1]; 188 psi.push( 189 Math.atan2( 190 delta_y[k] * cosine - delta_x[k] * sine, 191 delta_x[k] * cosine + delta_y[k] * sine 192 ) 193 ); 194 } 195 k++; 196 s = t; 197 if (s === q) { 198 n = k; 199 } 200 if (k >= n && s.ltype !== this.MP_END_CYCLE) { 201 break; 202 } 203 } 204 if (k === n) { 205 psi.push(0); 206 } else { 207 psi.push(psi[1]); 208 } 209 210 // Remove open types at the breakpoints 211 // MP 303 212 if (q.ltype === this.MP_OPEN) { 213 delx = q.rx - q.x; 214 dely = q.ry - q.y; 215 if (delx * delx + dely * dely < this.EPS_SQ) { 216 q.ltype = this.MP_CURL; 217 q.set_left_curl(this.UNITY); 218 } else { 219 q.ltype = this.MP_GIVEN; 220 q.set_left_given(Math.atan2(dely, delx)); 221 } 222 } 223 if (p.rtype === this.MP_OPEN && p.ltype === this.MP_EXPLICIT) { 224 delx = p.x - p.lx; 225 dely = p.y - p.ly; 226 if (delx * delx + dely * dely < this.EPS_SQ) { 227 p.rtype = this.MP_CURL; 228 p.set_right_curl(this.UNITY); 229 } else { 230 p.rtype = this.MP_GIVEN; 231 p.set_right_given(Math.atan2(dely, delx)); 232 } 233 } 234 this.mp_solve_choices(p, q, n, delta_x, delta_y, delta, psi); 235 } else if (p.rtype === this.MP_ENDPOINT) { 236 // MP 294 237 p.rx = p.x; 238 p.ry = p.y; 239 q.lx = q.x; 240 q.ly = q.y; 241 } 242 p = q; 243 244 if (p === h) { 245 break; 246 } 247 } 248 }, 249 250 /** 251 * Implements solve_choices form metapost 252 * MP 305 253 * @private 254 */ 255 mp_solve_choices: function (p, q, n, delta_x, delta_y, delta, psi) { 256 var aa, acc, vv, bb, ldelta, ee, k, 257 s, ww, uu, lt, r, t, ff, 258 theta, rt, dd, cc, ct_st, 259 ct, st, cf_sf, cf, sf, i, k_idx, 260 endless = true; 261 262 ldelta = delta.length + 1; 263 uu = new Array(ldelta); 264 ww = new Array(ldelta); 265 vv = new Array(ldelta); 266 theta = new Array(ldelta); 267 for (i = 0; i < ldelta; i++) { 268 theta[i] = vv[i] = ww[i] = uu[i] = 0; 269 } 270 k = 0; 271 s = p; 272 r = 0; 273 while (endless) { 274 t = s.next; 275 if (k === 0) { 276 // MP 306 277 if (s.rtype === this.MP_GIVEN) { 278 // MP 314 279 if (t.ltype === this.MP_GIVEN) { 280 aa = Math.atan2(delta_y[0], delta_x[0]); 281 ct_st = this.mp_n_sin_cos(p.right_given() - aa); 282 ct = ct_st[0]; 283 st = ct_st[1]; 284 cf_sf = this.mp_n_sin_cos(q.left_given() - aa); 285 cf = cf_sf[0]; 286 sf = cf_sf[1]; 287 this.mp_set_controls(p, q, delta_x[0], delta_y[0], st, ct, -sf, cf); 288 return; 289 } 290 vv[0] = s.right_given() - Math.atan2(delta_y[0], delta_x[0]); 291 vv[0] = this.reduce_angle(vv[0]); 292 uu[0] = 0; 293 ww[0] = 0; 294 } else if (s.rtype === this.MP_CURL) { 295 // MP 315 296 if (t.ltype === this.MP_CURL) { 297 p.rtype = this.MP_EXPLICIT; 298 q.ltype = this.MP_EXPLICIT; 299 lt = Math.abs(q.left_tension()); 300 rt = Math.abs(p.right_tension()); 301 ff = this.UNITY / (3.0 * rt); 302 p.rx = p.x + delta_x[0] * ff; 303 p.ry = p.y + delta_y[0] * ff; 304 ff = this.UNITY / (3.0 * lt); 305 q.lx = q.x - delta_x[0] * ff; 306 q.ly = q.y - delta_y[0] * ff; 307 return; 308 } 309 cc = s.right_curl(); 310 lt = Math.abs(t.left_tension()); 311 rt = Math.abs(s.right_tension()); 312 uu[0] = this.mp_curl_ratio(cc, rt, lt); 313 vv[0] = -psi[1] * uu[0]; 314 ww[0] = 0; 315 } else { 316 if (s.rtype === this.MP_OPEN) { 317 uu[0] = 0; 318 vv[0] = 0; 319 ww[0] = this.FRACTION_ONE; 320 } 321 } 322 } else { 323 if (s.ltype === this.MP_END_CYCLE || s.ltype === this.MP_OPEN) { 324 // MP 308 325 aa = this.UNITY / (3.0 * Math.abs(r.right_tension()) - this.UNITY); 326 dd = 327 delta[k] * 328 (this.FRACTION_THREE - this.UNITY / Math.abs(r.right_tension())); 329 bb = this.UNITY / (3 * Math.abs(t.left_tension()) - this.UNITY); 330 ee = 331 delta[k - 1] * 332 (this.FRACTION_THREE - this.UNITY / Math.abs(t.left_tension())); 333 cc = this.FRACTION_ONE - uu[k - 1] * aa; 334 dd = dd * cc; 335 lt = Math.abs(s.left_tension()); 336 rt = Math.abs(s.right_tension()); 337 if (lt < rt) { 338 dd *= Math.pow(lt / rt, 2); 339 } else { 340 if (lt > rt) { 341 ee *= Math.pow(rt / lt, 2); 342 } 343 } 344 ff = ee / (ee + dd); 345 uu[k] = ff * bb; 346 acc = -psi[k + 1] * uu[k]; 347 if (r.rtype === this.MP_CURL) { 348 ww[k] = 0; 349 vv[k] = acc - psi[1] * (this.FRACTION_ONE - ff); 350 } else { 351 ff = (this.FRACTION_ONE - ff) / cc; 352 acc = acc - psi[k] * ff; 353 ff = ff * aa; 354 vv[k] = acc - vv[k - 1] * ff; 355 ww[k] = -ww[k - 1] * ff; 356 } 357 if (s.ltype === this.MP_END_CYCLE) { 358 aa = 0; 359 bb = this.FRACTION_ONE; 360 while (endless) { 361 k -= 1; 362 if (k === 0) { 363 k = n; 364 } 365 aa = vv[k] - aa * uu[k]; 366 bb = ww[k] - bb * uu[k]; 367 if (k === n) { 368 break; 369 } 370 } 371 aa = aa / (this.FRACTION_ONE - bb); 372 theta[n] = aa; 373 vv[0] = aa; 374 // k_val = range(1, n); 375 // for (k_idx in k_val) { 376 // k = k_val[k_idx]; 377 for (k_idx = 1; k_idx < n; k_idx++) { 378 vv[k_idx] = vv[k_idx] + aa * ww[k_idx]; 379 } 380 break; 381 } 382 } else { 383 if (s.ltype === this.MP_CURL) { 384 cc = s.left_curl(); 385 lt = Math.abs(s.left_tension()); 386 rt = Math.abs(r.right_tension()); 387 ff = this.mp_curl_ratio(cc, lt, rt); 388 theta[n] = -(vv[n - 1] * ff) / (this.FRACTION_ONE - ff * uu[n - 1]); 389 break; 390 } 391 if (s.ltype === this.MP_GIVEN) { 392 theta[n] = s.left_given() - Math.atan2(delta_y[n - 1], delta_x[n - 1]); 393 theta[n] = this.reduce_angle(theta[n]); 394 break; 395 } 396 } 397 } 398 r = s; 399 s = t; 400 k += 1; 401 } 402 403 // MP 318 404 for (k = n - 1; k > -1; k--) { 405 theta[k] = vv[k] - theta[k + 1] * uu[k]; 406 } 407 408 s = p; 409 k = 0; 410 while (endless) { 411 t = s.next; 412 ct_st = this.mp_n_sin_cos(theta[k]); 413 ct = ct_st[0]; 414 st = ct_st[1]; 415 cf_sf = this.mp_n_sin_cos(-psi[k + 1] - theta[k + 1]); 416 cf = cf_sf[0]; 417 sf = cf_sf[1]; 418 this.mp_set_controls(s, t, delta_x[k], delta_y[k], st, ct, sf, cf); 419 k++; 420 s = t; 421 if (k === n) { 422 break; 423 } 424 } 425 }, 426 427 /** 428 * @private 429 */ 430 mp_n_sin_cos: function (z) { 431 return [Math.cos(z), Math.sin(z)]; 432 }, 433 434 /** 435 * @private 436 */ 437 mp_set_controls: function (p, q, delta_x, delta_y, st, ct, sf, cf) { 438 var rt, ss, lt, sine, rr; 439 lt = Math.abs(q.left_tension()); 440 rt = Math.abs(p.right_tension()); 441 rr = this.mp_velocity(st, ct, sf, cf, rt); 442 ss = this.mp_velocity(sf, cf, st, ct, lt); 443 444 // console.log('lt rt rr ss', lt, rt, rr, ss); 445 if (p.right_tension() < 0 || q.left_tension() < 0) { 446 if ((st >= 0 && sf >= 0) || (st <= 0 && sf <= 0)) { 447 sine = Math.abs(st) * cf + Math.abs(sf) * ct; 448 if (sine > 0) { 449 sine *= 1.00024414062; 450 if (p.right_tension() < 0) { 451 if (this.mp_ab_vs_cd(Math.abs(sf), this.FRACTION_ONE, rr, sine) < 0) { 452 rr = Math.abs(sf) / sine; 453 } 454 } 455 if (q.left_tension() < 0) { 456 if (this.mp_ab_vs_cd(Math.abs(st), this.FRACTION_ONE, ss, sine) < 0) { 457 ss = Math.abs(st) / sine; 458 } 459 } 460 } 461 } 462 } 463 p.rx = p.x + (delta_x * ct - delta_y * st) * rr; 464 p.ry = p.y + (delta_y * ct + delta_x * st) * rr; 465 q.lx = q.x - (delta_x * cf + delta_y * sf) * ss; 466 q.ly = q.y - (delta_y * cf - delta_x * sf) * ss; 467 p.rtype = this.MP_EXPLICIT; 468 q.ltype = this.MP_EXPLICIT; 469 }, 470 471 /** 472 * @private 473 */ 474 mp_pyth_add: function (a, b) { 475 return Math.sqrt(a * a + b * b); 476 }, 477 478 /** 479 * 480 * @private 481 */ 482 mp_curl_ratio: function (gamma, a_tension, b_tension) { 483 var alpha = 1.0 / a_tension, 484 beta = 1.0 / b_tension; 485 486 return Math.min( 487 4.0, 488 ((3.0 - alpha) * alpha * alpha * gamma + beta * beta * beta) / 489 (alpha * alpha * alpha * gamma + (3.0 - beta) * beta * beta) 490 ); 491 }, 492 493 /** 494 * @private 495 */ 496 mp_ab_vs_cd: function (a, b, c, d) { 497 if (a * b === c * d) { 498 return 0; 499 } 500 if (a * b > c * d) { 501 return 1; 502 } 503 return -1; 504 }, 505 506 /** 507 * @private 508 */ 509 mp_velocity: function (st, ct, sf, cf, t) { 510 return Math.min( 511 4.0, 512 (2.0 + Math.sqrt(2) * (st - sf / 16.0) * (sf - st / 16.0) * (ct - cf)) / 513 (1.5 * t * (2 + (Math.sqrt(5) - 1) * ct + (3 - Math.sqrt(5)) * cf)) 514 ); 515 }, 516 517 /** 518 * @private 519 * @param {Number} A 520 */ 521 reduce_angle: function (A) { 522 if (Math.abs(A) > this.ONE_EIGHTY_DEG) { 523 if (A > 0) { 524 A -= this.THREE_SIXTY_DEG; 525 } else { 526 A += this.THREE_SIXTY_DEG; 527 } 528 } 529 return A; 530 }, 531 532 /** 533 * 534 * @private 535 * @param {Array} p 536 * @param {Number} tension 537 * @param {Boolean} cycle 538 */ 539 makeknots: function (p, tension, cycle) { 540 var i, 541 len, 542 knots = []; 543 544 tension = tension || 1; 545 546 len = p.length; 547 for (i = 0; i < len; i++) { 548 knots.push({ 549 x: p[i][0], 550 y: p[i][1], 551 ltype: this.MP_OPEN, 552 rtype: this.MP_OPEN, 553 ly: tension, 554 ry: tension, 555 lx: tension, 556 rx: tension, 557 left_curl: function () { 558 return this.lx || 0; 559 }, 560 right_curl: function () { 561 return this.rx || 0; 562 }, 563 left_tension: function () { 564 if (!this.ly) { 565 this.ly = 1; 566 } 567 return this.ly; 568 }, 569 right_tension: function () { 570 if (!this.ry) { 571 this.ry = 1; 572 } 573 return this.ry; 574 }, 575 set_right_curl: function (x) { 576 this.rx = x || 0; 577 }, 578 set_left_curl: function (x) { 579 this.lx = x || 0; 580 } 581 }); 582 } 583 len = knots.length; 584 for (i = 0; i < len; i++) { 585 knots[i].next = knots[i + 1] || knots[i]; 586 knots[i].set_right_given = knots[i].set_right_curl; 587 knots[i].set_left_given = knots[i].set_left_curl; 588 knots[i].right_given = knots[i].right_curl; 589 knots[i].left_given = knots[i].left_curl; 590 } 591 knots[len - 1].next = knots[0]; 592 593 if (!cycle) { 594 knots[len - 1].rtype = this.MP_ENDPOINT; 595 596 knots[len - 1].ltype = this.MP_CURL; 597 knots[0].rtype = this.MP_CURL; 598 } 599 600 return knots; 601 }, 602 603 /** 604 * 605 * @param {Array} point_list 606 * @param {Object} controls 607 * 608 * @returns {Array} 609 */ 610 curve: function (point_list, controls) { 611 var knots, 612 len, 613 i, 614 val, 615 x = [], 616 y = []; 617 618 controls = controls || { 619 tension: 1, 620 direction: {}, 621 curl: {}, 622 isClosed: false 623 }; 624 625 knots = this.makeknots(point_list, Type.evaluate(controls.tension), controls.isClosed); 626 627 len = knots.length; 628 for (i in controls.direction) { 629 if (controls.direction.hasOwnProperty(i)) { 630 val = Type.evaluate(controls.direction[i]); 631 if (Type.isArray(val)) { 632 if (val[0] !== false) { 633 knots[i].lx = (val[0] * Math.PI) / 180; 634 knots[i].ltype = this.MP_GIVEN; 635 } 636 if (val[1] !== false) { 637 knots[i].rx = (val[1] * Math.PI) / 180; 638 knots[i].rtype = this.MP_GIVEN; 639 } 640 } else { 641 knots[i].lx = (val * Math.PI) / 180; 642 knots[i].rx = (val * Math.PI) / 180; 643 knots[i].ltype = knots[i].rtype = this.MP_GIVEN; 644 } 645 } 646 } 647 for (i in controls.curl) { 648 if (controls.curl.hasOwnProperty(i)) { 649 val = Type.evaluate(controls.curl[i]); 650 if (parseInt(i, 10) === 0) { 651 knots[i].rtype = this.MP_CURL; 652 knots[i].set_right_curl(val); 653 } else if (parseInt(i, 10) === len - 1) { 654 knots[i].ltype = this.MP_CURL; 655 knots[i].set_left_curl(val); 656 } 657 } 658 } 659 660 this.make_choices(knots); 661 662 for (i = 0; i < len - 1; i++) { 663 x.push(knots[i].x); 664 x.push(knots[i].rx); 665 x.push(knots[i + 1].lx); 666 y.push(knots[i].y); 667 y.push(knots[i].ry); 668 y.push(knots[i + 1].ly); 669 } 670 x.push(knots[len - 1].x); 671 y.push(knots[len - 1].y); 672 673 if (controls.isClosed) { 674 x.push(knots[len - 1].rx); 675 y.push(knots[len - 1].ry); 676 x.push(knots[0].lx); 677 y.push(knots[0].ly); 678 x.push(knots[0].x); 679 y.push(knots[0].y); 680 } 681 682 return [x, y]; 683 } 684 }; 685 686 export default Mat.Metapost; 687