1 /* 2 Copyright 2008-2022 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 <http://www.gnu.org/licenses/> 29 and <http://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 (http://pyx.sourceforge.net/) 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 (http://tug.org/metapost). 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 http://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 /* depends: 63 utils/type 64 math/math 65 */ 66 67 /** 68 * @fileoverview In this file the namespace Math.Metapost is defined which holds algorithms translated from Metapost 69 * by D.E. Knuth and J.D. Hobby. 70 */ 71 72 define(['utils/type', 'math/math'], function (Type, Mat) { 73 74 "use strict"; 75 76 /** 77 * The JXG.Math.Metapost namespace holds algorithms translated from Metapost 78 * by D.E. Knuth and J.D. Hobby. 79 * 80 * @name JXG.Math.Metapost 81 * @exports Mat.Metapost as JXG.Math.Metapost 82 * @namespace 83 */ 84 Mat.Metapost = { 85 MP_ENDPOINT: 0, 86 MP_EXPLICIT: 1, 87 MP_GIVEN: 2, 88 MP_CURL: 3, 89 MP_OPEN: 4, 90 MP_END_CYCLE: 5, 91 92 UNITY: 1.0, 93 // two: 2, 94 // fraction_half: 0.5, 95 FRACTION_ONE: 1.0, 96 FRACTION_THREE: 3.0, 97 ONE_EIGHTY_DEG: Math.PI, 98 THREE_SIXTY_DEG: 2 * Math.PI, 99 // EPSILON: 1e-5, 100 EPS_SQ: 1e-5 * 1e-5, 101 102 /** 103 * @private 104 */ 105 make_choices: function (knots) { 106 var dely, h, k, delx, n, 107 q, p, s, cosine, t, sine, 108 delta_x, delta_y, delta, psi; 109 110 p = knots[0]; 111 do { 112 if (!p) { 113 break; 114 } 115 q = p.next; 116 117 // Join two identical knots by setting the control points to the same 118 // coordinates. 119 // MP 291 120 if (p.rtype > this.MP_EXPLICIT && 121 ((p.x - q.x) * (p.x - q.x) + (p.y - q.y) * (p.y - q.y) < this.EPS_SQ)) { 122 123 p.rtype = this.MP_EXPLICIT; 124 if (p.ltype === this.MP_OPEN) { 125 p.ltype = this.MP_CURL; 126 p.set_left_curl(this.UNITY); 127 } 128 129 q.ltype = this.MP_EXPLICIT; 130 if (q.rtype === this.MP_OPEN) { 131 q.rtype = this.MP_CURL; 132 q.set_right_curl(this.UNITY); 133 } 134 135 p.rx = p.x; 136 q.lx = p.x; 137 p.ry = p.y; 138 q.ly = p.y; 139 } 140 p = q; 141 } while (p !== knots[0]); 142 143 // Find the first breakpoint, h, on the path 144 // MP 292 145 h = knots[0]; 146 while (true) { 147 if (h.ltype !== this.MP_OPEN || h.rtype !== this.MP_OPEN) { 148 break; 149 } 150 h = h.next; 151 if (h === knots[0]) { 152 h.ltype = this.MP_END_CYCLE; 153 break; 154 } 155 } 156 157 p = h; 158 while (true) { 159 if (!p) { 160 break; 161 } 162 163 // Fill in the control points between p and the next breakpoint, 164 // then advance p to that breakpoint 165 // MP 299 166 q = p.next; 167 if (p.rtype >= this.MP_GIVEN) { 168 while (q.ltype === this.MP_OPEN && q.rtype === this.MP_OPEN) { 169 q = q.next; 170 } 171 172 // Calculate the turning angles psi_ k and the distances d_{k,k+1}; 173 // set n to the length of the path 174 // MP 302 175 k = 0; 176 s = p; 177 n = knots.length; 178 179 delta_x = []; 180 delta_y = []; 181 delta = []; 182 psi = [null]; 183 184 // tuple([]) = tuple([[], [], [], [null]]); 185 while (true) { 186 t = s.next; 187 // None; 188 delta_x.push(t.x - s.x); 189 delta_y.push(t.y - s.y); 190 delta.push( this.mp_pyth_add(delta_x[k], delta_y[k]) ); 191 if (k > 0) { 192 sine = delta_y[k - 1] / delta[k - 1]; 193 cosine = delta_x[k - 1] / delta[k - 1]; 194 psi.push( 195 Math.atan2( 196 delta_y[k] * cosine - delta_x[k] * sine, 197 delta_x[k] * cosine + delta_y[k] * sine 198 ) 199 ); 200 } 201 k++; 202 s = t; 203 if (s === q) { 204 n = k; 205 } 206 if (k >= n && s.ltype !== this.MP_END_CYCLE) { 207 break; 208 } 209 } 210 if (k === n) { 211 psi.push(0); 212 } else { 213 psi.push(psi[1]); 214 } 215 216 // Remove open types at the breakpoints 217 // MP 303 218 if (q.ltype === this.MP_OPEN) { 219 delx = (q.rx - q.x); 220 dely = (q.ry - q.y); 221 if (delx * delx + dely * dely < this.EPS_SQ) { 222 q.ltype = this.MP_CURL; 223 q.set_left_curl(this.UNITY); 224 } else { 225 q.ltype = this.MP_GIVEN; 226 q.set_left_given(Math.atan2(dely, delx)); 227 } 228 } 229 if (p.rtype === this.MP_OPEN && p.ltype === this.MP_EXPLICIT) { 230 delx = (p.x - p.lx); 231 dely = (p.y - p.ly); 232 if ( delx * delx + dely * dely < this.EPS_SQ) { 233 p.rtype = this.MP_CURL; 234 p.set_right_curl(this.UNITY); 235 } else { 236 p.rtype = this.MP_GIVEN; 237 p.set_right_given(Math.atan2(dely, delx)); 238 } 239 } 240 this.mp_solve_choices(p, q, n, delta_x, delta_y, delta, psi); 241 } else if (p.rtype === this.MP_ENDPOINT) { 242 // MP 294 243 p.rx = p.x; 244 p.ry = p.y; 245 q.lx = q.x; 246 q.ly = q.y; 247 } 248 p = q; 249 250 if (p === h) { 251 break; 252 } 253 } 254 }, 255 256 /** 257 * Implements solve_choices form metapost 258 * MP 305 259 * @private 260 */ 261 mp_solve_choices: function (p, q, n, delta_x, delta_y, delta, psi) { 262 var aa, acc, vv, bb, ldelta, ee, k, s, 263 ww, uu, lt, r, t, ff, theta, rt, dd, cc, 264 ct_st, ct, st, cf_sf, cf, sf, i, 265 k_idx; 266 267 ldelta = delta.length + 1; 268 uu = new Array(ldelta); 269 ww = new Array(ldelta); 270 vv = new Array(ldelta); 271 theta = new Array(ldelta); 272 for (i = 0; i < ldelta; i++) { 273 theta[i] = vv[i] = ww[i] = uu[i] = 0; 274 } 275 k = 0; 276 s = p; 277 r = 0; 278 while (true) { 279 t = s.next; 280 if (k === 0) { 281 // MP 306 282 if (s.rtype === this.MP_GIVEN) { 283 // MP 314 284 if (t.ltype === this.MP_GIVEN) { 285 aa = Math.atan2(delta_y[0], delta_x[0]); 286 ct_st = this.mp_n_sin_cos(p.right_given() - aa); 287 ct = ct_st[0]; 288 st = ct_st[1]; 289 cf_sf = this.mp_n_sin_cos(q.left_given() - aa); 290 cf = cf_sf[0]; 291 sf = cf_sf[1]; 292 this.mp_set_controls(p, q, delta_x[0], delta_y[0], st, ct, -sf, cf); 293 return; 294 } 295 vv[0] = s.right_given() - Math.atan2(delta_y[0], delta_x[0]); 296 vv[0] = this.reduce_angle(vv[0]); 297 uu[0] = 0; 298 ww[0] = 0; 299 300 } else if (s.rtype === this.MP_CURL) { 301 // MP 315 302 if (t.ltype === this.MP_CURL) { 303 p.rtype = this.MP_EXPLICIT; 304 q.ltype = this.MP_EXPLICIT; 305 lt = Math.abs(q.left_tension()); 306 rt = Math.abs(p.right_tension()); 307 ff = this.UNITY / (3.0 * rt); 308 p.rx = p.x + delta_x[0] * ff; 309 p.ry = p.y + delta_y[0] * ff; 310 ff = this.UNITY / (3.0 * lt); 311 q.lx = q.x - delta_x[0] * ff; 312 q.ly = q.y - delta_y[0] * ff; 313 return; 314 } 315 cc = s.right_curl(); 316 lt = Math.abs(t.left_tension()); 317 rt = Math.abs(s.right_tension()); 318 uu[0] = this.mp_curl_ratio(cc, rt, lt); 319 vv[0] = -psi[1] * uu[0]; 320 ww[0] = 0; 321 } else { 322 if (s.rtype === this.MP_OPEN) { 323 uu[0] = 0; 324 vv[0] = 0; 325 ww[0] = this.FRACTION_ONE; 326 } 327 } 328 } else { 329 if (s.ltype === this.MP_END_CYCLE || s.ltype === this.MP_OPEN) { 330 // MP 308 331 aa = this.UNITY / (3.0 * Math.abs(r.right_tension()) - this.UNITY); 332 dd = delta[k] * (this.FRACTION_THREE - this.UNITY / Math.abs(r.right_tension())); 333 bb = this.UNITY / (3 * Math.abs(t.left_tension()) - this.UNITY); 334 ee = delta[k - 1] * (this.FRACTION_THREE - this.UNITY / Math.abs(t.left_tension())); 335 cc = this.FRACTION_ONE - uu[k - 1] * aa; 336 dd = dd * cc; 337 lt = Math.abs(s.left_tension()); 338 rt = Math.abs(s.right_tension()); 339 if (lt < rt) { 340 dd *= Math.pow(lt / rt, 2); 341 } else { 342 if (lt > rt) { 343 ee *= Math.pow(rt / lt, 2); 344 } 345 } 346 ff = ee / (ee + dd); 347 uu[k] = ff * bb; 348 acc = -psi[k + 1] * uu[k]; 349 if (r.rtype === this.MP_CURL) { 350 ww[k] = 0; 351 vv[k] = acc - psi[1] * (this.FRACTION_ONE - ff); 352 } else { 353 ff = (this.FRACTION_ONE - ff) / cc; 354 acc = acc - psi[k] * ff; 355 ff = ff * aa; 356 vv[k] = acc - vv[k - 1] * ff; 357 ww[k] = -ww[k - 1] * ff; 358 } 359 if (s.ltype === this.MP_END_CYCLE) { 360 aa = 0; 361 bb = this.FRACTION_ONE; 362 while (true) { 363 k -= 1; 364 if (k === 0) { 365 k = n; 366 } 367 aa = vv[k] - aa * uu[k]; 368 bb = ww[k] - bb * uu[k]; 369 if (k === n) { 370 break; 371 } 372 } 373 aa = aa / (this.FRACTION_ONE - bb); 374 theta[n] = aa; 375 vv[0] = aa; 376 // k_val = range(1, n); 377 // for (k_idx in k_val) { 378 // k = k_val[k_idx]; 379 for (k_idx = 1; k_idx < n; k_idx++) { 380 vv[k_idx] = vv[k_idx] + aa * ww[k_idx]; 381 } 382 break; 383 } 384 } else { 385 if (s.ltype === this.MP_CURL) { 386 cc = s.left_curl(); 387 lt = Math.abs(s.left_tension()); 388 rt = Math.abs(r.right_tension()); 389 ff = this.mp_curl_ratio(cc, lt, rt); 390 theta[n] = -(vv[n - 1] * ff) / (this.FRACTION_ONE - ff * uu[n - 1]); 391 break; 392 } 393 if (s.ltype === this.MP_GIVEN) { 394 theta[n] = s.left_given() - Math.atan2(delta_y[n - 1], delta_x[n - 1]); 395 theta[n] = this.reduce_angle(theta[n]); 396 break; 397 } 398 } 399 } 400 r = s; 401 s = t; 402 k += 1; 403 } 404 405 // MP 318 406 for (k = n-1; k > -1; k--) { 407 theta[k] = vv[k] - theta[k + 1] * uu[k]; 408 } 409 410 s = p; 411 k = 0; 412 while (true) { 413 t = s.next; 414 ct_st = this.mp_n_sin_cos(theta[k]); 415 ct = ct_st[0]; 416 st = ct_st[1]; 417 cf_sf = this.mp_n_sin_cos((-(psi[k + 1]) - theta[k + 1])); 418 cf = cf_sf[0]; 419 sf = cf_sf[1]; 420 this.mp_set_controls(s, t, delta_x[k], delta_y[k], st, ct, sf, cf); 421 k++; 422 s = t; 423 if (k === n) { 424 break; 425 } 426 } 427 }, 428 429 /** 430 * @private 431 */ 432 mp_n_sin_cos: function (z) { 433 return [Math.cos(z), Math.sin(z)]; 434 }, 435 436 /** 437 * @private 438 */ 439 mp_set_controls: function (p, q, delta_x, delta_y, st, ct, sf, cf) { 440 var rt, ss, lt, sine, rr; 441 lt = Math.abs(q.left_tension()); 442 rt = Math.abs(p.right_tension()); 443 rr = this.mp_velocity(st, ct, sf, cf, rt); 444 ss = this.mp_velocity(sf, cf, st, ct, lt); 445 446 // console.log('lt rt rr ss', lt, rt, rr, ss); 447 if (p.right_tension() < 0 || q.left_tension() < 0) { 448 if ((st >= 0 && sf >= 0) || (st <= 0 && sf <= 0)) { 449 sine = Math.abs(st) * cf + Math.abs(sf) * ct; 450 if (sine > 0) { 451 sine *= 1.00024414062; 452 if (p.right_tension() < 0) { 453 if (this.mp_ab_vs_cd(Math.abs(sf), this.FRACTION_ONE, rr, sine) < 0) { 454 rr = Math.abs(sf) / sine; 455 } 456 } 457 if (q.left_tension() < 0) { 458 if (this.mp_ab_vs_cd(Math.abs(st), this.FRACTION_ONE, ss, sine) < 0) { 459 ss = Math.abs(st) / sine; 460 } 461 } 462 } 463 } 464 } 465 p.rx = p.x + (delta_x * ct - delta_y * st) * rr; 466 p.ry = p.y + (delta_y * ct + delta_x * st) * rr; 467 q.lx = q.x - (delta_x * cf + delta_y * sf) * ss; 468 q.ly = q.y - (delta_y * cf - delta_x * sf) * ss; 469 p.rtype = this.MP_EXPLICIT; 470 q.ltype = this.MP_EXPLICIT; 471 }, 472 473 /** 474 * @private 475 */ 476 mp_pyth_add: function (a, b) { 477 return Math.sqrt((a * a + b * b)); 478 }, 479 480 /** 481 * 482 * @private 483 */ 484 mp_curl_ratio: function (gamma, a_tension, b_tension) { 485 var alpha = 1.0 / a_tension, 486 beta = 1.0 / b_tension; 487 488 return Math.min (4.0, 489 ((3.0 - alpha) * alpha * alpha * gamma + beta * beta * beta) / 490 (alpha * alpha * alpha * gamma + (3.0 - beta) * beta * beta) 491 ); 492 }, 493 494 /** 495 * @private 496 */ 497 mp_ab_vs_cd: function (a, b, c, d) { 498 if (a * b === c * d) { 499 return 0; 500 } 501 if (a * b > c * d) { 502 return 1; 503 } 504 return -1; 505 }, 506 507 /** 508 * @private 509 */ 510 mp_velocity: function (st, ct, sf, cf, t) { 511 return Math.min (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, len, 541 knots = []; 542 543 tension = tension || 1; 544 545 len = p.length; 546 for (i = 0; i < len; i++) { 547 knots.push({ 548 x: p[i][0], 549 y: p[i][1], 550 ltype: this.MP_OPEN, 551 rtype: this.MP_OPEN, 552 ly: tension, 553 ry: tension, 554 lx: tension, 555 rx: tension, 556 left_curl: function() { return this.lx || 0; }, 557 right_curl: function() { return this.rx || 0; }, 558 left_tension: function() { 559 if (!this.ly) { this.ly = 1; } 560 return this.ly; 561 }, 562 right_tension: function() { 563 if (!this.ry) { this.ry = 1; } 564 return this.ry; 565 }, 566 set_right_curl: function(x) { this.rx = x || 0; }, 567 set_left_curl: function(x) { this.lx = x || 0; } 568 }); 569 } 570 len = knots.length; 571 for (i = 0; i < len; i++) { 572 knots[i].next = knots[i+1] || knots[i]; 573 knots[i].set_right_given = knots[i].set_right_curl; 574 knots[i].set_left_given = knots[i].set_left_curl; 575 knots[i].right_given = knots[i].right_curl; 576 knots[i].left_given = knots[i].left_curl; 577 } 578 knots[len - 1].next = knots[0]; 579 580 if (!cycle) { 581 knots[len - 1].rtype = this.MP_ENDPOINT; 582 583 knots[len - 1].ltype = this.MP_CURL; 584 knots[0].rtype = this.MP_CURL; 585 } 586 587 return knots; 588 }, 589 590 /** 591 * 592 * @param {Array} point_list 593 * @param {Object} controls 594 * 595 * @returns {Array} 596 */ 597 curve: function(point_list, controls) { 598 var knots, len, i, val, 599 x = [], 600 y = []; 601 602 controls = controls || { 603 tension: 1, 604 direction: {}, 605 curl: {}, 606 isClosed: false 607 }; 608 609 knots = this.makeknots(point_list, Type.evaluate(controls.tension), controls.isClosed); 610 611 len = knots.length; 612 for (i in controls.direction) { 613 if (controls.direction.hasOwnProperty(i)) { 614 val = Type.evaluate(controls.direction[i]); 615 if (Type.isArray(val)) { 616 if (val[0] !== false) { 617 knots[i].lx = val[0] * Math.PI / 180; 618 knots[i].ltype = this.MP_GIVEN; 619 } 620 if (val[1] !== false) { 621 knots[i].rx = val[1] * Math.PI / 180; 622 knots[i].rtype = this.MP_GIVEN; 623 } 624 } else { 625 knots[i].lx = val * Math.PI / 180; 626 knots[i].rx = val * Math.PI / 180; 627 knots[i].ltype = knots[i].rtype = this.MP_GIVEN; 628 } 629 } 630 } 631 for (i in controls.curl) { 632 if (controls.curl.hasOwnProperty(i)) { 633 val = Type.evaluate(controls.curl[i]); 634 if (parseInt(i, 10) === 0) { 635 knots[i].rtype = this.MP_CURL; 636 knots[i].set_right_curl(val); 637 } else if (parseInt(i, 10) === len - 1) { 638 knots[i].ltype = this.MP_CURL; 639 knots[i].set_left_curl(val); 640 } 641 } 642 } 643 644 this.make_choices(knots); 645 646 for (i = 0; i < len - 1; i++) { 647 x.push(knots[i].x); 648 x.push(knots[i].rx); 649 x.push(knots[i + 1].lx); 650 y.push(knots[i].y); 651 y.push(knots[i].ry); 652 y.push(knots[i + 1].ly); 653 } 654 x.push(knots[len - 1].x); 655 y.push(knots[len - 1].y); 656 657 if (controls.isClosed) { 658 x.push(knots[len - 1].rx); 659 y.push(knots[len - 1].ry); 660 x.push(knots[0].lx); 661 y.push(knots[0].ly); 662 x.push(knots[0].x); 663 y.push(knots[0].y); 664 } 665 666 return [x, y]; 667 } 668 669 }; 670 671 return Mat.Metapost; 672 }); 673