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