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 origpath
 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 class smoothed(deformer): # <<<
 486
 487     """Bends corners in a path.
 488
 489     This decorator replaces corners in a path with bezier curves. There are two cases:
 490     - If the corner lies between two lines, _two_ bezier curves will be used
 491       that are highly optimized to look good (their curvature is to be zero at the ends
 492       and has to have zero derivative in the middle).
 493       Additionally, it can controlled by the softness-parameter.
 494     - If the corner lies between curves then _one_ bezier is used that is (except in some
 495       special cases) uniquely determined by the tangents and curvatures at its end-points.
 496       In some cases it is necessary to use only the absolute value of the curvature to avoid a
 497       cusp-shaped connection of the new bezier to the old path. In this case the use of
 498       "obeycurv=0" allows the sign of the curvature to switch.
 499     - The radius argument gives the arclength-distance of the corner to the points where the
 500       old path is cut and the beziers are inserted.
 501     - Path elements that are too short (shorter than the radius) are skipped
 502     """
 503
 504     def __init__(self, radius, softness=1, obeycurv=0):
 505         self.radius = radius
 506         self.softness = softness
 507         self.obeycurv = obeycurv
 508
 509     def __call__(self, radius=None, softness=None, obeycurv=None):
 510         if radius is None:
 511             radius = self.radius
 512         if softness is None:
 513             softness = self.softness
 514         if obeycurv is None:
 515             obeycurv = self.obeycurv
 516         return smoothed(radius=radius, softness=softness, obeycurv=obeycurv)
 517
 518     def deform(self, basepath):
 519         basepath = basepath.normpath()
 520         smoothpath = path.path()
 521
 522         for sp in basepath.normsubpaths:
 523             smoothpath += self.deformsubpath(sp)
 524
 525         return smoothpath
 526
 527     def deformsubpath(self, normsubpath):
 528
 529         radius = unit.topt(self.radius)
 530         epsilon = normsubpath.epsilon
 531
 532         # 1. Build up a list of all relevant normsubpathitems
 533         #    and the lengths where they will be cut (length with respect to the normsubpath)
 534         all_npitems = normsubpath.normsubpathitems
 535         rel_npitems, arclengths = [], []
 536         for npitem in all_npitems:
 537
 538             arclength = npitem.arclen_pt(epsilon)
 539
 540             # items that should be totally skipped:
 541             # (we throw away the possible ending "closepath" piece)
 542             if (arclength > radius):
 543                 rel_npitems.append(npitem)
 544                 arclengths.append(arclength)
 545             else:
 546                 warnings.warn("smoothed is skipping a too short normsubpathitem")
 547
 548         # 2. Find the parameters, points,
 549         #    and calculate tangents and curvatures
 550         params, points, tangents, curvatures = [], [], [], []
 551         for npitem, arclength in zip(rel_npitems, arclengths):
 552
 553             # find the parameter(s): either one or two
 554             # for short items we squeeze the radius
 555             if arclength > 2 * radius:
 556                 cut_alengths = [radius, arclength - radius]
 557             else:
 558                 cut_alengths = [0.5 * radius]
 559
 560             # get the parameters
 561             pars = npitem._arclentoparam_pt(cut_alengths, epsilon)[0]
 562
 563             # the endpoints of an open subpath must be handled specially
 564             if not normsubpath.closed:
 565                 if npitem is rel_npitems[0]:
 566                     pars[0] = 0
 567                 if npitem is rel_npitems[-1]:
 568                     pars[-1] = 1
 569
 570             # find points, tangents and curvatures
 571             ts,cs,ps = [],[],[]
 572             for par in pars:
 573                 thetrafo = npitem.trafo([par])[0]
 574                 p = thetrafo._apply(0,0)
 575                 t = thetrafo._apply(1,0)
 576                 ps.append(p)
 577                 ts.append((t[0]-p[0], t[1]-p[1]))
 578                 r = npitem.curveradius_pt([par])[0]
 579                 if r is None:
 580                     cs.append(0)
 581                 else:
 582                     cs.append(1.0 / r)
 583
 584             params.append(pars)
 585             points.append(ps)
 586             tangents.append(ts)
 587             curvatures.append(cs)
 588
 589
 590         # create an empty path to collect pathitems
 591         # this will be returned as normpath, later
 592         smoothpath = path.path()
 593         do_moveto = 1 # we do not know yet where to moveto
 594
 595         # 3. Do the splitting for the first to the last element,
 596         #    a closed path must be closed later
 597         #
 598         for i in range(len(rel_npitems)):
 599
 600             this = i
 601             next = (i+1) % len(rel_npitems)
 602             thisnpitem = rel_npitems[this]
 603             nextnpitem = rel_npitems[next]
 604
 605             # split thisnpitem apart and take the middle piece
 606             # We start the new path with the middle piece of the first path-elem
 607             if len(points[this]) == 2:
 608
 609                 middlepiece = thisnpitem.segments(params[this])[0]
 610
 611                 if do_moveto:
 612                     smoothpath.append(path.moveto_pt(*middlepiece.atbegin_pt()))
 613                     do_moveto = 0
 614
 615                 if isinstance(middlepiece, path.normline_pt):
 616                     smoothpath.append(path.lineto_pt(*middlepiece.atend_pt()))
 617                 elif isinstance(middlepiece, path.normcurve_pt):
 618                     smoothpath.append(path.curveto_pt(
 619                       middlepiece.x1_pt, middlepiece.y1_pt,
 620                       middlepiece.x2_pt, middlepiece.y2_pt,
 621                       middlepiece.x3_pt, middlepiece.y3_pt))
 622
 623             if (not normsubpath.closed) and (thisnpitem is rel_npitems[-1]):
 624                 continue
 625
 626             # add the curve(s) replacing the corner
 627             if isinstance(thisnpitem, path.normline_pt) and \
 628                isinstance(nextnpitem, path.normline_pt) and \
 629                epsilon > math.hypot(thisnpitem.atend_pt()[0] - nextnpitem.atbegin_pt()[0],
 630                                     thisnpitem.atend_pt()[0] - nextnpitem.atbegin_pt()[0]):
 631
 632                 d1,g1,f1,e,f2,g2,d2 = curvescontrols_from_endlines_pt(
 633                     thisnpitem.atend_pt(), tangents[this][-1], tangents[next][0],
 634                     math.hypot(points[this][-1][0] - thisnpitem.atend_pt()[0], points[this][-1][1] - thisnpitem.atend_pt()[1]),
 635                     math.hypot(points[next][0][0] - nextnpitem.atbegin_pt()[0], points[next][0][1] - nextnpitem.atbegin_pt()[1]),
 636                     softness=self.softness)
 637
 638                 if do_moveto:
 639                     smoothpath.append(path.moveto_pt(*d1))
 640                     do_moveto = 0
 641
 642                 smoothpath.append(path.curveto_pt(*(g1 + f1 + e)))
 643                 smoothpath.append(path.curveto_pt(*(f2 + g2 + d2)))
 644
 645             else:
 646
 647                 A, D = points[this][-1], points[next][0]
 648                 tangA, tangD = tangents[this][-1], tangents[next][0]
 649                 curvA, curvD = curvatures[this][-1], curvatures[next][0]
 650                 if not self.obeycurv:
 651                     # do not obey the sign of the curvature but
 652                     # make the sign such that the curve smoothly passes to the next point
 653                     # this results in a discontinuous curvature
 654                     # (but the absolute value is still continuous)
 655                     sA = sign1(tangA[0] * (D[1]-A[1]) - tangA[1] * (D[0]-A[0]))
 656                     sD = sign1(tangD[0] * (D[1]-A[1]) - tangD[1] * (D[0]-A[0]))
 657                     curvA = sA * abs(curvA)
 658                     curvD = sD * abs(curvD)
 659
 660                 # get the length of the control "arms"
 661                 a, d = controldists_from_endpoints_pt (
 662                     A, D, tangA, tangD, curvA, curvD,
 663                     epsilon=epsilon)
 664
 665                 # avoid overshooting at the corners:
 666                 # this changes not only the sign of the curvature
 667                 # but also the magnitude
 668                 if not self.obeycurv:
 669                     t, s = intersection(A, D, tangA, tangD)
 670                     if t is None or t < 0:
 671                         a = None
 672                     else:
 673                         a = min(a, t)
 674
 675                     if s is None or s > 0:
 676                         d = None
 677                     else:
 678                         d = min(d, -s)
 679
 680                 # if there is no useful result:
 681                 # take arbitrary smoothing curve that does not obey
 682                 # the curvature constraints
 683                 if a is None or d is None:
 684                     dist = math.hypot(A[0] - D[0], A[1] - D[1])
 685                     a = dist / (3.0 * math.hypot(*tangA))
 686                     d = dist / (3.0 * math.hypot(*tangD))
 687                     #warnings.warn("The connecting bezier cannot be found. Using a simple fallback.")
 688
 689                 # calculate the two missing control points
 690                 B = A[0] + a * tangA[0], A[1] + a * tangA[1]
 691                 C = D[0] - d * tangD[0], D[1] - d * tangD[1]
 692
 693                 if do_moveto:
 694                     smoothpath.append(path.moveto_pt(*A))
 695                     do_moveto = 0
 696
 697                 smoothpath.append(path.curveto_pt(*(B + C + D)))
 698
 699
 700         # 4. Second part of extra handling of closed paths
 701         if normsubpath.closed:
 702             if do_moveto:
 703                 smoothpath.append(path.moveto_pt(*dp.strokepath.atbegin()))
 704                 warnings.warn("The whole subpath has been smoothed away -- sorry")
 705             smoothpath.append(path.closepath())
 706
 707         return smoothpath
 708 # >>>
 709
 710 smoothed.clear = attr.clearclass(smoothed)
 711
 712 _base = unit.v_cm
 713 smoothed.SHARP = smoothed(radius=_base/math.sqrt(64))
 714 smoothed.SHARp = smoothed(radius=_base/math.sqrt(32))
 715 smoothed.SHArp = smoothed(radius=_base/math.sqrt(16))
 716 smoothed.SHarp = smoothed(radius=_base/math.sqrt(8))
 717 smoothed.Sharp = smoothed(radius=_base/math.sqrt(4))
 718 smoothed.sharp = smoothed(radius=_base/math.sqrt(2))
 719 smoothed.normal = smoothed(radius=_base)
 720 smoothed.round = smoothed(radius=_base*math.sqrt(2))
 721 smoothed.Round = smoothed(radius=_base*math.sqrt(4))
 722 smoothed.ROund = smoothed(radius=_base*math.sqrt(8))
 723 smoothed.ROUnd = smoothed(radius=_base*math.sqrt(16))
 724 smoothed.ROUNd = smoothed(radius=_base*math.sqrt(32))
 725 smoothed.ROUND = smoothed(radius=_base*math.sqrt(64))
 726
 727 class parallel(deformer): # <<<
 728
 729     """creates a parallel path with constant distance to the original path
 730
 731     A positive 'distance' results in a curve left of the original one -- and a
 732     negative 'distance' in a curve at the right. Left/Right are understood in
 733     terms of the parameterization of the original curve.
 734     At corners, either a circular arc is drawn around the corner, or, if the
 735     curve is on the other side, the parallel curve also exhibits a corner.
 736
 737     For each path element a parallel curve/line is constructed. For curves, the
 738     accuracy can be adjusted with the parameter 'relerr', thus, relerr*distance
 739     is the maximum allowable error somewhere in the middle of the curve (at
 740     parameter value 0.5).
 741     'relerr' only applies for the 'expensive' mode where the parallel curve for
 742     a single curve items may be composed of several (many) curve items.
 743     """
 744
 745     # TODO:
 746     # - check for greatest curvature and introduce extra corners
 747     #   if a normcurve is too heavily curved
 748     # - do relerr-checks at better points than just at parameter 0.5
 749
 750     def __init__(self, distance, relerr=0.05, expensive=1):
 751         self.distance = distance
 752         self.relerr = relerr
 753         self.expensive = expensive
 754
 755     def __call__(self, distance=None, relerr=None, expensive=None):
 756         # returns a copy of the deformer with different parameters
 757         if distance is None:
 758             d = self.distance
 759         if relerr is None:
 760             r = self.relerr
 761         if expensive is None:
 762             e = self.expensive
 763
 764         return parallel(distance=d, relerr=r, expensive=e)
 765
 766     def deform(self, basepath):
 767         resultnormsubpaths = [self.deformsubpath(nsp) for nsp in basepath.normpath().normsubpaths]
 768         return path.normpath(resultnormsubpaths)
 769
 770     def deformsubpath(self, orig_nspath):
 771
 772         distance = unit.topt(self.distance)
 773         relerr = self.relerr
 774         expensive = self.expensive
 775         epsilon = orig_nspath.epsilon
 776
 777         new_nspath = path.normsubpath(epsilon=epsilon)
 778
 779         # 1. Store endpoints, tangents and curvatures for each element
 780         points, tangents, curvatures = [], [], []
 781         for npitem in orig_nspath:
 782
 783             ps,ts,cs = [],[],[]
 784             trafos = npitem.trafo([0,1])
 785             for t in trafos:
 786                 p = t._apply(0,0)
 787                 t = t._apply(1,0)
 788                 ps.append(p)
 789                 ts.append((t[0]-p[0], t[1]-p[1]))
 790
 791             rs = npitem.curveradius_pt([0,1])
 792             cs = []
 793             for r in rs:
 794                 if r is None:
 795                     cs.append(0)
 796                 else:
 797                     cs.append(1.0 / r)
 798
 799             points.append(ps)
 800             tangents.append(ts)
 801             curvatures.append(cs)
 802
 803         closeparallel = (tangents[-1][1][0]*tangents[0][0][1] - tangents[-1][1][1]*tangents[0][0][0])
 804
 805         # 2. append the parallel path for each element:
 806         for cur in range(len(orig_nspath)):
 807
 808             if cur == 0:
 809                 old = cur
 810                 OldEnd = points[old][0]
 811                 OldEndTang = tangents[old][0]
 812             else:
 813                 old = cur - 1
 814                 OldEnd = points[old][1]
 815                 OldEndTang = tangents[old][1]
 816
 817             CurBeg, CurEnd = points[cur]
 818             CurBegTang, CurEndTang = tangents[cur]
 819             CurBegCurv, CurEndCurv = curvatures[cur]
 820
 821             npitem = orig_nspath[cur]
 822
 823             # get the control points for the shifted pathelement
 824             if isinstance(npitem, path.normline_pt):
 825                 # get the points explicitly from the normal vector
 826                 A = CurBeg[0] - distance * CurBegTang[1], CurBeg[1] + distance * CurBegTang[0]
 827                 D = CurEnd[0] - distance * CurEndTang[1], CurEnd[1] + distance * CurEndTang[0]
 828                 new_npitems = [path.normline_pt(A[0], A[1], D[0], D[1])]
 829             elif isinstance(npitem, path.normcurve_pt):
 830                 # call a function to return a list of controlpoints
 831                 cpoints_list = parallel_curvespoints_pt(npitem, distance, expensive, relerr, epsilon)
 832                 new_npitems = []
 833                 for cpoints in cpoints_list:
 834                     A,B,C,D = cpoints
 835                     new_npitems.append(path.normcurve_pt(A[0],A[1], B[0],B[1], C[0],C[1], D[0],D[1]))
 836                 # we will need the starting point of the new normpath items
 837                 A = cpoints_list[0][0]
 838
 839
 840             # append the next piece of the path:
 841             # it might contain of an extra arc or must be intersected before appending
 842             parallel = (OldEndTang[0]*CurBegTang[1] - OldEndTang[1]*CurBegTang[0])
 843             if parallel*distance < -epsilon:
 844
 845                 # append an arc around the corner
 846                 # from the preceding piece to the current piece
 847                 # we can never get here for the first npitem! (because cur==old)
 848                 endpoint = new_nspath.atend_pt()
 849                 center = CurBeg
 850                 angle1 = math.atan2(endpoint[1] - center[1], endpoint[0] - center[0]) * 180.0 / math.pi
 851                 angle2 = math.atan2(A[1] - center[1], A[0] - center[0]) * 180.0 / math.pi
 852                 if parallel > 0:
 853                     arc_npath = path.path(path.arc_pt(
 854                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 855                 else:
 856                     arc_npath = path.path(path.arcn_pt(
 857                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 858
 859                 for new_npitem in arc_npath[0]:
 860                     new_nspath.append(new_npitem)
 861
 862
 863             elif parallel*distance > epsilon:
 864                 # intersect the extra piece of the path with the rest of the new path
 865                 # and throw away the void parts
 866                 #
 867                 # build a subpath for intersection
 868                 extra_nspath = path.normsubpath(normsubpathitems=new_npitems, epsilon=epsilon)
 869
 870                 intsparams = extra_nspath.intersect(new_nspath)
 871                 # [[a,b,c], [a,b,c]]
 872                 if intsparams:
 873                     # take the first intersection point:
 874                     extra_param, new_param = intsparams[0][0], intsparams[1][0]
 875                     new_nspath = new_nspath.segments([0, new_param])[0]
 876                     extra_nspath = extra_nspath.segments([extra_param, len(extra_nspath)])[0]
 877                     new_npitems = extra_nspath[:]
 878                     # in case the intersection was not sufficiently exact:
 879                     new_npitems[0].x0_pt, new_npitems[0].y0_pt = new_nspath.atend_pt()
 880                 else:
 881                     raise # how did we get here?
 882
 883
 884             # at the (possible) closing corner we may have to intersect another time
 885             # or add another corner:
 886             # the intersection must be done before appending the parallel piece
 887             if orig_nspath.closed and cur == len(orig_nspath) - 1:
 888                 if closeparallel * distance > epsilon:
 889                     intsparams = extra_nspath.intersect(new_nspath)
 890                     # [[a,b,c], [a,b,c]]
 891                     if intsparams:
 892                         # take the first intersection point:
 893                         extra_param, new_param = intsparams[0][-1], intsparams[1][-1]
 894                         new_nspath = new_nspath.segments([new_param, len(new_nspath)])[0]
 895                         extra_nspath = extra_nspath.segments([0, extra_param])[0]
 896                         new_npitems = extra_nspath[:]
 897                         # in case the intersection was not sufficiently exact:
 898                         if isinstance(new_npitems[0], path.normcurve_pt):
 899                             new_npitems[0].x3_pt, new_npitems[0].y3_pt = new_nspath.atend_pt()
 900                         elif isinstance(new_npitems[0], path.normline_pt):
 901                             new_npitems[0].x1_pt, new_npitems[0].y1_pt = new_nspath.atend_pt()
 902                     else:
 903                         raise # how did we get here?
 904
 905                     pass
 906
 907
 908             # append the parallel piece
 909             for new_npitem in new_npitems:
 910                 new_nspath.append(new_npitem)
 911
 912
 913         # the curve around the closing corner must be added at last:
 914         if orig_nspath.closed:
 915             if closeparallel * distance < -epsilon:
 916                 endpoint = new_nspath.atend_pt()
 917                 center = orig_nspath.atend_pt()
 918                 A = new_nspath.atbegin_pt()
 919                 angle1 = math.atan2(endpoint[1] - center[1], endpoint[0] - center[0]) * 180.0 / math.pi
 920                 angle2 = math.atan2(A[1] - center[1], A[0] - center[0]) * 180.0 / math.pi
 921                 if parallel > 0:
 922                     arc_npath = path.path(path.arc_pt(
 923                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 924                 else:
 925                     arc_npath = path.path(path.arcn_pt(
 926                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 927
 928                 for new_npitem in arc_npath[0]:
 929                     new_nspath.append(new_npitem)
 930
 931         # 3. extra handling of closed paths
 932         if orig_nspath.closed:
 933             new_nspath.close()
 934
 935         return new_nspath
 936 # >>>
 937
 938 # vim:foldmethod=marker:foldmarker=<<<,>>>