pyx/metapost/mp_path.py

   1 # -*- coding: ISO-8859-1 -*-
   2 #
   3 # Copyright (C) 2011 Michael Schindler <m-schindler@users.sourceforge.net>
   4 #
   5 # This file is part of PyX (http://pyx.sourceforge.net/).
   6 #
   7 # PyX is free software; you can redistribute it and/or modify
   8 # it under the terms of the GNU General Public License as published by
   9 # the Free Software Foundation; either version 2 of the License, or
  10 # (at your option) any later version.
  11 #
  12 # PyX is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU General Public License
  18 # along with PyX; if not, write to the Free Software
  19 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  20
  21 from math import atan2, sin, cos, sqrt, pi
  22
  23 ################################################################################
  24 # Internal functions of MetaPost
  25 #
  26 # This file re-implements some of the functionality of MetaPost
  27 # (http://tug.org/metapost). MetaPost was developed by John D. Hobby
  28 # and others.
  29 #
  30 # This file is based on the MetaPost version distributed by TeXLive:
  31 # svn://tug.org/texlive/trunk/Build/source/texk/web2c/mplibdir
  32 # revision 22737 (2011-05-31)
  33 ################################################################################
  34
  35 # from mplib.h:
  36 mp_endpoint = 0
  37 mp_explicit = 1
  38 mp_given = 2
  39 mp_curl = 3
  40 mp_open = 4
  41 mp_end_cycle = 5
  42
  43 # from mpmath.c:
  44 unity = 1.0
  45 two = 2.0
  46 fraction_half = 0.5
  47 fraction_one = 1.0
  48 fraction_three = 3.0
  49 one_eighty_deg = pi
  50 three_sixty_deg = 2*pi
  51
  52 def mp_make_choices(knots): # <<<
  53     """Implements mp_make_choices from metapost (mp.c)"""
  54     # 334: If consecutive knots are equal, join them explicitly
  55     p = knots
  56     while True:
  57         q = p.next
  58         # XXX float comparison: use epsilon
  59         if p.x_pt == q.x_pt and p.y_pt == q.y_pt and p.rtype > mp_explicit:
  60             p.rtype = mp_explicit
  61             if p.ltype == mp_open:
  62                 p.ltype = mp_curl
  63                 p.set_left_curl(unity)
  64             q.ltype = mp_explicit
  65             if q.rtype == mp_open:
  66                 q.rtype = mp_curl
  67                 q.set_right_curl(unity)
  68             p.rx_pt = p.x_pt
  69             q.lx_pt = p.x_pt
  70             p.ry_pt = p.y_pt
  71             q.ly_pt = p.y_pt
  72         p = q
  73         if p is knots:
  74             break
  75
  76     # 335:
  77     # If there are no breakpoints, it is necessary to compute the direction angles around an entire cycle.
  78     # In this case the mp left type of the first node is temporarily changed to end cycle.
  79     # Find the first breakpoint, h, on the path
  80     # insert an artificial breakpoint if the path is an unbroken cycle
  81     h = knots
  82     while True:
  83         if h.ltype != mp_open or h.rtype != mp_open:
  84             break
  85         h = h.next
  86         if h is knots:
  87             h.ltype = mp_end_cycle
  88             break
  89
  90     p = h
  91     while True:
  92         # 336:
  93         # Fill in the control points between p and the next breakpoint, then advance p to that breakpoint
  94         q = p.next
  95         if p.rtype >= mp_given:
  96             while q.ltype == mp_open and q.rtype == mp_open:
  97                 q = q.next
  98
  99             # 346:
 100             # Calculate the turning angles psi_k and the distances d(k, k+1)
 101             # set n to the length of the path
 102             k = 0
 103             s = p
 104             n = knots.linked_len()
 105             delta_x, delta_y, delta, psi = [], [], [], [None]
 106             while True:
 107                 t = s.next
 108                 assert len(delta_x) == k
 109                 delta_x.append(t.x_pt - s.x_pt)
 110                 delta_y.append(t.y_pt - s.y_pt)
 111                 delta.append(mp_pyth_add(delta_x[k], delta_y[k]))
 112                 if k > 0:
 113                     sine = mp_make_fraction(delta_y[k-1], delta[k-1])
 114                     cosine = mp_make_fraction(delta_x[k-1], delta[k-1])
 115                     psi.append(mp_n_arg(
 116                       mp_take_fraction(delta_x[k], cosine) + mp_take_fraction(delta_y[k], sine),
 117                       mp_take_fraction(delta_y[k], cosine) - mp_take_fraction(delta_x[k], sine)))
 118                 k += 1
 119                 s = t
 120                 if s == q:
 121                     n = k
 122                 if k >= n and s.ltype != mp_end_cycle:
 123                     break
 124             if k == n:
 125                 psi.append(0)
 126             else:
 127                 # for closed paths:
 128                 psi.append(psi[1])
 129
 130             # 347: Remove open types at the breakpoints
 131             if q.ltype == mp_open:
 132                 delx_pt = q.rx_pt - q.x_pt
 133                 dely_pt = q.ry_pt - q.y_pt
 134                 if delx_pt == 0 and dely_pt == 0: # XXX float comparison
 135                     q.ltype = mp_curl
 136                     q.set_left_curl(unity)
 137                 else:
 138                     q.ltype = mp_given
 139                     q.set_left_given(mp_n_arg(delx_pt, dely_pt))
 140
 141             if p.rtype == mp_open and p.ltype == mp_explicit:
 142                 delx_pt = p.x_pt - p.lx_pt
 143                 dely_pt = p.y_pt - p.ly_pt
 144                 if delx_pt == 0 and dely_pt == 0: # XXX float comparison
 145                     p.rtype = mp_curl
 146                     p.set_right_curl(unity)
 147                 else:
 148                     p.rtype = mp_given
 149                     p.set_right_given(mp_n_arg(delx_pt, dely_pt))
 150
 151             # call the internal solving routine
 152             mp_solve_choices(p, q, n, delta_x, delta_y, delta, psi)
 153
 154         elif p.rtype == mp_endpoint:
 155             # 337: Give reasonable values for the unused control points between p and q
 156             p.rx_pt = p.x_pt
 157             p.ry_pt = p.y_pt
 158             q.lx_pt = q.x_pt
 159             q.ly_pt = q.y_pt
 160
 161         p = q
 162         if p is h:
 163             break
 164 # >>>
 165 def mp_solve_choices(p, q, n, delta_x, delta_y, delta, psi): # <<<
 166     """Implements mp_solve_choices form metapost (mp.c)"""
 167     uu = [None]*(len(delta)+1)
 168     vv = [None]*len(uu) # angles
 169     ww = [None]*len(uu)
 170     theta = [None]*len(uu)
 171     # k current knot number
 172     # r, s, t registers for list traversal
 173     k = 0
 174     s = p
 175     r = 0
 176     while True:
 177         t = s.next
 178         if k == 0: # <<<
 179             # 354:
 180             # Get the linear equations started
 181             # or return with the control points in place, if linear equations needn't be solved
 182
 183             if s.rtype == mp_given: # <<<
 184                 if t.ltype == mp_given:
 185                     # 372: Reduce to simple case of two givens and return
 186                     aa = mp_n_arg(delta_x[0], delta_y[0])
 187                     ct, st = mp_n_sin_cos(p.right_given() - aa)
 188                     cf, sf = mp_n_sin_cos(q.left_given() - aa)
 189                     mp_set_controls(p, q, delta_x[0], delta_y[0], st, ct, -sf, cf)
 190                     return
 191                 else:
 192                     # 362:
 193                     vv[0] = s.right_given() - mp_n_arg(delta_x[0], delta_y[0])
 194                     vv[0] = reduce_angle(vv[0])
 195                     uu[0] = 0
 196                     ww[0] = 0
 197             # >>>
 198             elif s.rtype == mp_curl: # <<<
 199                 if t.ltype == mp_curl:
 200                     # 373: (mp.pdf) Reduce to simple case of straight line and return
 201                     p.rtype = mp_explicit
 202                     q.ltype = mp_explicit
 203                     lt = abs(q.left_tension())
 204                     rt = abs(p.right_tension())
 205                     if rt == unity:
 206                         if delta_x[0] >= 0:
 207                             p.rx_pt = p.x_pt + delta_x[0]/3.0
 208                         else:
 209                             p.rx_pt = p.x_pt + delta_x[0]/3.0
 210                         if delta_y[0] >= 0:
 211                             p.ry_pt = p.y_pt + delta_y[0]/3.0
 212                         else:
 213                             p.ry_pt = p.y_pt + delta_y[0]/3.0
 214                     else:
 215                         ff = mp_make_fraction(unity, 3.0*rt)
 216                         p.rx_pt = p.x_pt + mp_take_fraction(delta_x[0], ff)
 217                         p.ry_pt = p.y_pt + mp_take_fraction(delta_y[0], ff)
 218
 219                     if lt == unity:
 220                         if delta_x[0] >= 0:
 221                             q.lx_pt = q.x_pt - delta_x[0]/3.0
 222                         else:
 223                             q.lx_pt = q.x_pt - delta_x[0]/3.0
 224                         if delta_y[0] >= 0:
 225                             q.ly_pt = q.y_pt - delta_y[0]/3.0
 226                         else:
 227                             q.ly_pt = q.y_pt - delta_y[0]/3.0
 228                     else:
 229                         ff = mp_make_fraction(unity, 3.0*lt)
 230                         q.lx_pt = q.x_pt - mp_take_fraction(delta_x[0], ff)
 231                         q.ly_pt = q.y_pt - mp_take_fraction(delta_y[0], ff)
 232                     return
 233
 234                 else: # t.ltype != mp_curl
 235                     # 363:
 236                     cc = s.right_curl()
 237                     lt = abs(t.left_tension())
 238                     rt = abs(s.right_tension())
 239                     if rt == unity and lt == unity: # XXX float comparison
 240                         uu[0] = mp_make_fraction(cc + cc + unity, cc + two)
 241                     else:
 242                         uu[0] = mp_curl_ratio(cc, rt, lt)
 243                     vv[0] = -mp_take_fraction(psi[1], uu[0])
 244                     ww[0] = 0
 245             # >>>
 246             elif s.rtype == mp_open: # <<<
 247                 uu[0] = 0
 248                 vv[0] = 0
 249                 ww[0] = fraction_one
 250             # >>>
 251         # end of 354 >>>
 252         else: # k > 0 <<<
 253
 254             if s.ltype == mp_end_cycle or s.ltype == mp_open: # <<<
 255                 # 356: Set up equation to match mock curvatures at z_k;
 256                 #      then finish loop with theta_n adjusted to equal theta_0, if a
 257                 #      cycle has ended
 258
 259                 # 357: Calculate the values
 260                 #      aa = Ak/Bk, bb = Dk/Ck, dd = (3-alpha_{k-1})d(k,k+1),
 261                 #      ee = (3-beta_{k+1})d(k-1,k), cc=(Bk-uk-Ak)/Bk
 262                 if abs(r.right_tension()) == unity: # XXX float comparison
 263                     aa = fraction_half
 264                     dd = 2*delta[k]
 265                 else:
 266                     aa = mp_make_fraction(unity, 3*abs(r.right_tension()) - unity)
 267                     dd = mp_take_fraction(delta[k],
 268                                           fraction_three - mp_make_fraction(unity, abs(r.right_tension())))
 269                 if abs(t.left_tension()) == unity: # XXX float comparison
 270                     bb = fraction_half
 271                     ee = 2*delta[k-1]
 272                 else:
 273                     bb = mp_make_fraction(unity, 3*abs(t.left_tension()) - unity)
 274                     ee = mp_take_fraction(delta[k-1],
 275                                           fraction_three - mp_make_fraction(unity, abs(t.left_tension())))
 276                 cc = fraction_one - mp_take_fraction(uu[k-1], aa)
 277
 278                 # 358: Calculate the ratio ff = Ck/(Ck + Bk - uk-1Ak)
 279                 dd = mp_take_fraction(dd, cc)
 280                 lt = abs(s.left_tension())
 281                 rt = abs(s.right_tension())
 282                 if lt != rt:
 283                     if lt < rt:
 284                         ff = mp_make_fraction(lt, rt)
 285                         ff = mp_take_fraction(ff, ff)
 286                         dd = mp_take_fraction(dd, ff)
 287                     else:
 288                         ff = mp_make_fraction(rt, lt)
 289                         ff = mp_take_fraction(ff, ff)
 290                         ee = mp_take_fraction(ee, ff)
 291                 ff = mp_make_fraction(ee, ee + dd)
 292
 293                 uu[k] = mp_take_fraction(ff, bb)
 294
 295                 # 359: Calculate the values of vk and wk
 296                 acc = -mp_take_fraction(psi[k+1], uu[k])
 297                 if r.rtype == mp_curl:
 298                     ww[k] = 0
 299                     vv[k] = acc - mp_take_fraction(psi[1], fraction_one - ff)
 300                 else:
 301                     ff = mp_make_fraction(fraction_one - ff, cc)
 302                     acc = acc - mp_take_fraction(psi[k], ff)
 303                     ff = mp_take_fraction(ff, aa)
 304                     vv[k] = acc - mp_take_fraction(vv[k-1], ff)
 305                     if ww[k-1] == 0:
 306                         ww[k] = 0
 307                     else:
 308                         ww[k] = -mp_take_fraction(ww[k-1], ff)
 309
 310                 if s.ltype == mp_end_cycle:
 311                     # 360: Adjust theta_n to equal theta_0 and finish loop
 312
 313                     aa = 0
 314                     bb = fraction_one
 315                     while True:
 316                         k -= 1
 317                         if k == 0:
 318                             k = n
 319                         aa = vv[k] - mp_take_fraction(aa, uu[k])
 320                         bb = ww[k] - mp_take_fraction(bb, uu[k])
 321                         if k == n:
 322                             break
 323                     aa = mp_make_fraction(aa, fraction_one - bb)
 324                     theta[n] = aa
 325                     vv[0] = aa
 326                     for k in range(1, n):
 327                         vv[k] = vv[k] + mp_take_fraction(aa, ww[k])
 328                     break
 329             # >>>
 330             elif s.ltype == mp_curl: # <<<
 331                 # 364:
 332                 cc = s.left_curl()
 333                 lt = abs(s.left_tension())
 334                 rt = abs(r.right_tension())
 335                 if rt == unity and lt == unity:
 336                     ff = mp_make_fraction(cc + cc + unity, cc + two)
 337                 else:
 338                     ff = mp_curl_ratio(cc, lt, rt)
 339                 theta[n] = -mp_make_fraction(mp_take_fraction(vv[n-1], ff),
 340                                              fraction_one - mp_take_fraction(ff, uu[n-1]))
 341                 break
 342             # >>>
 343             elif s.ltype == mp_given: # <<<
 344                 # 361:
 345                 theta[n] = s.left_given() - mp_n_arg(delta_x[n-1], delta_y[n-1])
 346                 theta[n] = reduce_angle(theta[n])
 347                 break
 348             # >>>
 349
 350         # end of k == 0, k != 0 >>>
 351
 352         r = s
 353         s = t
 354         k += 1
 355
 356     # 367:
 357     # Finish choosing angles and assigning control points
 358     for k in range(n-1, -1, -1):
 359         theta[k] = vv[k] - mp_take_fraction(theta[k+1], uu[k])
 360     s = p
 361     k = 0
 362     while True:
 363         t = s.next
 364         ct, st = mp_n_sin_cos(theta[k])
 365         cf, sf = mp_n_sin_cos(-psi[k+1]-theta[k+1])
 366         mp_set_controls(s, t, delta_x[k], delta_y[k], st, ct, sf, cf)
 367         k += 1
 368         s = t
 369         if k == n:
 370             break
 371 # >>>
 372 def mp_n_arg(x, y): # <<<
 373     return atan2(y, x)
 374 # >>>
 375 def mp_n_sin_cos(z): # <<<
 376     """Given an integer z that is 2**20 times an angle theta in degrees, the
 377     purpose of n sin cos(z) is to set x = r cos theta and y = r sin theta
 378     (approximately), for some rather large number r. The maximum of x and y
 379     will be between 2**28 and 2**30, so that there will be hardly any loss of
 380     accuracy. Then x and y are divided by r."""
 381     # 67: mpmath.pdf
 382     return cos(z), sin(z)
 383 # >>>
 384 def mp_set_controls(p, q, delta_x, delta_y, st, ct, sf, cf): # <<<
 385     """The set controls routine actually puts the control points into a pair of
 386     consecutive nodes p and q. Global variables are used to record the values
 387     of sin(theta), cos(theta), sin(phi), and cos(phi) needed in this
 388     calculation.
 389
 390     See mp.pdf, item 370"""
 391     lt = abs(q.left_tension())
 392     rt = abs(p.right_tension())
 393     rr = mp_velocity(st, ct, sf, cf, rt)
 394     ss = mp_velocity(sf, cf, st, ct, lt)
 395     if p.right_tension() < 0 or q.left_tension() < 0:
 396         # 371: Decrease the velocities, if necessary, to stay inside the bounding triangle
 397         # this is the only place where the sign of the tension counts
 398         if (st >= 0 and sf >= 0) or (st <= 0 and sf <= 0):
 399             sine = mp_take_fraction(abs(st), cf) + mp_take_fraction(abs(sf), ct) # sin(theta+phi)
 400             if sine > 0:
 401                 #sine = mp_take_fraction(sine, fraction_one + unity) # safety factor
 402                 sine *= 1.00024414062 # safety factor
 403                 if p.right_tension() < 0:
 404                     if mp_ab_vs_cd(abs(sf), fraction_one, rr, sine) < 0:
 405                         rr = mp_make_fraction(abs(sf), sine)
 406                 if q.left_tension() < 0:
 407                     if mp_ab_vs_cd(abs(st), fraction_one, ss, sine) < 0:
 408                         ss = mp_make_fraction(abs(st), sine)
 409
 410     p.rx_pt = p.x_pt + mp_take_fraction(mp_take_fraction(delta_x, ct) - mp_take_fraction(delta_y, st), rr)
 411     p.ry_pt = p.y_pt + mp_take_fraction(mp_take_fraction(delta_y, ct) + mp_take_fraction(delta_x, st), rr)
 412     q.lx_pt = q.x_pt - mp_take_fraction(mp_take_fraction(delta_x, cf) + mp_take_fraction(delta_y, sf), ss)
 413     q.ly_pt = q.y_pt - mp_take_fraction(mp_take_fraction(delta_y, cf) - mp_take_fraction(delta_x, sf), ss)
 414     p.rtype = mp_explicit
 415     q.ltype = mp_explicit
 416 # >>>
 417 def mp_make_fraction(p, q): # <<<
 418     # 17: mpmath.pdf
 419     """The make fraction routine produces the fraction equivalent of p/q, given
 420     integers p and q; it computes the integer f = floor(2**28 p/q + 1/2), when
 421     p and q are positive.
 422
 423     In machine language this would simply be (2**28*p)div q"""
 424     return p/q
 425 # >>>
 426 def mp_take_fraction(q, f): # <<<
 427     # 20: mpmath.pdf
 428     """The dual of make fraction is take fraction, which multiplies a given
 429     integer q by a fraction f. When the operands are positive, it computes
 430     p = floor(q*f/2**28 + 1/2), a symmetric function of q and f."""
 431     return q*f
 432 # >>>
 433 def mp_pyth_add(a, b): # <<<
 434     # 44: mpmath.pdf
 435     """Pythagorean addition sqrt(a**2 + b**2) is implemented by an elegant
 436     iterative scheme due to Cleve Moler and Donald Morrison [IBM Journal of
 437     Research and Development 27 (1983), 577-581]. It modifies a and b in such a
 438     way that their Pythagorean sum remains invariant, while the smaller
 439     argument decreases."""
 440     return sqrt(a*a + b*b)
 441 # >>>
 442 def mp_curl_ratio(gamma, a_tension, b_tension): # <<<
 443     """The curl ratio subroutine has three arguments, which our previous
 444     notation encourages us to call gamma, 1/alpha, and 1/beta. It is a somewhat
 445     tedious program to calculate
 446       [(3-alpha)alpha^2 gamma + beta^3] / [alpha^3 gamma + (3-beta)beta^2],
 447     with the result reduced to 4 if it exceeds 4. (This reduction of curl is
 448     necessary only if the curl and tension are both large.) The values of alpha
 449     and beta will be at most 4/3.
 450
 451     See mp.pdf (items 365, 366)."""
 452     alpha = 1.0/a_tension
 453     beta = 1.0/b_tension
 454     return min(4.0, ((3.0-alpha)*alpha**2*gamma + beta**3) /
 455                     (alpha**3*gamma + (3.0-beta)*beta**2))
 456 # >>>
 457 def mp_ab_vs_cd(a, b, c, d): # <<<
 458     """Tests rigorously if ab is greater than, equal to, or less than cd, given
 459     integers (a, b, c, d). In most cases a quick decision is reached. The
 460     result is +1, 0, or -1 in the three respective cases.
 461     See mpmath.pdf (item 33)."""
 462     if a*b == c*d:
 463         return 0
 464     if a*b > c*d:
 465         return 1
 466     return -1
 467 # >>>
 468 def mp_velocity(st, ct, sf, cf, t): # <<<
 469     """Metapost's standard velocity subroutine for cubic Bezier curves.
 470     See mpmath.pdf (item 30) and mp.pdf (item 339)."""
 471     return min(4.0, (2.0 + sqrt(2)*(st-sf/16.0)*(sf-st/16.0)*(ct-cf)) /
 472                     (1.5*t*(2+(sqrt(5)-1)*ct + (3-sqrt(5))*cf)))
 473 # >>>
 474 def reduce_angle(A): # <<<
 475     """A macro in mp.c"""
 476     if abs(A) > one_eighty_deg:
 477         if A > 0:
 478             A -= three_sixty_deg
 479         else:
 480             A += three_sixty_deg
 481     return A
 482 # >>>
 483
 484 # vim:foldmethod=marker:foldmarker=<<<,>>>