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