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
  28 def sign1(x):
  29     return (x >= 0) and 1 or -1
  30
  31 def curvescontrols_from_endlines_pt (B, tangent1, tangent2, r1, r2, softness): # <<<
  32     # calculates the parameters for two bezier curves connecting two lines (curvature=0)
  33     # starting at B - r1*tangent1
  34     # ending at   B + r2*tangent2
  35     #
  36     # Takes the corner B
  37     # and two tangent vectors heading to and from B
  38     # and two radii r1 and r2:
  39     # All arguments must be in Points
  40     # Returns the seven control points of the two bezier curves:
  41     #  - start d1
  42     #  - control points g1 and f1
  43     #  - midpoint e
  44     #  - control points f2 and g2
  45     #  - endpoint d2
  46
  47     # make direction vectors d1: from B to A
  48     #                        d2: from B to C
  49     d1 = -tangent1[0] / math.hypot(*tangent1), -tangent1[1] / math.hypot(*tangent1)
  50     d2 =  tangent2[0] / math.hypot(*tangent2),  tangent2[1] / math.hypot(*tangent2)
  51
  52     # 0.3192 has turned out to be the maximum softness available
  53     # for straight lines ;-)
  54     f = 0.3192 * softness
  55     g = (15.0 * f + math.sqrt(-15.0*f*f + 24.0*f))/12.0
  56
  57     # make the control points of the two bezier curves
  58     f1 = B[0] + f * r1 * d1[0], B[1] + f * r1 * d1[1]
  59     f2 = B[0] + f * r2 * d2[0], B[1] + f * r2 * d2[1]
  60     g1 = B[0] + g * r1 * d1[0], B[1] + g * r1 * d1[1]
  61     g2 = B[0] + g * r2 * d2[0], B[1] + g * r2 * d2[1]
  62     d1 = B[0] +     r1 * d1[0], B[1] +     r1 * d1[1]
  63     d2 = B[0] +     r2 * d2[0], B[1] +     r2 * d2[1]
  64     e  = 0.5 * (f1[0] + f2[0]), 0.5 * (f1[1] + f2[1])
  65
  66     return (d1, g1, f1, e, f2, g2, d2)
  67 # >>>
  68
  69 def controldists_from_endpoints_pt (A, B, tangA, tangB, curvA, curvB): # <<<
  70
  71     """distances for a curve given by tangents and curvatures at the endpoints
  72
  73     This helper routine returns the two distances between the endpoints and the
  74     corresponding control points of a (cubic) bezier curve that has
  75     prescribed tangents tangentA, tangentB and curvatures curvA, curvB at the
  76     end points.
  77     """
  78
  79     # some shortcuts
  80     T = tangB[0] * tangA[1] - tangB[1] * tangA[0]
  81     D =   tangA[0] * (B[1]-A[1]) - tangA[1] * (B[0]-A[0])
  82     E = - tangB[0] * (B[1]-A[1]) + tangB[1] * (B[0]-A[0])
  83     # the variables: \dot X(0) = 3 * a * tangA
  84     #                \dot X(1) = 3 * b * tangB
  85     a, b = None, None
  86
  87
  88     # try some special cases where the equations decouple
  89     try:
  90         1.0 / T
  91     except ZeroDivisionError:
  92         try:
  93             a = math.sqrt(2.0 * D / (3.0 * curvA))
  94             b = math.sqrt(2.0 * E / (3.0 * curvB))
  95         except ZeroDivisionError:
  96             a = b = None
  97         except ValueError:
  98             raise # ???
  99     else:
 100         try:
 101             1.0 / curvA
 102         except ZeroDivisionError:
 103             b = -D / T
 104             a = (1.5*curvB*b*b - E) / T
 105         else:
 106             try:
 107                 1.0 / curvB
 108             except ZeroDivisionError:
 109                 a = -E / T
 110                 b = (1.5*curvA*a*a - D) / T
 111             else:
 112                 a, b = None, None
 113
 114
 115     # else find a solution for the full problem
 116     if a is None:
 117         # we first try to find all the zeros of the polynomials for a or b (4th order)
 118         # this needs Numeric and LinearAlgebra
 119
 120         #     0 = Ga(a,b) = 0.5 a |a| curvA + b * T - D
 121         #     0 = Gb(a,b) = 0.5 b |b| curvB + a * T + E
 122
 123         coeffs_a = (3.375*curvA*curvA*curvB, 0, -4.5*curvA*curvB*D, -T**3,  1.5*curvB*D*D - T*T*E)
 124         coeffs_b = (3.375*curvA*curvB*curvB, 0, -4.5*curvA*curvB*E, -T**3,  1.5*curvA*E*E - T*T*D)
 125
 126         # First try the equation for a
 127         cands_a = [cand for cand in helper.realpolyroots(coeffs_a) if cand > 0]
 128
 129         if cands_a:
 130             a = min(cands_a)
 131             b = (1.5*curvA*a*a - D) / T
 132         else:
 133             # then, try the equation for b
 134             cands_b = [cand for cand in helper.realpolyroots(coeffs_b) if cand > 0]
 135             if cands_b:
 136                 b = min(cands_b)
 137                 a = (1.5*curvB*b*b - E) / T
 138             else:
 139                 a = b = None
 140
 141     if a < 0 or b < 0:
 142         a = b = None
 143
 144     return a, b
 145 # >>>
 146
 147 def intersection (A, D, tangA, tangD): # <<<
 148
 149     """returns the intersection parameters of two evens
 150
 151     they are defined by:
 152       x(t) = A + t * tangA
 153       x(s) = D + s * tangD
 154     """
 155     det = -tangA[0] * tangD[1] + tangA[1] * tangD[0]
 156     try:
 157         1.0 / det
 158     except ArithmeticError:
 159         return None, None
 160
 161     DA = D[0] - A[0], D[1] - A[1]
 162
 163     t = (-tangD[1]*DA[0] + tangD[0]*DA[1]) / det
 164     s = (-tangA[1]*DA[0] + tangA[0]*DA[1]) / det
 165
 166     return t, s
 167 # >>>
 168
 169 def parallel_curvespoints_pt (orig_ncurve, shift, expensive=0, relerr=0.05, epsilon=1e-5, counter=1): # <<<
 170
 171     A = orig_ncurve.x0_pt, orig_ncurve.y0_pt
 172     B = orig_ncurve.x1_pt, orig_ncurve.y1_pt
 173     C = orig_ncurve.x2_pt, orig_ncurve.y2_pt
 174     D = orig_ncurve.x3_pt, orig_ncurve.y3_pt
 175
 176     # non-normalized tangential vector
 177     # from begin/end point to the corresponding controlpoint
 178     tangA = (B[0] - A[0], B[1] - A[1])
 179     tangD = (D[0] - C[0], D[1] - C[1])
 180
 181     # normalized normal vectors
 182     # turned to the left (+90 degrees) from the tangents
 183     NormA = (-tangA[1] / math.hypot(*tangA), tangA[0] / math.hypot(*tangA))
 184     NormD = (-tangD[1] / math.hypot(*tangD), tangD[0] / math.hypot(*tangD))
 185
 186     # radii of curvature
 187     radiusA, radiusD = orig_ncurve.curveradius_pt([0,1])
 188
 189     # get the new begin/end points
 190     A = A[0] + shift * NormA[0], A[1] + shift * NormA[1]
 191     D = D[0] + shift * NormD[0], D[1] + shift * NormD[1]
 192
 193     try:
 194         if radiusA is None:
 195             curvA = 0
 196         else:
 197             curvA = 1.0 / (radiusA - shift)
 198         if radiusD is None:
 199             curvD = 0
 200         else:
 201             curvD = 1.0 / (radiusD - shift)
 202     except ZeroDivisionError:
 203         raise
 204     else:
 205         a, d = controldists_from_endpoints_pt (A, D, tangA, tangD, curvA, curvD)
 206
 207         if a is None or d is None:
 208             # fallback heuristic
 209             a = (radiusA - shift) / radiusA
 210             d = (radiusD - shift) / radiusD
 211
 212         B = A[0] + a * tangA[0], A[1] + a * tangA[1]
 213         C = D[0] - d * tangD[0], D[1] - d * tangD[1]
 214
 215     controlpoints = [(A,B,C,D)]
 216
 217     # check if the distance is really the wanted distance
 218     if expensive and counter < 10:
 219         # measure the distance in the "middle" of the original curve
 220         trafo = orig_ncurve.trafo([0.5])[0]
 221         M = trafo.apply_pt(0,0)
 222         NormM = trafo.apply_pt(0,1)
 223         NormM = NormM[0] - M[0], NormM[1] - M[1]
 224
 225         nline = path.normline_pt (
 226           M[0] + (1.0 - 2*relerr) * shift * NormM[0],
 227           M[1] + (1.0 - 2*relerr) * shift * NormM[1],
 228           M[0] + (1.0 + 2*relerr) * shift * NormM[0],
 229           M[1] + (1.0 + 2*relerr) * shift * NormM[1])
 230
 231         new_ncurve = path.normcurve_pt(A[0],A[1], B[0],B[1], C[0],C[1], D[0],D[1])
 232
 233         #cutparams = nline.intersect(orig_ncurve, epsilon)
 234         cutparams = new_ncurve.intersect(nline, epsilon)
 235         if cutparams:
 236             cutpoints = nline.at_pt(cutparams[0])
 237         else:
 238             cutpoints = []
 239         good = 0
 240         for cutpoint in cutpoints:
 241             if cutpoint is not None:
 242                 dist = math.hypot(M[0] - cutpoint[0], M[1] - cutpoint[1])
 243                 if abs(dist - abs(shift)) < relerr * abs(shift):
 244                     good = 1
 245
 246         if not good:
 247             first, second = orig_ncurve.segments([0,0.5,1])
 248             controlpoints = \
 249               parallel_curvespoints_pt (first,  shift, expensive, relerr, epsilon, counter+1) + \
 250               parallel_curvespoints_pt (second, shift, expensive, relerr, epsilon, counter+1)
 251
 252
 253
 254     # TODO:
 255     # too big curvatures: intersect curves
 256     # there is something wrong with the recursion
 257     return controlpoints
 258 # >>>
 259
 260
 261 class deformer(attr.attr):
 262
 263     def deform (self, basepath):
 264         return basepath
 265
 266 class cycloid(deformer): # <<<
 267     """Wraps a cycloid around a path.
 268
 269     The outcome looks like a metal spring with the originalpath as the axis.
 270     radius: radius of the cycloid
 271     loops:  number of loops from beginning to end of the original path
 272     skipfirst/skiplast: undeformed end lines of the original path
 273
 274     """
 275
 276     def __init__(self, radius=0.5*unit.t_cm, halfloops=10,
 277     skipfirst=1*unit.t_cm, skiplast=1*unit.t_cm, curvesperhloop=3, sign=1, turnangle=45):
 278         self.skipfirst = skipfirst
 279         self.skiplast = skiplast
 280         self.radius = radius
 281         self.halfloops = halfloops
 282         self.curvesperhloop = curvesperhloop
 283         self.sign = sign
 284         self.turnangle = turnangle
 285
 286     def __call__(self, radius=None, halfloops=None,
 287     skipfirst=None, skiplast=None, curvesperhloop=None, sign=None, turnangle=None):
 288         if radius is None:
 289             radius = self.radius
 290         if halfloops is None:
 291             halfloops = self.halfloops
 292         if skipfirst is None:
 293             skipfirst = self.skipfirst
 294         if skiplast is None:
 295             skiplast = self.skiplast
 296         if curvesperhloop is None:
 297             curvesperhloop = self.curvesperhloop
 298         if sign is None:
 299             sign = self.sign
 300         if turnangle is None:
 301             turnangle = self.turnangle
 302
 303         return cycloid(radius=radius, halfloops=halfloops, skipfirst=skipfirst, skiplast=skiplast,
 304                        curvesperhloop=curvesperhloop, sign=sign, turnangle=turnangle)
 305
 306     def deform(self, basepath):
 307         resultnormsubpaths = [self.deformsubpath(nsp) for nsp in basepath.normpath().normsubpaths]
 308         return path.normpath(resultnormsubpaths)
 309
 310     def deformsubpath(self, normsubpath):
 311
 312         skipfirst = abs(unit.topt(self.skipfirst))
 313         skiplast = abs(unit.topt(self.skiplast))
 314         radius = abs(unit.topt(self.radius))
 315         turnangle = self.turnangle * math.pi / 180.0
 316
 317         cosTurn = math.cos(turnangle)
 318         sinTurn = math.sin(turnangle)
 319
 320         # make list of the lengths and parameters at points on normsubpath where we will add cycloid-points
 321         totlength = normsubpath.arclen_pt()
 322         if totlength <= skipfirst + skiplast + 2*radius*sinTurn:
 323             warnings.warn("normsubpath is too short for deformation with cycloid -- skipping...")
 324             return normsubpath
 325
 326         # parametrisation is in rotation-angle around the basepath
 327         # differences in length, angle ... between two basepoints
 328         # and between basepoints and controlpoints
 329         Dphi = math.pi / self.curvesperhloop
 330         phis = [i * Dphi for i in range(self.halfloops * self.curvesperhloop + 1)]
 331         DzDphi = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * math.pi * cosTurn)
 332         # Dz = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * self.curvesperhloop * cosTurn)
 333         # zs = [i * Dz for i in range(self.halfloops * self.curvesperhloop + 1)]
 334         # from path._arctobcurve:
 335         # optimal relative distance along tangent for second and third control point
 336         L = 4 * radius * (1 - math.cos(Dphi/2)) / (3 * math.sin(Dphi/2))
 337
 338         # Now the transformation of z into the turned coordinate system
 339         Zs = [ skipfirst + radius*sinTurn # here the coordinate z starts
 340              - sinTurn*radius*math.cos(phi) + cosTurn*DzDphi*phi # the transformed z-coordinate
 341              for phi in phis]
 342         params = normsubpath._arclentoparam_pt(Zs)[0]
 343
 344         # get the positions of the splitpoints in the cycloid
 345         points = []
 346         for phi, param in zip(phis, params):
 347             # the cycloid is a circle that is stretched along the normsubpath
 348             # here are the points of that circle
 349             basetrafo = normsubpath.trafo([param])[0]
 350
 351             # The point on the cycloid, in the basepath's local coordinate system
 352             baseZ, baseY = 0, radius*math.sin(phi)
 353
 354             # The tangent there, also in local coords
 355             tangentX = -cosTurn*radius*math.sin(phi) + sinTurn*DzDphi
 356             tangentY = radius*math.cos(phi)
 357             tangentZ = sinTurn*radius*math.sin(phi) + DzDphi*cosTurn
 358             norm = math.sqrt(tangentX*tangentX + tangentY*tangentY + tangentZ*tangentZ)
 359             tangentY, tangentZ = tangentY/norm, tangentZ/norm
 360
 361             # Respect the curvature of the basepath for the cycloid's curvature
 362             # XXX this is only a heuristic, not a "true" expression for
 363             #     the curvature in curved coordinate systems
 364             pathradius = normsubpath.curveradius_pt([param])[0]
 365             if pathradius is not None:
 366                 factor = (pathradius - baseY) / pathradius
 367                 factor = abs(factor)
 368             else:
 369                 factor = 1
 370             l = L * factor
 371
 372             # The control points prior and after the point on the cycloid
 373             preeZ, preeY = baseZ - l * tangentZ, baseY - l * tangentY
 374             postZ, postY = baseZ + l * tangentZ, baseY + l * tangentY
 375
 376             # Now put everything at the proper place
 377             points.append(basetrafo.apply_pt(preeZ, self.sign * preeY) +
 378                           basetrafo.apply_pt(baseZ, self.sign * baseY) +
 379                           basetrafo.apply_pt(postZ, self.sign * postY))
 380
 381         if len(points) <= 1:
 382             warnings.warn("normsubpath is too short for deformation with cycloid -- skipping...")
 383             return normsubpath
 384
 385         # Build the path from the pointlist
 386         # containing (control x 2,  base x 2, control x 2)
 387         if skipfirst > normsubpath.epsilon:
 388             normsubpathitems = normsubpath.segments([0, params[0]])[0]
 389             normsubpathitems.append(path.normcurve_pt(*(points[0][2:6] + points[1][0:4])))
 390         else:
 391             normsubpathitems = [path.normcurve_pt(*(points[0][2:6] + points[1][0:4]))]
 392         for i in range(1, len(points)-1):
 393             normsubpathitems.append(path.normcurve_pt(*(points[i][2:6] + points[i+1][0:4])))
 394         if skiplast > normsubpath.epsilon:
 395             for nsp in normsubpath.segments([params[-1], len(normsubpath)]):
 396                 normsubpathitems.extend(nsp.normsubpathitems)
 397
 398         # That's it
 399         return path.normsubpath(normsubpathitems, epsilon=normsubpath.epsilon)
 400 # >>>
 401
 402 cycloid.clear = attr.clearclass(cycloid)
 403
 404 class smoothed(deformer): # <<<
 405
 406     """Bends corners in a path.
 407
 408     This decorator replaces corners in a path with bezier curves. There are two cases:
 409     - If the corner lies between two lines, _two_ bezier curves will be used
 410       that are highly optimized to look good (their curvature is to be zero at the ends
 411       and has to have zero derivative in the middle).
 412       Additionally, it can controlled by the softness-parameter.
 413     - If the corner lies between curves then _one_ bezier is used that is (except in some
 414       special cases) uniquely determined by the tangents and curvatures at its end-points.
 415       In some cases it is necessary to use only the absolute value of the curvature to avoid a
 416       cusp-shaped connection of the new bezier to the old path. In this case the use of
 417       "obeycurv=0" allows the sign of the curvature to switch.
 418     - The radius argument gives the arclength-distance of the corner to the points where the
 419       old path is cut and the beziers are inserted.
 420     - Path elements that are too short (shorter than the radius) are skipped
 421     """
 422
 423     def __init__(self, radius, softness=1, obeycurv=0, relskipthres=0.01):
 424         self.radius = radius
 425         self.softness = softness
 426         self.obeycurv = obeycurv
 427         self.relskipthres = relskipthres
 428
 429     def __call__(self, radius=None, softness=None, obeycurv=None, relskipthres=None):
 430         if radius is None:
 431             radius = self.radius
 432         if softness is None:
 433             softness = self.softness
 434         if obeycurv is None:
 435             obeycurv = self.obeycurv
 436         if relskipthres is None:
 437             relskipthres = self.relskipthres
 438         return smoothed(radius=radius, softness=softness, obeycurv=obeycurv, relskipthres=relskipthres)
 439
 440     def deform(self, basepath):
 441         return path.normpath([self.deformsubpath(normsubpath)
 442                               for normsubpath in basepath.normpath().normsubpaths])
 443
 444     def deformsubpath(self, normsubpath):
 445         radius_pt = unit.topt(self.radius)
 446         epsilon = normsubpath.epsilon
 447
 448         # remove too short normsubpath items (shorter than self.relskipthres*radius_pt or epsilon)
 449         pertinentepsilon = max(epsilon, self.relskipthres*radius_pt)
 450         pertinentnormsubpath = path.normsubpath(normsubpath.normsubpathitems,
 451                                                 epsilon=pertinentepsilon)
 452         pertinentnormsubpath.flushskippedline()
 453         pertinentnormsubpathitems = pertinentnormsubpath.normsubpathitems
 454
 455         # calculate the splitting parameters for the pertinentnormsubpathitems
 456         arclens_pt = []
 457         params = []
 458         for pertinentnormsubpathitem in pertinentnormsubpathitems:
 459             arclen_pt = pertinentnormsubpathitem.arclen_pt(epsilon)
 460             arclens_pt.append(arclen_pt)
 461             l1_pt = min(radius_pt, 0.5*arclen_pt)
 462             l2_pt = max(0.5*arclen_pt, arclen_pt - radius_pt)
 463             params.append(pertinentnormsubpathitem.arclentoparam_pt([l1_pt, l2_pt], epsilon))
 464
 465         # handle the first and last pertinentnormsubpathitems for a non-closed normsubpath
 466         if not normsubpath.closed:
 467             l1_pt = 0
 468             l2_pt = max(0, arclens_pt[0] - radius_pt)
 469             params[0] = pertinentnormsubpathitem.arclentoparam_pt([l1_pt, l2_pt], epsilon)
 470             l1_pt = min(radius_pt, arclens_pt[-1])
 471             l2_pt = arclens_pt[-1]
 472             params[-1] = pertinentnormsubpathitem.arclentoparam_pt([l1_pt, l2_pt], epsilon)
 473
 474         newnormsubpath = path.normsubpath(epsilon=normsubpath.epsilon)
 475         for i in range(len(pertinentnormsubpathitems)):
 476             this = i
 477             next = (i+1) % len(pertinentnormsubpathitems)
 478             thisparams = params[this]
 479             nextparams = params[next]
 480             thisnormsubpathitem = pertinentnormsubpathitems[this]
 481             nextnormsubpathitem = pertinentnormsubpathitems[next]
 482             thisarclen_pt = arclens_pt[this]
 483             nextarclen_pt = arclens_pt[next]
 484
 485             # insert the middle segment
 486             newnormsubpath.append(thisnormsubpathitem.segments(thisparams)[0])
 487
 488             # insert replacement curves for the corners
 489             if next or normsubpath.closed:
 490
 491                 t1 = thisnormsubpathitem.rotation([thisparams[1]])[0].apply_pt(1, 0)
 492                 t2 = nextnormsubpathitem.rotation([nextparams[0]])[0].apply_pt(1, 0)
 493
 494                 if (isinstance(thisnormsubpathitem, path.normline_pt) and
 495                     isinstance(nextnormsubpathitem, path.normline_pt)):
 496
 497                     # case of two lines -> replace by two curves
 498                     d1, g1, f1, e, f2, g2, d2 = curvescontrols_from_endlines_pt(
 499                         thisnormsubpathitem.atend_pt(), t1, t2,
 500                         thisarclen_pt*(1-thisparams[1]), nextarclen_pt*(nextparams[0]), softness=self.softness)
 501
 502                     p1 = thisnormsubpathitem.at_pt([thisparams[1]])[0]
 503                     p2 = nextnormsubpathitem.at_pt([nextparams[0]])[0]
 504
 505                     newnormsubpath.append(path.normcurve_pt(*(d1 + g1 + f1 + e)))
 506                     newnormsubpath.append(path.normcurve_pt(*(e + f2 + g2 + d2)))
 507
 508                 else:
 509
 510                     # generic case -> replace by a single curve with prescribed tangents and curvatures
 511                     p1 = thisnormsubpathitem.at_pt([thisparams[1]])[0]
 512                     p2 = nextnormsubpathitem.at_pt([nextparams[0]])[0]
 513
 514                     # XXX supply curvature_pt methods in path module or transfere algorithms to work with curveradii
 515                     def curvature(normsubpathitem, param):
 516                         r = normsubpathitem.curveradius_pt([param])[0]
 517                         if r is None:
 518                             return 0
 519                         return 1.0 / r
 520
 521                     c1 = curvature(thisnormsubpathitem, thisparams[1])
 522                     c2 = curvature(nextnormsubpathitem, nextparams[0])
 523
 524                     if not self.obeycurv:
 525                         # do not obey the sign of the curvature but
 526                         # make the sign such that the curve smoothly passes to the next point
 527                         # this results in a discontinuous curvature
 528                         # (but the absolute value is still continuous)
 529                         s1 = sign1(t1[0] * (p2[1]-p1[1]) - t1[1] * (p2[0]-p1[0]))
 530                         s2 = sign1(t2[0] * (p2[1]-p1[1]) - t2[1] * (p2[0]-p1[0]))
 531                         c1 = s1 * abs(c1)
 532                         c2 = s2 * abs(c2)
 533
 534                     # get the length of the control "arms"
 535                     a, d = controldists_from_endpoints_pt(p1, p2, t1, t2, c1, c2)
 536
 537                     # avoid overshooting at the corners:
 538                     # this changes not only the sign of the curvature
 539                     # but also the magnitude
 540                     if not self.obeycurv:
 541                         t, s = intersection(p1, p2, t1, t2)
 542                         if t is None or t < 0:
 543                             a = None
 544                         else:
 545                             a = min(a, t)
 546
 547                         if s is None or s > 0:
 548                             d = None
 549                         else:
 550                             d = min(d, -s)
 551
 552                     # if there is no useful result:
 553                     # take arbitrary smoothing curve that does not obey
 554                     # the curvature constraints
 555                     if a is None or d is None:
 556                         dist = math.hypot(p1[0] - p2[0], p1[1] - p2[1])
 557                         a = dist / (3.0 * math.hypot(*t1))
 558                         d = dist / (3.0 * math.hypot(*t2))
 559
 560                     # calculate the two missing control points
 561                     q1 = p1[0] + a * t1[0], p1[1] + a * t1[1]
 562                     q2 = p2[0] - d * t2[0], p2[1] - d * t2[1]
 563
 564                     newnormsubpath.append(path.normcurve_pt(*(p1 + q1 + q2 + p2)))
 565
 566         if normsubpath.closed:
 567             newnormsubpath.close()
 568         return newnormsubpath
 569
 570 # >>>
 571
 572 smoothed.clear = attr.clearclass(smoothed)
 573
 574 class parallel(deformer): # <<<
 575
 576     """creates a parallel path with constant distance to the original path
 577
 578     A positive 'distance' results in a curve left of the original one -- and a
 579     negative 'distance' in a curve at the right. Left/Right are understood in
 580     terms of the parameterization of the original curve.
 581     At corners, either a circular arc is drawn around the corner, or, if the
 582     curve is on the other side, the parallel curve also exhibits a corner.
 583
 584     For each path element a parallel curve/line is constructed. For curves, the
 585     accuracy can be adjusted with the parameter 'relerr', thus, relerr*distance
 586     is the maximum allowable error somewhere in the middle of the curve (at
 587     parameter value 0.5).
 588     'relerr' only applies for the 'expensive' mode where the parallel curve for
 589     a single curve items may be composed of several (many) curve items.
 590     """
 591
 592     # TODO:
 593     # - check for greatest curvature and introduce extra corners
 594     #   if a normcurve is too heavily curved
 595     # - do relerr-checks at better points than just at parameter 0.5
 596
 597     def __init__(self, distance, relerr=0.05, expensive=1):
 598         self.distance = distance
 599         self.relerr = relerr
 600         self.expensive = expensive
 601
 602     def __call__(self, distance=None, relerr=None, expensive=None):
 603         # returns a copy of the deformer with different parameters
 604         if distance is None:
 605             d = self.distance
 606         if relerr is None:
 607             r = self.relerr
 608         if expensive is None:
 609             e = self.expensive
 610
 611         return parallel(distance=d, relerr=r, expensive=e)
 612
 613     def deform(self, basepath):
 614         resultnormsubpaths = [self.deformsubpath(nsp) for nsp in basepath.normpath().normsubpaths]
 615         return path.normpath(resultnormsubpaths)
 616
 617     def deformsubpath(self, orig_nspath):
 618
 619         distance = unit.topt(self.distance)
 620         relerr = self.relerr
 621         expensive = self.expensive
 622         epsilon = orig_nspath.epsilon
 623
 624         new_nspath = path.normsubpath(epsilon=epsilon)
 625
 626         # 1. Store endpoints, tangents and curvatures for each element
 627         points, tangents, curvatures = [], [], []
 628         for npitem in orig_nspath:
 629
 630             ps,ts,cs = [],[],[]
 631             trafos = npitem.trafo([0,1])
 632             for t in trafos:
 633                 p = t.apply_pt(0,0)
 634                 t = t.apply_pt(1,0)
 635                 ps.append(p)
 636                 ts.append((t[0]-p[0], t[1]-p[1]))
 637
 638             rs = npitem.curveradius_pt([0,1])
 639             cs = []
 640             for r in rs:
 641                 if r is None:
 642                     cs.append(0)
 643                 else:
 644                     cs.append(1.0 / r)
 645
 646             points.append(ps)
 647             tangents.append(ts)
 648             curvatures.append(cs)
 649
 650         closeparallel = (tangents[-1][1][0]*tangents[0][0][1] - tangents[-1][1][1]*tangents[0][0][0])
 651
 652         # 2. append the parallel path for each element:
 653         for cur in range(len(orig_nspath)):
 654
 655             if cur == 0:
 656                 old = cur
 657                 # OldEnd = points[old][0]
 658                 OldEndTang = tangents[old][0]
 659             else:
 660                 old = cur - 1
 661                 # OldEnd = points[old][1]
 662                 OldEndTang = tangents[old][1]
 663
 664             CurBeg, CurEnd = points[cur]
 665             CurBegTang, CurEndTang = tangents[cur]
 666             CurBegCurv, CurEndCurv = curvatures[cur]
 667
 668             npitem = orig_nspath[cur]
 669
 670             # get the control points for the shifted pathelement
 671             if isinstance(npitem, path.normline_pt):
 672                 # get the points explicitly from the normal vector
 673                 A = CurBeg[0] - distance * CurBegTang[1], CurBeg[1] + distance * CurBegTang[0]
 674                 D = CurEnd[0] - distance * CurEndTang[1], CurEnd[1] + distance * CurEndTang[0]
 675                 new_npitems = [path.normline_pt(A[0], A[1], D[0], D[1])]
 676             elif isinstance(npitem, path.normcurve_pt):
 677                 # call a function to return a list of controlpoints
 678                 cpoints_list = parallel_curvespoints_pt(npitem, distance, expensive, relerr, epsilon)
 679                 new_npitems = []
 680                 for cpoints in cpoints_list:
 681                     A,B,C,D = cpoints
 682                     new_npitems.append(path.normcurve_pt(A[0],A[1], B[0],B[1], C[0],C[1], D[0],D[1]))
 683                 # we will need the starting point of the new normpath items
 684                 A = cpoints_list[0][0]
 685
 686
 687             # append the next piece of the path:
 688             # it might contain of an extra arc or must be intersected before appending
 689             parallel = (OldEndTang[0]*CurBegTang[1] - OldEndTang[1]*CurBegTang[0])
 690             if parallel*distance < -epsilon:
 691
 692                 # append an arc around the corner
 693                 # from the preceding piece to the current piece
 694                 # we can never get here for the first npitem! (because cur==old)
 695                 endpoint = new_nspath.atend_pt()
 696                 center = CurBeg
 697                 angle1 = math.atan2(endpoint[1] - center[1], endpoint[0] - center[0]) * 180.0 / math.pi
 698                 angle2 = math.atan2(A[1] - center[1], A[0] - center[0]) * 180.0 / math.pi
 699                 if parallel > 0:
 700                     arc_npath = path.path(path.arc_pt(
 701                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 702                 else:
 703                     arc_npath = path.path(path.arcn_pt(
 704                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 705
 706                 for new_npitem in arc_npath[0]:
 707                     new_nspath.append(new_npitem)
 708
 709
 710             elif parallel*distance > epsilon:
 711                 # intersect the extra piece of the path with the rest of the new path
 712                 # and throw away the void parts
 713                 #
 714                 # build a subpath for intersection
 715                 extra_nspath = path.normsubpath(normsubpathitems=new_npitems, epsilon=epsilon)
 716
 717                 intsparams = extra_nspath.intersect(new_nspath)
 718                 # [[a,b,c], [a,b,c]]
 719                 if intsparams:
 720                     # take the first intersection point:
 721                     extra_param, new_param = intsparams[0][0], intsparams[1][0]
 722                     new_nspath = new_nspath.segments([0, new_param])[0]
 723                     extra_nspath = extra_nspath.segments([extra_param, len(extra_nspath)])[0]
 724                     new_npitems = extra_nspath.normsubpathitems
 725                     # in case the intersection was not sufficiently exact:
 726                     # CAREFUL! because we newly created all the new_npitems and
 727                     # the items in extra_nspath, we may in-place change the starting point
 728                     new_npitems[0] = new_npitems[0].modifiedbegin_pt(*new_nspath.atend_pt())
 729                 else:
 730                     raise # how did we get here?
 731
 732
 733             # at the (possible) closing corner we may have to intersect another time
 734             # or add another corner:
 735             # the intersection must be done before appending the parallel piece
 736             if orig_nspath.closed and cur == len(orig_nspath) - 1:
 737                 if closeparallel * distance > epsilon:
 738                     intsparams = extra_nspath.intersect(new_nspath)
 739                     # [[a,b,c], [a,b,c]]
 740                     if intsparams:
 741                         # take the last intersection point:
 742                         extra_param, new_param = intsparams[0][-1], intsparams[1][-1]
 743                         new_nspath = new_nspath.segments([new_param, len(new_nspath)])[0]
 744                         extra_nspath = extra_nspath.segments([0, extra_param])[0]
 745                         new_npitems = extra_nspath.normsubpathitems
 746                         # in case the intersection was not sufficiently exact:
 747                         # CAREFUL! because we newly created all the new_npitems and
 748                         # the items in extra_nspath, we may in-place change the end point
 749                         new_npitems[-1] = new_npitems[-1].modifiedend_pt(*new_nspath.atbegin_pt())
 750                     else:
 751                         raise # how did we get here?
 752
 753                     pass
 754
 755
 756             # append the parallel piece
 757             for new_npitem in new_npitems:
 758                 new_nspath.append(new_npitem)
 759
 760
 761         # the curve around the closing corner must be added at last:
 762         if orig_nspath.closed:
 763             if closeparallel * distance < -epsilon:
 764                 endpoint = new_nspath.atend_pt()
 765                 center = orig_nspath.atend_pt()
 766                 A = new_nspath.atbegin_pt()
 767                 angle1 = math.atan2(endpoint[1] - center[1], endpoint[0] - center[0]) * 180.0 / math.pi
 768                 angle2 = math.atan2(A[1] - center[1], A[0] - center[0]) * 180.0 / math.pi
 769                 if parallel > 0:
 770                     arc_npath = path.path(path.arc_pt(
 771                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 772                 else:
 773                     arc_npath = path.path(path.arcn_pt(
 774                       center[0], center[1], abs(distance), angle1, angle2)).normpath()
 775
 776                 for new_npitem in arc_npath[0]:
 777                     new_nspath.append(new_npitem)
 778
 779         # 3. extra handling of closed paths
 780         if orig_nspath.closed:
 781             new_nspath.close()
 782
 783         return new_nspath
 784 # >>>
 785
 786 parallel.clear = attr.clearclass(parallel)
 787
 788 # vim:foldmethod=marker:foldmarker=<<<,>>>