1 /*
  2     Copyright 2008-2021
  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