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