pyx/deformer.py

   1 #!/usr/bin/env python
   2 # -*- coding: ISO-8859-1 -*-
   3 #
   4 #
   5 # Copyright (C) 2003-2005 Michael Schindler <m-schindler@users.sourceforge.net>
   6 # Copyright (C) 2003-2004 André Wobst <wobsta@users.sourceforge.net>
   7 #
   8 # This file is part of PyX (http://pyx.sourceforge.net/).
   9 #
  10 # PyX is free software; you can redistribute it and/or modify
  11 # it under the terms of the GNU General Public License as published by
  12 # the Free Software Foundation; either version 2 of the License, or
  13 # (at your option) any later version.
  14 #
  15 # PyX is distributed in the hope that it will be useful,
  16 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18 # GNU General Public License for more details.
  19 #
  20 # You should have received a copy of the GNU General Public License
  21 # along with PyX; if not, write to the Free Software
  22 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  23
  24
  25 import math, warnings
  26 import attr, color, helper, path, style, trafo, unit
  27 try:
  28   import Numeric, LinearAlgebra
  29   have_la_packages = 1
  30
  31   def realpolyroots(coeffs, epsilon=1e-5): # <<<
  32
  33       """returns the roots of a polynom with given coefficients
  34
  35       This helper routine uses the package Numeric to find the roots
  36       of the polynomial with coefficients given in coeffs:
  37         0 = \sum_{i=0}^N x^{N-i} coeffs[i]
  38       The solution is found via an equivalent eigenvalue problem
  39       """
  40
  41       try:
  42           1.0 / coeffs[0]
  43       except:
  44           return realpolyroots(coeffs[1:], epsilon=epsilon)
  45       else:
  46
  47           N = len(coeffs)
  48           # build the Matrix of the polynomial problem
  49           mat = Numeric.zeros((N, N), Numeric.Float)
  50           for i in range(N-1):
  51               mat[i+1][i] = 1
  52           for i in range(N-1):
  53               mat[0][i] = -coeffs[i+1]/coeffs[0]
  54           # find the eigenvalues of the matrix (== the zeros of the polynomial)
  55           zeros = [complex(zero) for zero in LinearAlgebra.eigenvalues(mat)]
  56           # take only the real zeros
  57           zeros = [zero.real for zero in zeros if -epsilon < zero.imag < epsilon]
  58
  59           ## check if the zeros are really zeros!
  60           #for zero in zeros:
  61           #    p = 0
  62           #    for i in range(N):
  63           #        p += coeffs[i] * zero**(N-i)
  64           #    if abs(p) > epsilon:
  65           #        raise Exception("value %f instead of 0" % p)
  66
  67       return zeros
  68   # >>>
  69
  70 except ImportError:
  71   have_la_packages = 0
  72
  73
  74 def sign1(x):
  75     return (x >= 0) and 1 or -1
  76
  77 def curvescontrols_from_endlines_pt (B, tangent1, tangent2, r1, r2, softness): # <<<
  78     # calculates the parameters for two bezier curves connecting two lines (curvature=0)
  79     # starting at B - r1*tangent1
  80     # ending at   B + r2*tangent2
  81     #
  82     # Takes the corner B
  83     # and two tangent vectors heading to and from B
  84     # and two radii r1 and r2:
  85     # All arguments must be in Points
  86     # Returns the seven control points of the two bezier curves:
  87     #  - start d1
  88     #  - control points g1 and f1
  89     #  - midpoint e
  90     #  - control points f2 and g2
  91     #  - endpoint d2
  92
  93     # make direction vectors d1: from B to A
  94     #                        d2: from B to C
  95     d1 = -tangent1[0] / math.hypot(*tangent1), -tangent1[1] / math.hypot(*tangent1)
  96     d2 =  tangent2[0] / math.hypot(*tangent2),  tangent2[1] / math.hypot(*tangent2)
  97
  98     # 0.3192 has turned out to be the maximum softness available
  99     # for straight lines ;-)
 100     f = 0.3192 * softness
 101     g = (15.0 * f + math.sqrt(-15.0*f*f + 24.0*f))/12.0
 102
 103     # make the control points of the two bezier curves
 104     f1 = B[0] + f * r1 * d1[0], B[1] + f * r1 * d1[1]
 105     f2 = B[0] + f * r2 * d2[0], B[1] + f * r2 * d2[1]
 106     g1 = B[0] + g * r1 * d1[0], B[1] + g * r1 * d1[1]
 107     g2 = B[0] + g * r2 * d2[0], B[1] + g * r2 * d2[1]
 108     d1 = B[0] +     r1 * d1[0], B[1] +     r1 * d1[1]
 109     d2 = B[0] +     r2 * d2[0], B[1] +     r2 * d2[1]
 110     e  = 0.5 * (f1[0] + f2[0]), 0.5 * (f1[1] + f2[1])
 111
 112     return (d1, g1, f1, e, f2, g2, d2)
 113 # >>>
 114
 115 def controldists_from_endpoints_pt (A, B, tangentA, tangentB, curvA, curvB, epsilon=1e-5): # <<<
 116
 117     """distances for a curve given by tangents and curvatures at the endpoints
 118
 119     This helper routine returns the two distances between the endpoints and the
 120     corresponding control points of a (cubic) bezier curve that has
 121     prescribed tangents tangentA, tangentB and curvatures curvA, curvB at the
 122     end points.
 123     """
 124
 125     # normalise the tangent vectors
 126     normA = math.hypot(*tangentA)
 127     tangA = (tangentA[0] / normA, tangentA[1] / normA)
 128     normB = math.hypot(*tangentB)
 129     tangB = (tangentB[0] / normB, tangentB[1] / normB)
 130     # some shortcuts
 131     T = tangB[0] * tangA[1] - tangB[1] * tangA[0]
 132     D =   tangA[0] * (B[1]-A[1]) - tangA[1] * (B[0]-A[0])
 133     E = - tangB[0] * (B[1]-A[1]) + tangB[1] * (B[0]-A[0])
 134     # the variables: \dot X(0) = 3 * a * tangA
 135     #                \dot X(1) = 3 * b * tangB
 136     a, b = None, None
 137
 138
 139     # try some special cases where the equations decouple
 140     try:
 141         1.0 / T
 142     except ZeroDivisionError:
 143         try:
 144             a = math.sqrt(2.0 * D / (3.0 * curvA))
 145             b = math.sqrt(2.0 * E / (3.0 * curvB))
 146         except ZeroDivisionError:
 147             a = b = None
 148         except ValueError:
 149             raise # ???
 150     else:
 151         try:
 152             1.0 / curvA
 153         except ZeroDivisionError:
 154             b = -D / T
 155             a = (1.5*curvB*b*b - E) / T
 156         else:
 157             try:
 158                 1.0 / curvB
 159             except ZeroDivisionError:
 160                 a = -E / T
 161                 b = (1.5*curvA*a*a - D) / T
 162             else:
 163                 a, b = None, None
 164
 165
 166     # else find a solution for the full problem
 167     if a is None:
 168         if have_la_packages:
 169             # we first try to find all the zeros of the polynomials for a or b (4th order)
 170             # this needs Numeric and LinearAlgebra
 171
 172             coeffs_a = (3.375*curvA*curvA*curvB, 0, -4.5*curvA*curvB*D, -T**3,  1.5*curvB*D*D - T*T*E)
 173             coeffs_b = (3.375*curvA*curvB*curvB, 0, -4.5*curvA*curvB*E, -T**3,  1.5*curvA*E*E - T*T*D)
 174
 175             # First try the equation for a
 176             cands_a = [cand for cand in realpolyroots(coeffs_a) if cand > 0]
 177
 178             if cands_a:
 179                 a = min(cands_a)
 180                 b = (1.5*curvA*a*a - D) / T
 181             else:
 182                 # then, try the equation for b
 183                 cands_b = [cand for cand in realpolyroots(coeffs_b) if cand > 0]
 184                 if cands_b:
 185                     b = min(cands_b)
 186                     a = (1.5*curvB*b*b - E) / T
 187                 else:
 188                     a = b = None
 189
 190         else:
 191             # if the Numeric modules are not available:
 192             # solve the coupled system by Newton iteration
 193             #     0 = Ga(a,b) = 0.5 a |a| curvA + b * T - D
 194             #     0 = Gb(a,b) = 0.5 b |b| curvB + a * T + E
 195             # this system is equivalent to the geometric contraints:
 196             #     the curvature and the normal tangent vectors
 197             #     at parameters 0 and 1 are to be continuous
 198             # the system is solved by 2-dim Newton-Iteration
 199             # (a,b)^{i+1} = (a,b)^i - (DG)^{-1} (Ga(a^i,b^i), Gb(a^i,b^i))
 200             a = 1.0 / abs(curvA)
 201             b = 1.0 / abs(curvB)
 202             eps = 1.0 # stepwith for the Newton iteration
 203             da = db = 2*epsilon
 204             counter = 0
 205             while max(abs(da),abs(db)) > epsilon and counter < 1000:
 206
 207                 Ga = eps * (1.5*curvA*a*a - T*b - D)
 208                 Gb = eps * (1.5*curvB*b*b - T*a - E)
 209
 210                 detDG = 9.0*a*b*curvA*curvB - T*T
 211                 invDG = ((3.0*curvB*b/detDG, T/detDG), (T/detDG, 3.0*curvA*a/detDG))
 212
 213                 da = invDG[0][0] * Ga + invDG[0][1] * Gb
 214                 db = invDG[1][0] * Ga + invDG[1][1] * Gb
 215
 216                 a -= da
 217                 b -= db
 218
 219                 counter += 1
 220
 221     if a < 0 or b < 0:
 222         a = b = None
 223
 224     if a is not None: a /= normA
 225     if b is not None: b /= normB
 226
 227     return a, b
 228 # >>>
 229
 230 def intersection (A, D, tangA, tangD): # <<<
 231
 232     """returns the intersection parameters of two evens
 233
 234     they are defined by:
 235       x(t) = A + t * tangA
 236       x(s) = D + s * tangD
 237     """
 238     det = -tangA[0] * tangD[1] + tangA[1] * tangD[0]
 239     try:
 240         1.0 / det
 241     except ArithmeticError:
 242         return None, None
 243
 244     DA = D[0] - A[0], D[1] - A[1]
 245
 246     t = (-tangD[1]*DA[0] + tangD[0]*DA[1]) / det
 247     s = (-tangA[1]*DA[0] + tangA[0]*DA[1]) / det
 248
 249     return t, s
 250 # >>>
 251
 252 def parallel_curvespoints_pt (orig_ncurve, shift, expensive=0, relerr=0.05, epsilon=1e-5, counter=1): # <<<
 253
 254     A = orig_ncurve.x0_pt, orig_ncurve.y0_pt
 255     B = orig_ncurve.x1_pt, orig_ncurve.y1_pt
 256     C = orig_ncurve.x2_pt, orig_ncurve.y2_pt
 257     D = orig_ncurve.x3_pt, orig_ncurve.y3_pt
 258
 259     # non-normalized tangential vector
 260     # from begin/end point to the corresponding controlpoint
 261     tangA = (B[0] - A[0], B[1] - A[1])
 262     tangD = (D[0] - C[0], D[1] - C[1])
 263
 264     # normalized normal vectors
 265     # turned to the left (+90 degrees) from the tangents
 266     NormA = (-tangA[1] / math.hypot(*tangA), tangA[0] / math.hypot(*tangA))
 267     NormD = (-tangD[1] / math.hypot(*tangD), tangD[0] / math.hypot(*tangD))
 268
 269     # radii of curvature
 270     radiusA, radiusD = orig_ncurve.curveradius_pt([0,1])
 271
 272     # get the new begin/end points
 273     A = A[0] + shift * NormA[0], A[1] + shift * NormA[1]
 274     D = D[0] + shift * NormD[0], D[1] + shift * NormD[1]
 275
 276     try:
 277         if radiusA is None:
 278             curvA = 0
 279         else:
 280             curvA = 1.0 / (radiusA - shift)
 281         if radiusD is None:
 282             curvD = 0
 283         else:
 284             curvD = 1.0 / (radiusD - shift)
 285     except ZeroDivisionError:
 286         raise
 287     else:
 288         a, d = controldists_from_endpoints_pt (A, D, tangA, tangD, curvA, curvD, epsilon=epsilon)
 289
 290         if a is None or d is None:
 291             # fallback heuristic
 292             a = (radiusA - shift) / radiusA
 293             d = (radiusD - shift) / radiusD
 294
 295         B = A[0] + a * tangA[0], A[1] + a * tangA[1]
 296         C = D[0] - d * tangD[0], D[1] - d * tangD[1]
 297
 298     controlpoints = [(A,B,C,D)]
 299
 300     # check if the distance is really the wanted distance
 301     if expensive and counter < 10:
 302         # measure the distance in the "middle" of the original curve
 303         trafo = orig_ncurve.trafo([0.5])[0]
 304         M = trafo._apply(0,0)
 305         NormM = trafo._apply(0,1)
 306         NormM = NormM[0] - M[0], NormM[1] - M[1]
 307
 308         nline = path.normline_pt (
 309           M[0] + (1.0 - 2*relerr) * shift * NormM[0],
 310           M[1] + (1.0 - 2*relerr) * shift * NormM[1],
 311           M[0] + (1.0 + 2*relerr) * shift * NormM[0],
 312           M[1] + (1.0 + 2*relerr) * shift * NormM[1])
 313
 314         new_ncurve = path.normcurve_pt(A[0],A[1], B[0],B[1], C[0],C[1], D[0],D[1])
 315
 316         #cutparams = nline.intersect(orig_ncurve, epsilon)
 317         cutparams = new_ncurve.intersect(nline, epsilon)
 318         if cutparams:
 319             cutpoints = nline.at_pt(cutparams[0])
 320         else:
 321             cutpoints = []
 322         good = 0
 323         for cutpoint in cutpoints:
 324             if cutpoint is not None:
 325                 dist = math.hypot(M[0] - cutpoint[0], M[1] - cutpoint[1])
 326                 if abs(dist - abs(shift)) < relerr * abs(shift):
 327                     good = 1
 328
 329         if not good:
 330             first, second = orig_ncurve.segments([0,0.5,1])
 331             controlpoints = \
 332               parallel_curvespoints_pt (first,  shift, expensive, relerr, epsilon, counter+1) + \
 333               parallel_curvespoints_pt (second, shift, expensive, relerr, epsilon, counter+1)
 334
 335
 336
 337     # TODO:
 338     # too big curvatures: intersect curves
 339     # there is something wrong with the recursion
 340     return controlpoints
 341 # >>>
 342
 343
 344 class deformer(attr.attr):
 345
 346     def deform (self, basepath):
 347         return basepath
 348
 349 class cycloid(deformer): # <<<
 350     """Wraps a cycloid around a path.
 351
 352     The outcome looks like a metal spring with the originalpath as the axis.
 353     radius: radius of the cycloid
 354     loops:  number of loops from beginning to end of the original path
 355     skipfirst/skiplast: undeformed end lines of the original path
 356
 357     """
 358
 359     def __init__(self, radius=0.5*unit.t_cm, halfloops=10,
 360     skipfirst=1*unit.t_cm, skiplast=1*unit.t_cm, curvesperhloop=3, sign=1, turnangle=45):
 361         self.skipfirst = skipfirst
 362         self.skiplast = skiplast
 363         self.radius = radius
 364         self.halfloops = halfloops
 365         self.curvesperhloop = curvesperhloop
 366         self.sign = sign
 367         self.turnangle = turnangle
 368
 369     def __call__(self, radius=None, halfloops=None,
 370     skipfirst=None, skiplast=None, curvesperhloop=None, sign=None, turnangle=None):
 371         if radius is None:
 372             radius = self.radius
 373         if halfloops is None:
 374             halfloops = self.halfloops
 375         if skipfirst is None:
 376             skipfirst = self.skipfirst
 377         if skiplast is None:
 378             skiplast = self.skiplast
 379         if curvesperhloop is None:
 380             curvesperhloop = self.curvesperhloop
 381         if sign is None:
 382             sign = self.sign
 383         if turnangle is None:
 384             turnangle = self.turnangle
 385
 386         return cycloid(radius=radius, halfloops=halfloops, skipfirst=skipfirst, skiplast=skiplast,
 387                        curvesperhloop=curvesperhloop, sign=sign, turnangle=turnangle)
 388
 389     def deform(self, basepath):
 390         resultnormsubpaths = [self.deformsubpath(nsp) for nsp in basepath.normpath().normsubpaths]
 391         return path.normpath(resultnormsubpaths)
 392
 393     def deformsubpath(self, normsubpath):
 394
 395         skipfirst = abs(unit.topt(self.skipfirst))
 396         skiplast = abs(unit.topt(self.skiplast))
 397         radius = abs(unit.topt(self.radius))
 398         turnangle = self.turnangle * math.pi / 180.0
 399
 400         cosTurn = math.cos(turnangle)
 401         sinTurn = math.sin(turnangle)
 402
 403         # make list of the lengths and parameters at points on normsubpath where we will add cycloid-points
 404         totlength = normsubpath.arclen_pt()
 405         if totlength <= skipfirst + skiplast + 2*radius*sinTurn:
 406             warnings.warn("normsubpath is too short for deformation with cycloid -- skipping...")
 407             return normsubpath
 408
 409         # parametrisation is in rotation-angle around the basepath
 410         # differences in length, angle ... between two basepoints
 411         # and between basepoints and controlpoints
 412         Dphi = math.pi / self.curvesperhloop
 413         phis = [i * Dphi for i in range(self.halfloops * self.curvesperhloop + 1)]
 414         DzDphi = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * math.pi * cosTurn)
 415         # Dz = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * self.curvesperhloop * cosTurn)
 416         # zs = [i * Dz for i in range(self.halfloops * self.curvesperhloop + 1)]
 417         # from path._arctobcurve:
 418         # optimal relative distance along tangent for second and third control point
 419         L = 4 * radius * (1 - math.cos(Dphi/2)) / (3 * math.sin(Dphi/2))
 420
 421         # Now the transformation of z into the turned coordinate system
 422         Zs = [ skipfirst + radius*sinTurn # here the coordinate z starts
 423              - sinTurn*radius*math.cos(phi) + cosTurn*DzDphi*phi # the transformed z-coordinate
 424              for phi in phis]
 425         params = normsubpath._arclentoparam_pt(Zs)[0]
 426
 427         # get the positions of the splitpoints in the cycloid
 428         points = []
 429         for phi, param in zip(phis, params):
 430             # the cycloid is a circle that is stretched along the normsubpath
 431             # here are the points of that circle
 432             basetrafo = normsubpath.trafo([param])[0]
 433
 434             # The point on the cycloid, in the basepath's local coordinate system
 435             baseZ, baseY = 0, radius*math.sin(phi)
 436
 437             # The tangent there, also in local coords
 438             tangentX = -cosTurn*radius*math.sin(phi) + sinTurn*DzDphi
 439             tangentY = radius*math.cos(phi)
 440             tangentZ = sinTurn*radius*math.sin(phi) + DzDphi*cosTurn
 441             norm = math.sqrt(tangentX*tangentX + tangentY*tangentY + tangentZ*tangentZ)
 442             tangentY, tangentZ = tangentY/norm, tangentZ/norm
 443
 444             # Respect the curvature of the basepath for the cycloid's curvature
 445             # XXX this is only a heuristic, not a "true" expression for
 446             #     the curvature in curved coordinate systems
 447             pathradius = normsubpath.curveradius_pt([param])[0]
 448             if pathradius is not None:
 449                 factor = (pathradius - baseY) / pathradius
 450                 factor = abs(factor)
 451             else:
 452                 factor = 1
 453             l = L * factor
 454
 455             # The control points prior and after the point on the cycloid
 456             preeZ, preeY = baseZ - l * tangentZ, baseY - l * tangentY
 457             postZ, postY = baseZ + l * tangentZ, baseY + l * tangentY
 458
 459             # Now put everything at the proper place
 460             points.append(basetrafo._apply(preeZ, self.sign * preeY) +
 461                           basetrafo._apply(baseZ, self.sign * baseY) +
 462                           basetrafo._apply(postZ, self.sign * postY))
 463
 464         if len(points) <= 1:
 465             warnings.warn("normsubpath is too short for deformation with cycloid -- skipping...")
 466             return normsubpath
 467
 468         # Build the path from the pointlist
 469         # containing (control x 2,  base x 2, control x 2)
 470         if skipfirst > normsubpath.epsilon:
 471             normsubpathitems = normsubpath.segments([0, params[0]])[0]
 472             normsubpathitems.append(path.normcurve_pt(*(points[0][2:6] + points[1][0:4])))
 473         else:
 474             normsubpathitems = [path.normcurve_pt(*(points[0][2:6] + points[1][0:4]))]
 475         for i in range(1, len(points)-1):
 476             normsubpathitems.append(path.normcurve_pt(*(points[i][2:6] + points[i+1][0:4])))
 477         if skiplast > normsubpath.epsilon:
 478             for nsp in normsubpath.segments([params[-1], len(normsubpath)]):
 479                 normsubpathitems.extend(nsp.normsubpathitems)
 480
 481         # That's it
 482         return path.normsubpath(normsubpathitems, epsilon=normsubpath.epsilon)
 483 # >>>
 484
 485 cycloid.clear = attr.clearclass(cycloid)
 486
 487 class smoothed(deformer): # <<<
 488
 489     """Bends corners in a path.
 490
 491     This decorator replaces corners in a path with bezier curves. There are two cases:
 492     - If the corner lies between two lines, _two_ bezier curves will be used
 493       that are highly optimized to look good (their curvature is to be zero at the ends
 494       and has to have zero derivative in the middle).
 495       Additionally, it can controlled by the softness-parameter.
 496     - If the corner lies between curves then _one_ bezier is used that is (except in some
 497       special cases) uniquely determined by the tangents and curvatures at its end-points.
 498       In some cases it is necessary to use only the absolute value of the curvature to avoid a
 499       cusp-shaped connection of the new bezier to the old path. In this case the use of
 500       "obeycurv=0" allows the sign of the curvature to switch.
 501     - The radius argument gives the arclength-distance of the corner to the points where the
 502       old path is cut and the beziers are inserted.
 503     - Path elements that are too short (shorter than the radius) are skipped
 504     """
 505
 506     def __init__(self, radius, softness=1, obeycurv=0):
 507         self.radius = radius
 508         self.softness = softness
 509         self.obeycurv = obeycurv
 510
 511     def __call__(self, radius=None, softness=None, obeycurv=None):
 512         if radius is None:
 513             radius = self.radius
 514         if softness is None:
 515             softness = self.softness
 516         if obeycurv is None:
 517             obeycurv = self.obeycurv
 518         return smoothed(radius=radius, softness=softness, obeycurv=obeycurv)
 519
 520     def deform(self, basepath):
 521         basepath = basepath.normpath()
 522         smoothpath = path.path()
 523
 524         for sp in basepath.normsubpaths:
 525             smoothpath += self.deformsubpath(sp)
 526
 527         return smoothpath
 528
 529     def deformsubpath(self, normsubpath):
 530
 531         radius = unit.topt(self.radius)
 532         epsilon = normsubpath.epsilon
 533
 534         # 1. Build up a list of all relevant normsubpathitems
 535         #    and the lengths where they will be cut (length with respect to the normsubpath)
 536         rel_npitems, arclengths = [], []
 537         for npitem in normsubpath.normsubpathitems:
 538
 539             arclength = npitem.arclen_pt(epsilon)
 540
 541             # items that should be totally skipped:
 542             # (we throw away the possible ending "closepath" piece of zero length)
 543             if (arclength > radius):
 544                 rel_npitems.append(npitem)
 545                 arclengths.append(arclength)
 546             else:
 547                 warnings.warn("smoothed is skipping a too short normsubpathitem")
 548
 549         # 2. Find the parameters, points,
 550         #    and calculate tangents and curvatures
 551         params, points, tangents, curvatures = [], [], [], []
 552         for npitem, arclength in zip(rel_npitems, arclengths):
 553
 554             # find the parameter(s): either one or two
 555             # for short items we squeeze the radius
 556             if arclength > 2 * radius:
 557                 cut_alengths = [radius, arclength - radius]
 558             else:
 559                 cut_alengths = [0.5 * radius]
 560
 561             # get the parameters
 562             pars = npitem._arclentoparam_pt(cut_alengths, epsilon)[0]
 563
 564             # the endpoints of an open subpath must be handled specially
 565             if not normsubpath.closed:
 566                 if npitem is rel_npitems[0]:
 567                     pars[0] = 0
 568                 if npitem is rel_npitems[-1]:
 569                     pars[-1] = 1
 570
 571             # find points, tangents and curvatures
 572             ts,cs,ps = [],[],[]
 573             for par in pars:
 574                 thetrafo = npitem.trafo([par])[0]
 575                 p = thetrafo._apply(0,0)
 576                 t = thetrafo._apply(1,0)
 577                 ps.append(p)
 578                 ts.append((t[0]-p[0], t[1]-p[1]))
 579                 r = npitem.curveradius_pt([par])[0]
 580                 if r is None:
 581                     cs.append(0)
 582                 else:
 583                     cs.append(1.0 / r)
 584
 585             params.append(pars)
 586             points.append(ps)
 587             tangents.append(ts)
 588             curvatures.append(cs)
 589
 590
 591         # create an empty path to collect pathitems
 592         # this will be returned as normpath, later
 593         smoothpath = path.path()
 594         do_moveto = 1 # we do not know yet where to moveto
 595
 596         # 3. Do the splitting for the first to the last element,
 597         #    a closed path must be closed later
 598         #
 599         for i in range(len(rel_npitems)):
 600
 601             this = i
 602             next = (i+1) % len(rel_npitems)
 603             thisnpitem = rel_npitems[this]
 604             nextnpitem = rel_npitems[next]
 605
 606             # split thisnpitem apart and take the middle piece
 607             # We start the new path with the middle piece of the first path-elem
 608             if len(points[this]) == 2:
 609
 610                 middlepiece = thisnpitem.segments(params[this])[0]
 611
 612                 if do_moveto:
 613                     smoothpath.append(path.moveto_pt(*middlepiece.atbegin_pt()))
 614                     do_moveto = 0
 615
 616                 if isinstance(middlepiece, path.normline_pt):
 617                     smoothpath.append(path.lineto_pt(*middlepiece.atend_pt()))
 618                 elif isinstance(middlepiece, path.normcurve_pt):
 619                     smoothpath.append(path.curveto_pt(
 620                       middlepiece.x1_pt, middlepiece.y1_pt,
 621                       middlepiece.x2_pt, middlepiece.y2_pt,
 622                       middlepiece.x3_pt, middlepiece.y3_pt))
 623
 624             if (not normsubpath.closed) and (thisnpitem is rel_npitems[-1]):
 625                 continue
 626
 627             # add the curve(s) replacing the corner
 628             if isinstance(thisnpitem, path.normline_pt) and \
 629                isinstance(nextnpitem, path.normline_pt) and \
 630                epsilon > math.hypot(thisnpitem.atend_pt()[0] - nextnpitem.atbegin_pt()[0],
 631                                     thisnpitem.atend_pt()[0] - nextnpitem.atbegin_pt()[0]):
 632
 633                 d1,g1,f1,e,f2,g2,d2 = curvescontrols_from_endlines_pt(
 634                     thisnpitem.atend_pt(), tangents[this][-1], tangents[next][0],
 635                     math.hypot(points[this][-1][0] - thisnpitem.atend_pt()[0], points[this][-1][1] - thisnpitem.atend_pt()[1]),
 636                     math.hypot(points[next][0][0] - nextnpitem.atbegin_pt()[0], points[next][0][1] - nextnpitem.atbegin_pt()[1]),
 637                     softness=self.softness)
 638
 639                 if do_moveto:
 640                     smoothpath.append(path.moveto_pt(*d1))
 641                     do_moveto = 0
 642
 643                 smoothpath.append(path.curveto_pt(*(g1 + f1 + e)))
 644                 smoothpath.append(path.curveto_pt(*(f2 + g2 + d2)))
 645
 646             else:
 647
 648                 A, D = points[this][-1], points[next][0]
 649                 tangA, tangD = tangents[this][-1], tangents[next][0]
 650                 curvA, curvD = curvatures[this][-1], curvatures[next][0]
 651                 if not self.obeycurv:
 652                     # do not obey the sign of the curvature but
 653                     # make the sign such that the curve smoothly passes to the next point
 654                     # this results in a discontinuous curvature
 655                     # (but the absolute value is still continuous)
 656                     sA = sign1(tangA[0] * (D[1]-A[1]) - tangA[1] * (D[0]-A[0]))
 657                     sD = sign1(tangD[0] * (D[1]-A[1]) - tangD[1] * (D[0]-A[0]))
 658                     curvA = sA * abs(curvA)
 659                     curvD = sD * abs(curvD)
 660
 661                 # get the length of the control "arms"
 662                 a, d = controldists_from_endpoints_pt (
 663                     A, D, tangA, tangD, curvA, curvD,
 664                     epsilon=epsilon)
 665
 666                 # avoid overshooting at the corners:
 667                 # this changes not only the sign of the curvature
 668                 # but also the magnitude
 669                 if not self.obeycurv:
 670                     t, s = intersection(A, D, tangA, tangD)
 671                     if t is None or t < 0:
 672                         a = None
 673                     else:
 674                         a = min(a, t)
 675
 676                     if s is None or s > 0:
 677                         d = None
 678                     else:
 679                         d = min(d, -s)
 680
 681                 # if there is no useful result:
 682                 # take arbitrary smoothing curve that does not obey
 683                 # the curvature constraints
 684                 if a is None or d is None:
 685                     dist = math.hypot(A[0] - D[0], A[1] - D[1])
 686                     a = dist / (3.0 * math.hypot(*tangA))
 687                     d = dist / (3.0 * math.hypot(*tangD))
 688                     #warnings.warn("The connecting bezier cannot be found. Using a simple fallback.")
 689
 690                 # calculate the two missing control points
 691                 B = A[0] + a * tangA[0], A[1] + a * tangA[1]
 692                 C = D[0] - d * tangD[0], D[1] - d * tangD[1]
 693
 694                 if do_moveto:
 695                     smoothpath.append(path.moveto_pt(*A))
 696                     do_moveto = 0
 697
 698                 smoothpath.append(path.curveto_pt(*(B + C + D)))
 699
 700
 701         # 4. Second part of extra handling of closed paths
 702         if normsubpath.closed:
 703             if do_moveto:
 704                 # XXX the following does not work since dp is not defined
 705                 # probably, we just want a moveto_pt(*normsubpath.atbegin())
 706                 smoothpath.append(path.moveto_pt(*dp.strokepath.atbegin()))
 707                 warnings.warn("The whole subpath has been smoothed away -- sorry")
 708             smoothpath.append(path.closepath())
 709
 710         return smoothpath
 711 # >>>
 712
 713 smoothed.clear = attr.clearclass(smoothed)
 714
 715 class parallel(deformer): # <<<
 716
 717     """creates a parallel path with constant distance to the original path
 718
 719     A positive 'distance' results in a curve left of the original one -- and a
 720     negative 'distance' in a curve at the right. Left/Right are understood in
 721     terms of the parameterization of the original curve.
 722     At corners, either a circular arc is drawn around the corner, or, if the
 723     curve is on the other side, the parallel curve also exhibits a corner.
 724
 725     For each path element a parallel curve/line is constructed. For curves, the
 726     accuracy can be adjusted with the parameter 'relerr', thus, relerr*distance
 727     is the maximum allowable error somewhere in the middle of the curve (at
 728     parameter value 0.5).
 729     'relerr' only applies for the 'expensive' mode where the parallel curve for
 730     a single curve items may be composed of several (many) curve items.
 731     """
 732
 733     # TODO:
 734     # - check for greatest curvature and introduce extra corners
 735     #   if a normcurve is too heavily curved
 736     # - do relerr-checks at better points than just at parameter 0.5
 737
 738     def __init__(self, distance, relerr=0.05, expensive=1):
 739         self.distance = distance
 740         self.relerr = relerr
 741         self.expensive = expensive
 742
 743     def __call__(self, distance=None, relerr=None, expensive=None):
 744         # returns a copy of the deformer with different parameters
 745         if distance is None:
 746             d = self.distance
 747         if relerr is None:
 748             r = self.relerr
 749         if expensive is None:
 750             e = self.expensive
 751
 752         return parallel(distance=d, relerr=r, expensive=e)
 753
 754     def deform(self, basepath):
 755         resultnormsubpaths = [self.deformsubpath(nsp) for nsp in basepath.normpath().normsubpaths]
 756         return path.normpath(resultnormsubpaths)
 757
 758     def deformsubpath(self, orig_nspath):
 759
 760         distance = unit.topt(self.distance)
 761         relerr = self.relerr
 762         expensive = self.expensive
 763         epsilon = orig_nspath.epsilon
 764
 765         new_nspath = path.normsubpath(epsilon=epsilon)
 766
 767         # 1. Store endpoints, tangents and curvatures for each element
 768         points, tangents, curvatures = [], [], []
 769         for npitem in orig_nspath:
 770
 771             ps,ts,cs = [],[],[]
 772             trafos = npitem.trafo([0,1])
 773             for t in trafos:
 774                 p = t._apply(0,0)
 775                 t = t._apply(1,0)
 776                 ps.append(p)
 777                 ts.append((t[0]-p[0], t[1]-p[1]))
 778
 779             rs = npitem.curveradius_pt([0,1])
 780             cs = []
 781             for r in rs:
 782                 if r is None:
 783                     cs.append(0)
 784                 else:
 785                     cs.append(1.0 / r)
 786
 787             points.append(ps)
 788             tangents.append(ts)
 789             curvatures.append(cs)
 790
 791         closeparallel = (tangents[-1][1][0]*tangents[0][0][1] - tangents[-1][1][1]*tangents[0][0][0])
 792
 793         # 2. append the parallel path for each element:
 794         for cur in range(len(orig_nspath)):
 795
 796             if cur == 0:
 797                 old = cur
 798                 # OldEnd = points[old][0]
 799                 OldEndTang = tangents[old][0]
 800             else:
 801                 old = cur - 1
 802                 # OldEnd = points[old][1]
 803                 OldEndTang = tangents[old][1]
 804
 805             CurBeg, CurEnd = points[cur]
 806             CurBegTang, CurEndTang = tangents[cur]
 807             CurBegCurv, CurEndCurv = curvatures[cur]
 808
 809             npitem = orig_nspath[cur]
 810
 811             # get the control points for the shifted pathelement
 812             if isinstance(npitem, path.normline_pt):
 813                 # get the points explicitly from the normal vector
 814                 A = CurBeg[0] - distance * CurBegTang[1], CurBeg[1] + distance * CurBegTang[0]
 815                 D = CurEnd[0] - distance * CurEndTang[1], CurEnd[1] + distance * CurEndTang[0]
 816                 new_npitems = [path.normline_pt(A[0], A[1], D[0], D[1])]
 817             elif isinstance(npitem, path.normcurve_pt):
 818                 # call a function to return a list of controlpoints
 819                 cpoints_list = parallel_curvespoints_pt(npitem, distance, expensive, relerr, epsilon)
 820                 new_npitems = []
 821                 for cpoints in cpoints_list:
 822                     A,B,C,D = cpoints
 823                     new_npitems.append(path.normcurve_pt(A[0],A[1], B[0],B[1], C[0],C[1], D[0],D[1]))
 824                 # we will need the starting point of the new normpath items
 825                 A = cpoints_list[0][0]
 826
 827
 828             # append the next piece of the path:
 829             # it might contain of an extra arc or must be intersected before appending
 830             parallel = (OldEndTang[0]*CurBegTang[1] - OldEndTang[1]*CurBegTang[0])
 831             if parallel*distance < -epsilon:
 832
 833                 # append an arc around the corner
 834                 # from the preceding piece to the current piece
 835                 # we can never get here for the first npitem! (because cur==old)
 836                 endpoint = new_nspath.atend_pt()
 837                 center = CurBeg
 838                 angle1 = math.atan2(endpoint[1] - center[1], endpoint[0] - center[0]) * 180.0 / math.pi
 839                 angle2 = math.atan2(A[1] - center[1], A[0] - center[0]) * 180.0 / math.pi
 840                 if parallel > 0:
 841                     arc_npath = path.path(path.arc_pt(
 842                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 843                 else:
 844                     arc_npath = path.path(path.arcn_pt(
 845                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 846
 847                 for new_npitem in arc_npath[0]:
 848                     new_nspath.append(new_npitem)
 849
 850
 851             elif parallel*distance > epsilon:
 852                 # intersect the extra piece of the path with the rest of the new path
 853                 # and throw away the void parts
 854                 #
 855                 # build a subpath for intersection
 856                 extra_nspath = path.normsubpath(normsubpathitems=new_npitems, epsilon=epsilon)
 857
 858                 intsparams = extra_nspath.intersect(new_nspath)
 859                 # [[a,b,c], [a,b,c]]
 860                 if intsparams:
 861                     # take the first intersection point:
 862                     extra_param, new_param = intsparams[0][0], intsparams[1][0]
 863                     new_nspath = new_nspath.segments([0, new_param])[0]
 864                     extra_nspath = extra_nspath.segments([extra_param, len(extra_nspath)])[0]
 865                     new_npitems = extra_nspath.normsubpathitems
 866                     # in case the intersection was not sufficiently exact:
 867                     # CAREFUL! because we newly created all the new_npitems and
 868                     # the items in extra_nspath, we may in-place change the starting point
 869                     new_npitems[0] = new_npitems[0].modifiedbegin_pt(*new_nspath.atend_pt())
 870                 else:
 871                     raise # how did we get here?
 872
 873
 874             # at the (possible) closing corner we may have to intersect another time
 875             # or add another corner:
 876             # the intersection must be done before appending the parallel piece
 877             if orig_nspath.closed and cur == len(orig_nspath) - 1:
 878                 if closeparallel * distance > epsilon:
 879                     intsparams = extra_nspath.intersect(new_nspath)
 880                     # [[a,b,c], [a,b,c]]
 881                     if intsparams:
 882                         # take the last intersection point:
 883                         extra_param, new_param = intsparams[0][-1], intsparams[1][-1]
 884                         new_nspath = new_nspath.segments([new_param, len(new_nspath)])[0]
 885                         extra_nspath = extra_nspath.segments([0, extra_param])[0]
 886                         new_npitems = extra_nspath.normsubpathitems
 887                         # in case the intersection was not sufficiently exact:
 888                         # CAREFUL! because we newly created all the new_npitems and
 889                         # the items in extra_nspath, we may in-place change the end point
 890                         new_npitems[-1] = new_npitems[-1].modifiedend_pt(*new_nspath.atbegin_pt())
 891                     else:
 892                         raise # how did we get here?
 893
 894                     pass
 895
 896
 897             # append the parallel piece
 898             for new_npitem in new_npitems:
 899                 new_nspath.append(new_npitem)
 900
 901
 902         # the curve around the closing corner must be added at last:
 903         if orig_nspath.closed:
 904             if closeparallel * distance < -epsilon:
 905                 endpoint = new_nspath.atend_pt()
 906                 center = orig_nspath.atend_pt()
 907                 A = new_nspath.atbegin_pt()
 908                 angle1 = math.atan2(endpoint[1] - center[1], endpoint[0] - center[0]) * 180.0 / math.pi
 909                 angle2 = math.atan2(A[1] - center[1], A[0] - center[0]) * 180.0 / math.pi
 910                 if parallel > 0:
 911                     arc_npath = path.path(path.arc_pt(
 912                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 913                 else:
 914                     arc_npath = path.path(path.arcn_pt(
 915                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 916
 917                 for new_npitem in arc_npath[0]:
 918                     new_nspath.append(new_npitem)
 919
 920         # 3. extra handling of closed paths
 921         if orig_nspath.closed:
 922             new_nspath.close()
 923
 924         return new_nspath
 925 # >>>
 926
 927 parallel.clear = attr.clearclass(parallel)
 928
 929 # vim:foldmethod=marker:foldmarker=<<<,>>>