deformer.py

   1 # -*- encoding: utf-8 -*-
   2 #
   3 #
   4 # Copyright (C) 2003-2006 Michael Schindler <m-schindler@users.sourceforge.net>
   5 # Copyright (C) 2003-2005 André Wobst <wobsta@users.sourceforge.net>
   6 #
   7 # This file is part of PyX (http://pyx.sourceforge.net/).
   8 #
   9 # PyX is free software; you can redistribute it and/or modify
  10 # it under the terms of the GNU General Public License as published by
  11 # the Free Software Foundation; either version 2 of the License, or
  12 # (at your option) any later version.
  13 #
  14 # PyX is distributed in the hope that it will be useful,
  15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17 # GNU General Public License for more details.
  18 #
  19 # You should have received a copy of the GNU General Public License
  20 # along with PyX; if not, write to the Free Software
  21 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  22
  23 import math, warnings
  24 from . import attr, baseclasses, mathutils, path, normpath, unit, color
  25
  26 # specific exception for an invalid parameterization point
  27 # used in parallel
  28 class InvalidParamException(Exception):
  29
  30     def __init__(self, param):
  31         self.normsubpathitemparam = param
  32
  33 # None has a meaning in linesmoothed
  34 class _marker: pass
  35
  36 def curvescontrols_from_endlines_pt(B, tangent1, tangent2, r1, r2, softness): # <<<
  37     # calculates the parameters for two bezier curves connecting two lines (curvature=0)
  38     # starting at B - r1*tangent1
  39     # ending at   B + r2*tangent2
  40     #
  41     # Takes the corner B
  42     # and two tangent vectors heading to and from B
  43     # and two radii r1 and r2:
  44     # All arguments must be in Points
  45     # Returns the seven control points of the two bezier curves:
  46     #  - start d1
  47     #  - control points g1 and f1
  48     #  - midpoint e
  49     #  - control points f2 and g2
  50     #  - endpoint d2
  51
  52     # make direction vectors d1: from B to A
  53     #                        d2: from B to C
  54     d1 = -tangent1[0] / math.hypot(*tangent1), -tangent1[1] / math.hypot(*tangent1)
  55     d2 =  tangent2[0] / math.hypot(*tangent2),  tangent2[1] / math.hypot(*tangent2)
  56
  57     # 0.3192 has turned out to be the maximum softness available
  58     # for straight lines ;-)
  59     f = 0.3192 * softness
  60     g = (15.0 * f + math.sqrt(-15.0*f*f + 24.0*f))/12.0
  61
  62     # make the control points of the two bezier curves
  63     f1 = B[0] + f * r1 * d1[0], B[1] + f * r1 * d1[1]
  64     f2 = B[0] + f * r2 * d2[0], B[1] + f * r2 * d2[1]
  65     g1 = B[0] + g * r1 * d1[0], B[1] + g * r1 * d1[1]
  66     g2 = B[0] + g * r2 * d2[0], B[1] + g * r2 * d2[1]
  67     d1 = B[0] +     r1 * d1[0], B[1] +     r1 * d1[1]
  68     d2 = B[0] +     r2 * d2[0], B[1] +     r2 * d2[1]
  69     e  = 0.5 * (f1[0] + f2[0]), 0.5 * (f1[1] + f2[1])
  70
  71     return (d1, g1, f1, e, f2, g2, d2)
  72 # >>>
  73
  74 def controldists_from_endgeometry_pt(A, B, tangA, tangB, curvA, curvB, allownegative=0): # <<<
  75
  76     """For a curve with given tangents and curvatures at the endpoints this gives the distances between the controlpoints
  77
  78     This helper routine returns a list of two distances between the endpoints and the
  79     corresponding control points of a (cubic) bezier curve that has
  80     prescribed tangents tangentA, tangentB and curvatures curvA, curvB at the
  81     end points.
  82
  83     Note: The returned distances are not always positive.
  84           But only positive values are geometrically correct, so please check!
  85           The outcome is sorted so that the first entry is expected to be the
  86           most reasonable one
  87     """
  88     debug = 0
  89
  90     def test_divisions(T, D, E, AB, curvA, curvB, debug):# <<<
  91
  92         def is_zero(x):
  93             try:
  94                 1.0 / x
  95             except ZeroDivisionError:
  96                 return 1
  97             return 0
  98
  99         T_is_zero = is_zero(T)
 100         curvA_is_zero = is_zero(curvA)
 101         curvB_is_zero = is_zero(curvB)
 102
 103         if T_is_zero:
 104             if curvA_is_zero:
 105                 assert abs(D) < 1.0e-10
 106                 a = AB / 3.0
 107                 if curvB_is_zero:
 108                     assert abs(E) < 1.0e-10
 109                     b = AB / 3.0
 110                 else:
 111                     b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
 112             else:
 113                 a = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA)
 114                 if curvB_is_zero:
 115                     assert abs(E) < 1.0e-10
 116                     b = AB / 3.0
 117                 else:
 118                     b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
 119         else:
 120             if curvA_is_zero:
 121                 b = D / T
 122                 a = (E - 1.5*curvB*b*abs(b)) / T
 123             elif curvB_is_zero:
 124                 a = E / T
 125                 b = (D - 1.5*curvA*a*abs(a)) / T
 126             else:
 127                 return []
 128
 129         if debug:
 130             print("fallback with exact zero value")
 131         return [(a, b)]
 132     # >>>
 133     def fallback_smallT(T, D, E, AB, curvA, curvB, threshold, debug):# <<<
 134         a = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA)
 135         b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
 136         q1 = min(abs(1.5*a*a*curvA), abs(D))
 137         q2 = min(abs(1.5*b*b*curvB), abs(E))
 138         if (a >= 0 and b >= 0 and
 139             abs(b*T) < threshold * q1 and abs(1.5*a*abs(a)*curvA - D) < threshold * q1 and
 140             abs(a*T) < threshold * q2 and abs(1.5*b*abs(b)*curvB - E) < threshold * q2):
 141             if debug:
 142                 print("fallback with T approx 0")
 143             return [(a, b)]
 144         return []
 145     # >>>
 146     def fallback_smallcurv(T, D, E, AB, curvA, curvB, threshold, debug):# <<<
 147         result = []
 148
 149         # is curvB approx zero?
 150         a = E / T
 151         b = (D - 1.5*curvA*a*abs(a)) / T
 152         if (a >= 0 and b >= 0 and
 153             abs(1.5*b*b*curvB) < threshold * min(abs(a*T), abs(E)) and
 154             abs(a*T - E) < threshold * min(abs(a*T), abs(E))):
 155             if debug:
 156                 print("fallback with curvB approx 0")
 157             result.append((a, b))
 158
 159         # is curvA approx zero?
 160         b = D / T
 161         a = (E - 1.5*curvB*b*abs(b)) / T
 162         if (a >= 0 and b >= 0 and
 163             abs(1.5*a*a*curvA) < threshold * min(abs(b*T), abs(D)) and
 164             abs(b*T - D) < threshold * min(abs(b*T), abs(D))):
 165             if debug:
 166                 print("fallback with curvA approx 0")
 167             result.append((a, b))
 168
 169         return result
 170     # >>>
 171     def findnearest(x, ys): # <<<
 172         I = 0
 173         Y = ys[I]
 174         mindist = abs(x - Y)
 175
 176         # find the value in ys which is nearest to x
 177         for i, y in enumerate(ys[1:]):
 178             dist = abs(x - y)
 179             if dist < mindist:
 180                 I, Y, mindist = i, y, dist
 181
 182         return I, Y
 183     # >>>
 184
 185     # some shortcuts
 186     T = tangA[0] * tangB[1] - tangA[1] * tangB[0]
 187     D = tangA[0] * (B[1]-A[1]) - tangA[1] * (B[0]-A[0])
 188     E = tangB[0] * (A[1]-B[1]) - tangB[1] * (A[0]-B[0])
 189     AB = math.hypot(A[0] - B[0], A[1] - B[1])
 190
 191     # try if one of the prefactors is exactly zero
 192     testsols = test_divisions(T, D, E, AB, curvA, curvB, debug)
 193     if testsols:
 194         return testsols
 195
 196     # The general case:
 197     # we try to find all the zeros of the decoupled 4th order problem
 198     # for the combined problem:
 199     # The control points of a cubic Bezier curve are given by a, b:
 200     #     A, A + a*tangA, B - b*tangB, B
 201     # for the derivation see /design/beziers.tex
 202     #     0 = 1.5 a |a| curvA + b * T - D
 203     #     0 = 1.5 b |b| curvB + a * T - E
 204     # because of the absolute values we get several possibilities for the signs
 205     # in the equation. We test all signs, also the invalid ones!
 206     if allownegative:
 207         signs = [(+1, +1), (-1, +1), (+1, -1), (-1, -1)]
 208     else:
 209         signs = [(+1, +1)]
 210
 211     candidates_a = []
 212     candidates_b = []
 213     for sign_a, sign_b in signs:
 214         coeffs_a = (sign_b*3.375*curvA*curvA*curvB, 0.0, -sign_b*sign_a*4.5*curvA*curvB*D, T**3, sign_b*1.5*curvB*D*D - T*T*E)
 215         coeffs_b = (sign_a*3.375*curvA*curvB*curvB, 0.0, -sign_a*sign_b*4.5*curvA*curvB*E, T**3, sign_a*1.5*curvA*E*E - T*T*D)
 216         candidates_a += [root for root in mathutils.realpolyroots(*coeffs_a) if sign_a*root >= 0]
 217         candidates_b += [root for root in mathutils.realpolyroots(*coeffs_b) if sign_b*root >= 0]
 218     solutions = []
 219     if candidates_a and candidates_b:
 220         for a in candidates_a:
 221             i, b = findnearest((D - 1.5*curvA*a*abs(a))/T, candidates_b)
 222             solutions.append((a, b))
 223
 224     # try if there is an approximate solution
 225     for thr in [1.0e-2, 1.0e-1]:
 226         if not solutions:
 227             solutions = fallback_smallT(T, D, E, AB, curvA, curvB, thr, debug)
 228         if not solutions:
 229             solutions = fallback_smallcurv(T, D, E, AB, curvA, curvB, thr, debug)
 230
 231     # sort the solutions: the more reasonable values at the beginning
 232     def mycmp(x,y): # <<<
 233         # first the pairs that are purely positive, then all the pairs with some negative signs
 234         # inside the two sets: sort by magnitude
 235         sx = (x[0] > 0 and x[1] > 0)
 236         sy = (y[0] > 0 and y[1] > 0)
 237
 238         # experimental stuff:
 239         # what criterion should be used for sorting ?
 240         #
 241         #errx = abs(1.5*curvA*x[0]*abs(x[0]) + x[1]*T - D) + abs(1.5*curvB*x[1]*abs(x[1]) + x[0]*T - E)
 242         #erry = abs(1.5*curvA*y[0]*abs(y[0]) + y[1]*T - D) + abs(1.5*curvB*y[1]*abs(y[1]) + y[0]*T - E)
 243         # # For each equation, a value like
 244         # #   abs(1.5*curvA*y[0]*abs(y[0]) + y[1]*T - D) / abs(curvA*(D - y[1]*T))
 245         # # indicates how good the solution is. In order to avoid the division,
 246         # # we here multiply with all four denominators:
 247         # errx = max(abs( (1.5*curvA*y[0]*abs(y[0]) + y[1]*T - D) * (curvB*(E - y[0]*T))*(curvA*(D - x[1]*T))*(curvB*(E - x[0]*T)) ),
 248         #            abs( (1.5*curvB*y[1]*abs(y[1]) + y[0]*T - E) * (curvA*(D - y[1]*T))*(curvA*(D - x[1]*T))*(curvB*(E - x[0]*T)) ))
 249         # errx = max(abs( (1.5*curvA*x[0]*abs(x[0]) + x[1]*T - D) * (curvA*(D - y[1]*T))*(curvB*(E - y[0]*T))*(curvB*(E - x[0]*T)) ),
 250         #            abs( (1.5*curvB*x[1]*abs(x[1]) + x[0]*T - E) * (curvA*(D - y[1]*T))*(curvB*(E - y[0]*T))*(curvA*(D - x[1]*T)) ))
 251         #errx = (abs(curvA*x[0]) - 1.0)**2 + (abs(curvB*x[1]) - 1.0)**2
 252         #erry = (abs(curvA*y[0]) - 1.0)**2 + (abs(curvB*y[1]) - 1.0)**2
 253
 254         errx = x[0]**2 + x[1]**2
 255         erry = y[0]**2 + y[1]**2
 256
 257         if sx == 1 and sy == 1:
 258             # try to use longer solutions if there are any crossings in the control-arms
 259             # the following combination yielded fewest sorting errors in test_bezier.py
 260             t, s = intersection(A, B, tangA, tangB)
 261             t, s = abs(t), abs(s)
 262             if (t > 0 and t < x[0] and s > 0 and s < x[1]):
 263                 if (t > 0 and t < y[0] and s > 0 and s < y[1]):
 264                     # use the shorter one
 265                     return cmp(errx, erry)
 266                 else:
 267                     # use the longer one
 268                     return -1
 269             else:
 270                 if (t > 0 and t < y[0] and s > 0 and s < y[1]):
 271                     # use the longer one
 272                     return 1
 273                 else:
 274                     # use the shorter one
 275                     return cmp(errx, erry)
 276             #return cmp(x[0]**2 + x[1]**2, y[0]**2 + y[1]**2)
 277         else:
 278             return cmp(sy, sx)
 279     # >>>
 280     solutions.sort(mycmp)
 281
 282     return solutions
 283 # >>>
 284
 285 def normcurve_from_endgeometry_pt(A, B, tangA, tangB, curvA, curvB): # <<<
 286     a, b = controldists_from_endgeometry_pt(A, B, tangA, tangB, curvA, curvB)[0]
 287     return normpath.normcurve_pt(A[0], A[1],
 288         A[0] + a * tangA[0], A[1] + a * tangA[1],
 289         B[0] - b * tangB[0], B[1] - b * tangB[1], B[0], B[1])
 290     # >>>
 291
 292 def intersection(A, D, tangA, tangD): # <<<
 293
 294     """returns the intersection parameters of two evens
 295
 296     they are defined by:
 297       x(t) = A + t * tangA
 298       x(s) = D + s * tangD
 299     """
 300     det = -tangA[0] * tangD[1] + tangA[1] * tangD[0]
 301     try:
 302         1.0 / det
 303     except ArithmeticError:
 304         return None, None
 305
 306     DA = D[0] - A[0], D[1] - A[1]
 307
 308     t = (-tangD[1]*DA[0] + tangD[0]*DA[1]) / det
 309     s = (-tangA[1]*DA[0] + tangA[0]*DA[1]) / det
 310
 311     return t, s
 312 # >>>
 313
 314 class cycloid(baseclasses.deformer): # <<<
 315     """Wraps a cycloid around a path.
 316
 317     The outcome looks like a spring with the originalpath as the axis.
 318     radius: radius of the cycloid
 319     halfloops:  number of halfloops
 320     skipfirst/skiplast: undeformed end lines of the original path
 321     curvesperhloop:
 322     sign: start left (1) or right (-1) with the first halfloop
 323     turnangle: angle of perspective on a (3D) spring
 324                turnangle=0 will produce a sinus-like cycloid,
 325                turnangle=90 will procude a row of connected circles
 326
 327     """
 328
 329     def __init__(self, radius=0.5*unit.t_cm, halfloops=10,
 330     skipfirst=1*unit.t_cm, skiplast=1*unit.t_cm, curvesperhloop=3, sign=1, turnangle=45):
 331         self.skipfirst = skipfirst
 332         self.skiplast = skiplast
 333         self.radius = radius
 334         self.halfloops = halfloops
 335         self.curvesperhloop = curvesperhloop
 336         self.sign = sign
 337         self.turnangle = turnangle
 338
 339     def __call__(self, radius=None, halfloops=None,
 340     skipfirst=None, skiplast=None, curvesperhloop=None, sign=None, turnangle=None):
 341         if radius is None:
 342             radius = self.radius
 343         if halfloops is None:
 344             halfloops = self.halfloops
 345         if skipfirst is None:
 346             skipfirst = self.skipfirst
 347         if skiplast is None:
 348             skiplast = self.skiplast
 349         if curvesperhloop is None:
 350             curvesperhloop = self.curvesperhloop
 351         if sign is None:
 352             sign = self.sign
 353         if turnangle is None:
 354             turnangle = self.turnangle
 355
 356         return cycloid(radius=radius, halfloops=halfloops, skipfirst=skipfirst, skiplast=skiplast,
 357                        curvesperhloop=curvesperhloop, sign=sign, turnangle=turnangle)
 358
 359     def deform(self, basepath):
 360         resultnormsubpaths = [self.deformsubpath(nsp) for nsp in basepath.normpath().normsubpaths]
 361         return normpath.normpath(resultnormsubpaths)
 362
 363     def deformsubpath(self, normsubpath):
 364
 365         skipfirst = abs(unit.topt(self.skipfirst))
 366         skiplast = abs(unit.topt(self.skiplast))
 367         radius = abs(unit.topt(self.radius))
 368         turnangle = math.radians(self.turnangle)
 369         sign = mathutils.sign(self.sign)
 370
 371         cosTurn = math.cos(turnangle)
 372         sinTurn = math.sin(turnangle)
 373
 374         # make list of the lengths and parameters at points on normsubpath
 375         # where we will add cycloid-points
 376         totlength = normsubpath.arclen_pt()
 377         if totlength <= skipfirst + skiplast + 2*radius*sinTurn:
 378             warnings.warn("normsubpath is too short for deformation with cycloid -- skipping...")
 379             return normsubpath
 380
 381         # parameterization is in rotation-angle around the basepath
 382         # differences in length, angle ... between two basepoints
 383         # and between basepoints and controlpoints
 384         Dphi = math.pi / self.curvesperhloop
 385         phis = [i * Dphi for i in range(self.halfloops * self.curvesperhloop + 1)]
 386         DzDphi = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * math.pi * cosTurn)
 387         # Dz = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * self.curvesperhloop * cosTurn)
 388         # zs = [i * Dz for i in range(self.halfloops * self.curvesperhloop + 1)]
 389         # from path._arctobcurve:
 390         # optimal relative distance along tangent for second and third control point
 391         L = 4 * radius * (1 - math.cos(Dphi/2)) / (3 * math.sin(Dphi/2))
 392
 393         # Now the transformation of z into the turned coordinate system
 394         Zs = [ skipfirst + radius*sinTurn # here the coordinate z starts
 395              - sinTurn*radius*math.cos(phi) + cosTurn*DzDphi*phi # the transformed z-coordinate
 396              for phi in phis]
 397         params = normsubpath._arclentoparam_pt(Zs)[0]
 398
 399         # get the positions of the splitpoints in the cycloid
 400         points = []
 401         for phi, param in zip(phis, params):
 402             # the cycloid is a circle that is stretched along the normsubpath
 403             # here are the points of that circle
 404             basetrafo = normsubpath.trafo([param])[0]
 405
 406             # The point on the cycloid, in the basepath's local coordinate system
 407             baseZ, baseY = 0, radius*math.sin(phi)
 408
 409             # The tangent there, also in local coords
 410             tangentX = -cosTurn*radius*math.sin(phi) + sinTurn*DzDphi
 411             tangentY = radius*math.cos(phi)
 412             tangentZ = sinTurn*radius*math.sin(phi) + DzDphi*cosTurn
 413             norm = math.sqrt(tangentX*tangentX + tangentY*tangentY + tangentZ*tangentZ)
 414             tangentY, tangentZ = tangentY/norm, tangentZ/norm
 415
 416             # Respect the curvature of the basepath for the cycloid's curvature
 417             # XXX this is only a heuristic, not a "true" expression for
 418             #     the curvature in curved coordinate systems
 419             pathradius = normsubpath.curveradius_pt([param])[0]
 420             if pathradius is not normpath.invalid:
 421                 factor = (pathradius - baseY) / pathradius
 422                 factor = abs(factor)
 423             else:
 424                 factor = 1
 425             l = L * factor
 426
 427             # The control points prior and after the point on the cycloid
 428             preeZ, preeY = baseZ - l * tangentZ, baseY - l * tangentY
 429             postZ, postY = baseZ + l * tangentZ, baseY + l * tangentY
 430
 431             # Now put everything at the proper place
 432             points.append(basetrafo.apply_pt(preeZ, sign * preeY) +
 433                           basetrafo.apply_pt(baseZ, sign * baseY) +
 434                           basetrafo.apply_pt(postZ, sign * postY))
 435
 436         if len(points) <= 1:
 437             warnings.warn("normsubpath is too short for deformation with cycloid -- skipping...")
 438             return normsubpath
 439
 440         # Build the path from the pointlist
 441         # containing (control x 2,  base x 2, control x 2)
 442         if skipfirst > normsubpath.epsilon:
 443             normsubpathitems = normsubpath.segments([0, params[0]])[0]
 444             normsubpathitems.append(normpath.normcurve_pt(*(points[0][2:6] + points[1][0:4])))
 445         else:
 446             normsubpathitems = [normpath.normcurve_pt(*(points[0][2:6] + points[1][0:4]))]
 447         for i in range(1, len(points)-1):
 448             normsubpathitems.append(normpath.normcurve_pt(*(points[i][2:6] + points[i+1][0:4])))
 449         if skiplast > normsubpath.epsilon:
 450             for nsp in normsubpath.segments([params[-1], len(normsubpath)]):
 451                 normsubpathitems.extend(nsp.normsubpathitems)
 452
 453         # That's it
 454         return normpath.normsubpath(normsubpathitems, epsilon=normsubpath.epsilon)
 455 # >>>
 456
 457 cycloid.clear = attr.clearclass(cycloid)
 458
 459 class cornersmoothed(baseclasses.deformer): # <<<
 460
 461     """Bends corners in a normpath.
 462
 463     This decorator replaces corners in a normpath with bezier curves. There are two cases:
 464     - If the corner lies between two lines, _two_ bezier curves will be used
 465       that are highly optimized to look good (their curvature is to be zero at the ends
 466       and has to have zero derivative in the middle).
 467       Additionally, it can controlled by the softness-parameter.
 468     - If the corner lies between curves then _one_ bezier is used that is (except in some
 469       special cases) uniquely determined by the tangents and curvatures at its end-points.
 470       In some cases it is necessary to use only the absolute value of the curvature to avoid a
 471       cusp-shaped connection of the new bezier to the old path. In this case the use of
 472       "obeycurv=0" allows the sign of the curvature to switch.
 473     - The radius argument gives the arclength-distance of the corner to the points where the
 474       old path is cut and the beziers are inserted.
 475     - Path elements that are too short (shorter than the radius) are skipped
 476     """
 477
 478     def __init__(self, radius, softness=1, obeycurv=0, relskipthres=0.01):
 479         self.radius = radius
 480         self.softness = softness
 481         self.obeycurv = obeycurv
 482         self.relskipthres = relskipthres
 483
 484     def __call__(self, radius=None, softness=None, obeycurv=None, relskipthres=None):
 485         if radius is None:
 486             radius = self.radius
 487         if softness is None:
 488             softness = self.softness
 489         if obeycurv is None:
 490             obeycurv = self.obeycurv
 491         if relskipthres is None:
 492             relskipthres = self.relskipthres
 493         return cornersmoothed(radius=radius, softness=softness, obeycurv=obeycurv, relskipthres=relskipthres)
 494
 495     def deform(self, basepath):
 496         return normpath.normpath([self.deformsubpath(normsubpath)
 497                               for normsubpath in basepath.normpath().normsubpaths])
 498
 499     def deformsubpath(self, normsubpath):
 500         radius_pt = unit.topt(self.radius)
 501         epsilon = normsubpath.epsilon
 502
 503         # remove too short normsubpath items (shorter than self.relskipthres*radius_pt or epsilon)
 504         pertinentepsilon = max(epsilon, self.relskipthres*radius_pt)
 505         pertinentnormsubpath = normpath.normsubpath(normsubpath.normsubpathitems,
 506                                                 epsilon=pertinentepsilon)
 507         pertinentnormsubpath.flushskippedline()
 508         pertinentnormsubpathitems = pertinentnormsubpath.normsubpathitems
 509
 510         # calculate the splitting parameters for the pertinentnormsubpathitems
 511         arclens_pt = []
 512         params = []
 513         for pertinentnormsubpathitem in pertinentnormsubpathitems:
 514             arclen_pt = pertinentnormsubpathitem.arclen_pt(epsilon)
 515             arclens_pt.append(arclen_pt)
 516             l1_pt = min(radius_pt, 0.5*arclen_pt)
 517             l2_pt = max(0.5*arclen_pt, arclen_pt - radius_pt)
 518             params.append(pertinentnormsubpathitem.arclentoparam_pt([l1_pt, l2_pt], epsilon))
 519
 520         # handle the first and last pertinentnormsubpathitems for a non-closed normsubpath
 521         if not normsubpath.closed:
 522             l1_pt = 0
 523             l2_pt = max(0, arclens_pt[0] - radius_pt)
 524             params[0] = pertinentnormsubpathitems[0].arclentoparam_pt([l1_pt, l2_pt], epsilon)
 525             l1_pt = min(radius_pt, arclens_pt[-1])
 526             l2_pt = arclens_pt[-1]
 527             params[-1] = pertinentnormsubpathitems[-1].arclentoparam_pt([l1_pt, l2_pt], epsilon)
 528
 529         newnormsubpath = normpath.normsubpath(epsilon=normsubpath.epsilon)
 530         for i in range(len(pertinentnormsubpathitems)):
 531             this = i
 532             next = (i+1) % len(pertinentnormsubpathitems)
 533             thisparams = params[this]
 534             nextparams = params[next]
 535             thisnormsubpathitem = pertinentnormsubpathitems[this]
 536             nextnormsubpathitem = pertinentnormsubpathitems[next]
 537             thisarclen_pt = arclens_pt[this]
 538             nextarclen_pt = arclens_pt[next]
 539
 540             # insert the middle segment
 541             newnormsubpath.append(thisnormsubpathitem.segments(thisparams)[0])
 542
 543             # insert replacement curves for the corners
 544             if next or normsubpath.closed:
 545
 546                 t1 = thisnormsubpathitem.rotation([thisparams[1]])[0].apply_pt(1, 0)
 547                 t2 = nextnormsubpathitem.rotation([nextparams[0]])[0].apply_pt(1, 0)
 548                 # TODO: normpath.invalid
 549
 550                 if (isinstance(thisnormsubpathitem, normpath.normline_pt) and
 551                     isinstance(nextnormsubpathitem, normpath.normline_pt)):
 552
 553                     # case of two lines -> replace by two curves
 554                     d1, g1, f1, e, f2, g2, d2 = curvescontrols_from_endlines_pt(
 555                         thisnormsubpathitem.atend_pt(), t1, t2,
 556                         thisarclen_pt*(1-thisparams[1]), nextarclen_pt*(nextparams[0]), softness=self.softness)
 557
 558                     p1 = thisnormsubpathitem.at_pt([thisparams[1]])[0]
 559                     p2 = nextnormsubpathitem.at_pt([nextparams[0]])[0]
 560
 561                     newnormsubpath.append(normpath.normcurve_pt(*(d1 + g1 + f1 + e)))
 562                     newnormsubpath.append(normpath.normcurve_pt(*(e + f2 + g2 + d2)))
 563
 564                 else:
 565
 566                     # generic case -> replace by a single curve with prescribed tangents and curvatures
 567                     p1 = thisnormsubpathitem.at_pt([thisparams[1]])[0]
 568                     p2 = nextnormsubpathitem.at_pt([nextparams[0]])[0]
 569                     c1 = thisnormsubpathitem.curvature_pt([thisparams[1]])[0]
 570                     c2 = nextnormsubpathitem.curvature_pt([nextparams[0]])[0]
 571                     # TODO: normpath.invalid
 572
 573                     # TODO: more intelligent fallbacks:
 574                     #   circle -> line
 575                     #   circle -> circle
 576
 577                     if not self.obeycurv:
 578                         # do not obey the sign of the curvature but
 579                         # make the sign such that the curve smoothly passes to the next point
 580                         # this results in a discontinuous curvature
 581                         # (but the absolute value is still continuous)
 582                         s1 = +mathutils.sign(t1[0] * (p2[1]-p1[1]) - t1[1] * (p2[0]-p1[0]))
 583                         s2 = -mathutils.sign(t2[0] * (p2[1]-p1[1]) - t2[1] * (p2[0]-p1[0]))
 584                         c1 = s1 * abs(c1)
 585                         c2 = s2 * abs(c2)
 586
 587                     # get the length of the control "arms"
 588                     controldists = controldists_from_endgeometry_pt(p1, p2, t1, t2, c1, c2)
 589
 590                     if controldists and (controldists[0][0] >= 0 and controldists[0][1] >= 0):
 591                         # use the first entry in the controldists
 592                         # this should be the "smallest" pair
 593                         a, d = controldists[0]
 594                         # avoid curves with invalid parameterization
 595                         a = max(a, epsilon)
 596                         d = max(d, epsilon)
 597
 598                         # avoid overshooting at the corners:
 599                         # this changes not only the sign of the curvature
 600                         # but also the magnitude
 601                         if not self.obeycurv:
 602                             t, s = intersection(p1, p2, t1, t2)
 603                             if (t is not None and s is not None and
 604                                 t > 0 and s < 0):
 605                                 a = min(a, abs(t))
 606                                 d = min(d, abs(s))
 607
 608                     else:
 609                         # use a fallback
 610                         t, s = intersection(p1, p2, t1, t2)
 611                         if t is not None and s is not None:
 612                             a = 0.65 * abs(t)
 613                             d = 0.65 * abs(s)
 614                         else:
 615                             # if there is no useful result:
 616                             # take an arbitrary smoothing curve that does not obey
 617                             # the curvature constraints
 618                             dist = math.hypot(p1[0] - p2[0], p1[1] - p2[1])
 619                             a = dist / (3.0 * math.hypot(*t1))
 620                             d = dist / (3.0 * math.hypot(*t2))
 621
 622                     # calculate the two missing control points
 623                     q1 = p1[0] + a * t1[0], p1[1] + a * t1[1]
 624                     q2 = p2[0] - d * t2[0], p2[1] - d * t2[1]
 625
 626                     newnormsubpath.append(normpath.normcurve_pt(*(p1 + q1 + q2 + p2)))
 627
 628         if normsubpath.closed:
 629             newnormsubpath.close()
 630         return newnormsubpath
 631
 632 # >>>
 633
 634 cornersmoothed.clear = attr.clearclass(cornersmoothed)
 635 smoothed = cornersmoothed
 636 smoothed.clear = attr.clearclass(smoothed)
 637
 638 class parallel(baseclasses.deformer): # <<<
 639
 640     """creates a parallel normpath with constant distance to the original normpath
 641
 642     A positive 'distance' results in a curve left of the original one -- and a
 643     negative 'distance' in a curve at the right. Left/Right are understood in
 644     terms of the parameterization of the original curve. For each path element
 645     a parallel curve/line is constructed. At corners, either a circular arc is
 646     drawn around the corner, or, if possible, the parallel curve is cut in
 647     order to also exhibit a corner.
 648
 649     distance:            the distance of the parallel normpath
 650     relerr:              distance*relerr is the maximal allowed error in the parallel distance
 651     sharpoutercorners:   make the outer corners not round but sharp.
 652                          The inner corners (corners after inflection points) will stay round
 653     dointersection:      boolean for doing the intersection step (default: 1).
 654                          Set this value to 0 if you want the whole parallel path
 655     checkdistanceparams: a list of parameter values in the interval (0,1) where the
 656                          parallel distance is checked on each normpathitem
 657     lookforcurvatures:   number of points per normpathitem where is looked for
 658                          a critical value of the curvature
 659     """
 660
 661     # TODO:
 662     # * do testing for curv=0, T=0, D=0, E=0 cases
 663     # * do testing for several random curves
 664     #   -- does the recursive deformnicecurve converge?
 665
 666
 667     def __init__(self, distance, relerr=0.05, sharpoutercorners=0, dointersection=1,
 668                        checkdistanceparams=[0.5], lookforcurvatures=11, debug=None):
 669         self.distance = distance
 670         self.relerr = relerr
 671         self.sharpoutercorners = sharpoutercorners
 672         self.checkdistanceparams = checkdistanceparams
 673         self.lookforcurvatures = lookforcurvatures
 674         self.dointersection = dointersection
 675         self.debug = debug
 676
 677     def __call__(self, distance=None, relerr=None, sharpoutercorners=None, dointersection=None,
 678                        checkdistanceparams=None, lookforcurvatures=None, debug=None):
 679         # returns a copy of the deformer with different parameters
 680         if distance is None:
 681             distance = self.distance
 682         if relerr is None:
 683             relerr = self.relerr
 684         if sharpoutercorners is None:
 685             sharpoutercorners = self.sharpoutercorners
 686         if dointersection is None:
 687             dointersection = self.dointersection
 688         if checkdistanceparams is None:
 689             checkdistanceparams = self.checkdistanceparams
 690         if lookforcurvatures is None:
 691             lookforcurvatures = self.lookforcurvatures
 692         if debug is None:
 693             debug = self.debug
 694
 695         return parallel(distance=distance, relerr=relerr,
 696                         sharpoutercorners=sharpoutercorners,
 697                         dointersection=dointersection,
 698                         checkdistanceparams=checkdistanceparams,
 699                         lookforcurvatures=lookforcurvatures,
 700                         debug=debug)
 701
 702     def deform(self, basepath):
 703         self.dist_pt = unit.topt(self.distance)
 704         resultnormsubpaths = []
 705         for nsp in basepath.normpath().normsubpaths:
 706             parallel_normpath = self.deformsubpath(nsp)
 707             resultnormsubpaths += parallel_normpath.normsubpaths
 708         result = normpath.normpath(resultnormsubpaths)
 709         return result
 710
 711     def deformsubpath(self, orig_nsp): # <<<
 712
 713         """returns a list of normsubpaths building the parallel curve"""
 714
 715         dist = self.dist_pt
 716         epsilon = orig_nsp.epsilon
 717
 718         # avoid too small dists: we would run into instabilities
 719         if abs(dist) < abs(epsilon):
 720             return normpath.normpath([orig_nsp])
 721
 722         result = normpath.normpath()
 723
 724         # iterate over the normsubpath in the following manner:
 725         # * for each item first append the additional arc / intersect
 726         #   and then add the next parallel piece
 727         # * for the first item only add the parallel piece
 728         #   (because this is done for next_orig_nspitem, we need to start with next=0)
 729         for i in range(len(orig_nsp.normsubpathitems)):
 730             prev = i-1
 731             next = i
 732             prev_orig_nspitem = orig_nsp.normsubpathitems[prev]
 733             next_orig_nspitem = orig_nsp.normsubpathitems[next]
 734
 735             stepsize = 0.01
 736             prev_param, prev_rotation = self.valid_near_rotation(prev_orig_nspitem, 1, 0, stepsize, 0.5*epsilon)
 737             next_param, next_rotation = self.valid_near_rotation(next_orig_nspitem, 0, 1, stepsize, 0.5*epsilon)
 738             # TODO: eventually shorten next_orig_nspitem
 739
 740             prev_tangent = prev_rotation.apply_pt(1, 0)
 741             next_tangent = next_rotation.apply_pt(1, 0)
 742
 743             # get the next parallel piece for the normpath
 744             try:
 745                 next_parallel_normpath = self.deformsubpathitem(next_orig_nspitem, epsilon)
 746             except InvalidParamException as e:
 747                 invalid_nspitem_param = e.normsubpathitemparam
 748                 # split the nspitem apart and continue with pieces that do not contain
 749                 # the invalid point anymore. At the end, simply take one piece, otherwise two.
 750                 stepsize = 0.01
 751                 if self.length_pt(next_orig_nspitem, invalid_nspitem_param, 0) > epsilon:
 752                     if self.length_pt(next_orig_nspitem, invalid_nspitem_param, 1) > epsilon:
 753                         p1, foo = self.valid_near_rotation(next_orig_nspitem, invalid_nspitem_param, 0, stepsize, 0.5*epsilon)
 754                         p2, foo = self.valid_near_rotation(next_orig_nspitem, invalid_nspitem_param, 1, stepsize, 0.5*epsilon)
 755                         segments = next_orig_nspitem.segments([0, p1, p2, 1])
 756                         segments = segments[0], segments[2].modifiedbegin_pt(*(segments[0].atend_pt()))
 757                     else:
 758                         p1, foo = self.valid_near_rotation(next_orig_nspitem, invalid_nspitem_param, 0, stepsize, 0.5*epsilon)
 759                         segments = next_orig_nspitem.segments([0, p1])
 760                 else:
 761                     p2, foo = self.valid_near_rotation(next_orig_nspitem, invalid_nspitem_param, 1, stepsize, 0.5*epsilon)
 762                     segments = next_orig_nspitem.segments([p2, 1])
 763
 764                 next_parallel_normpath = self.deformsubpath(normpath.normsubpath(segments, epsilon=epsilon))
 765
 766             if not (next_parallel_normpath.normsubpaths and next_parallel_normpath[0].normsubpathitems):
 767                 continue
 768
 769             # this starts the whole normpath
 770             if not result.normsubpaths:
 771                 result = next_parallel_normpath
 772                 continue
 773
 774             # sinus of the angle between the tangents
 775             # sinangle > 0 for a left-turning nexttangent
 776             # sinangle < 0 for a right-turning nexttangent
 777             sinangle = prev_tangent[0]*next_tangent[1] - prev_tangent[1]*next_tangent[0]
 778             cosangle = prev_tangent[0]*next_tangent[0] + prev_tangent[1]*next_tangent[1]
 779             if cosangle < 0 or abs(dist*math.asin(sinangle)) >= epsilon:
 780                 if self.sharpoutercorners and dist*sinangle < 0:
 781                     A1, A2 = result.atend_pt(), next_parallel_normpath.atbegin_pt()
 782                     t1, t2 = intersection(A1, A2, prev_tangent, next_tangent)
 783                     B = A1[0] + t1 * prev_tangent[0], A1[1] + t1 * prev_tangent[1]
 784                     arc_normpath = normpath.normpath([normpath.normsubpath([
 785                         normpath.normline_pt(A1[0], A1[1], B[0], B[1]),
 786                         normpath.normline_pt(B[0], B[1], A2[0], A2[1])
 787                         ])])
 788                 else:
 789                     # We must append an arc around the corner
 790                     arccenter = next_orig_nspitem.atbegin_pt()
 791                     arcbeg = result.atend_pt()
 792                     arcend = next_parallel_normpath.atbegin_pt()
 793                     angle1 = math.atan2(arcbeg[1] - arccenter[1], arcbeg[0] - arccenter[0])
 794                     angle2 = math.atan2(arcend[1] - arccenter[1], arcend[0] - arccenter[0])
 795
 796                     # depending on the direction we have to use arc or arcn
 797                     if dist > 0:
 798                         arcclass = path.arcn_pt
 799                     else:
 800                         arcclass = path.arc_pt
 801                     arc_normpath = path.path(arcclass(
 802                       arccenter[0], arccenter[1], abs(dist),
 803                       math.degrees(angle1), math.degrees(angle2))).normpath(epsilon=epsilon)
 804
 805                 # append the arc to the parallel path
 806                 result.join(arc_normpath)
 807                 # append the next parallel piece to the path
 808                 result.join(next_parallel_normpath)
 809             else:
 810                 # The path is quite straight between prev and next item:
 811                 # normpath.normpath.join adds a straight line if necessary
 812                 result.join(next_parallel_normpath)
 813
 814
 815         # end here if nothing has been found so far
 816         if not (result.normsubpaths and result[-1].normsubpathitems):
 817             return result
 818
 819         # the curve around the closing corner may still be missing
 820         if orig_nsp.closed:
 821             # TODO: normpath.invalid
 822             stepsize = 0.01
 823             prev_param, prev_rotation = self.valid_near_rotation(result[-1][-1], 1, 0, stepsize, 0.5*epsilon)
 824             next_param, next_rotation = self.valid_near_rotation(result[0][0], 0, 1, stepsize, 0.5*epsilon)
 825             # TODO: eventually shorten next_orig_nspitem
 826
 827             prev_tangent = prev_rotation.apply_pt(1, 0)
 828             next_tangent = next_rotation.apply_pt(1, 0)
 829             sinangle = prev_tangent[0]*next_tangent[1] - prev_tangent[1]*next_tangent[0]
 830             cosangle = prev_tangent[0]*next_tangent[0] + prev_tangent[1]*next_tangent[1]
 831
 832             if cosangle < 0 or abs(dist*math.asin(sinangle)) >= epsilon:
 833                 # We must append an arc around the corner
 834                 # TODO: avoid the code dublication
 835                 if self.sharpoutercorners and dist*sinangle < 0:
 836                     A1, A2 = result.atend_pt(), result.atbegin_pt()
 837                     t1, t2 = intersection(A1, A2, prev_tangent, next_tangent)
 838                     B = A1[0] + t1 * prev_tangent[0], A1[1] + t1 * prev_tangent[1]
 839                     arc_normpath = normpath.normpath([normpath.normsubpath([
 840                         normpath.normline_pt(A1[0], A1[1], B[0], B[1]),
 841                         normpath.normline_pt(B[0], B[1], A2[0], A2[1])
 842                         ])])
 843                 else:
 844                     arccenter = orig_nsp.atend_pt()
 845                     arcbeg = result.atend_pt()
 846                     arcend = result.atbegin_pt()
 847                     angle1 = math.atan2(arcbeg[1] - arccenter[1], arcbeg[0] - arccenter[0])
 848                     angle2 = math.atan2(arcend[1] - arccenter[1], arcend[0] - arccenter[0])
 849
 850                     # depending on the direction we have to use arc or arcn
 851                     if dist > 0:
 852                         arcclass = path.arcn_pt
 853                     else:
 854                         arcclass = path.arc_pt
 855                     arc_normpath = path.path(arcclass(
 856                         arccenter[0], arccenter[1], abs(dist),
 857                         math.degrees(angle1), math.degrees(angle2))).normpath(epsilon=epsilon)
 858
 859                 # append the arc to the parallel path
 860                 if (result.normsubpaths and result[-1].normsubpathitems and
 861                     arc_normpath.normsubpaths and arc_normpath[-1].normsubpathitems):
 862                     result.join(arc_normpath)
 863
 864             if len(result) == 1:
 865                 result[0].close()
 866             else:
 867                 # if the parallel normpath is split into several subpaths anyway,
 868                 # then use the natural beginning and ending
 869                 # closing is not possible anymore
 870                 for nspitem in result[0]:
 871                     result[-1].append(nspitem)
 872                 result.normsubpaths = result.normsubpaths[1:]
 873
 874         if self.dointersection:
 875             result = self.rebuild_intersected_normpath(result, normpath.normpath([orig_nsp]), epsilon)
 876
 877         return result
 878         # >>>
 879     def deformsubpathitem(self, nspitem, epsilon): # <<<
 880
 881         """Returns a parallel normpath for a single normsubpathitem
 882
 883         Analyzes the curvature of a normsubpathitem and returns a normpath with
 884         the appropriate number of normsubpaths. This must be a normpath because
 885         a normcurve can be strongly curved, such that the parallel path must
 886         contain a hole"""
 887
 888         dist = self.dist_pt
 889
 890         # for a simple line we return immediately
 891         if isinstance(nspitem, normpath.normline_pt):
 892             normal = nspitem.rotation([0])[0].apply_pt(0, 1)
 893             start = nspitem.atbegin_pt()
 894             end = nspitem.atend_pt()
 895             return path.line_pt(
 896                 start[0] + dist * normal[0], start[1] + dist * normal[1],
 897                 end[0] + dist * normal[0], end[1] + dist * normal[1]).normpath(epsilon=epsilon)
 898
 899         # for a curve we have to check if the curvatures
 900         # cross the singular value 1/dist
 901         crossings = self.distcrossingparameters(nspitem, epsilon)
 902
 903         # depending on the number of crossings we must consider
 904         # three different cases:
 905         if crossings:
 906             # The curvature crosses the borderline 1/dist
 907             # the parallel curve contains points with infinite curvature!
 908             result = normpath.normpath()
 909
 910             # we need the endpoints of the nspitem
 911             if self.length_pt(nspitem, crossings[0], 0) > epsilon:
 912                 crossings.insert(0, 0)
 913             if self.length_pt(nspitem, crossings[-1], 1) > epsilon:
 914                 crossings.append(1)
 915
 916             for i in range(len(crossings) - 1):
 917                 middleparam = 0.5*(crossings[i] + crossings[i+1])
 918                 middlecurv = nspitem.curvature_pt([middleparam])[0]
 919                 if middlecurv is normpath.invalid:
 920                     raise InvalidParamException(middleparam)
 921                 # the radius is good if
 922                 #  - middlecurv and dist have opposite signs or
 923                 #  - middlecurv is "smaller" than 1/dist
 924                 if middlecurv*dist < 0 or abs(dist*middlecurv) < 1:
 925                     parallel_nsp = self.deformnicecurve(nspitem.segments(crossings[i:i+2])[0], epsilon)
 926                     # never append empty normsubpaths
 927                     if parallel_nsp.normsubpathitems:
 928                         result.append(parallel_nsp)
 929
 930             return result
 931
 932         else:
 933             # the curvature is either bigger or smaller than 1/dist
 934             middlecurv = nspitem.curvature_pt([0.5])[0]
 935             if dist*middlecurv < 0 or abs(dist*middlecurv) < 1:
 936                 # The curve is everywhere less curved than 1/dist
 937                 # We can proceed finding the parallel curve for the whole piece
 938                 parallel_nsp = self.deformnicecurve(nspitem, epsilon)
 939                 # never append empty normsubpaths
 940                 if parallel_nsp.normsubpathitems:
 941                     return normpath.normpath([parallel_nsp])
 942                 else:
 943                     return normpath.normpath()
 944             else:
 945                 # the curve is everywhere stronger curved than 1/dist
 946                 # There is nothing to be returned.
 947                 return normpath.normpath()
 948
 949         # >>>
 950     def deformnicecurve(self, normcurve, epsilon, startparam=0.0, endparam=1.0): # <<<
 951
 952         """Returns a parallel normsubpath for the normcurve.
 953
 954         This routine assumes that the normcurve is everywhere
 955         'less' curved than 1/dist and contains no point with an
 956         invalid parameterization
 957         """
 958         dist = self.dist_pt
 959         T_threshold = 1.0e-5
 960
 961         # normalized tangent directions
 962         tangA, tangD = normcurve.rotation([startparam, endparam])
 963         # if we find an unexpected normpath.invalid we have to
 964         # parallelise this normcurve on the level of split normsubpaths
 965         if tangA is normpath.invalid:
 966             raise InvalidParamException(startparam)
 967         if tangD is normpath.invalid:
 968             raise InvalidParamException(endparam)
 969         tangA = tangA.apply_pt(1, 0)
 970         tangD = tangD.apply_pt(1, 0)
 971
 972         # the new starting points
 973         orig_A, orig_D = normcurve.at_pt([startparam, endparam])
 974         A = orig_A[0] - dist * tangA[1], orig_A[1] + dist * tangA[0]
 975         D = orig_D[0] - dist * tangD[1], orig_D[1] + dist * tangD[0]
 976
 977         # we need to end this _before_ we will run into epsilon-problems
 978         # when creating curves we do not want to calculate the length of
 979         # or even split it for recursive calls
 980         if (math.hypot(A[0] - D[0], A[1] - D[1]) < epsilon and
 981             math.hypot(tangA[0] - tangD[0], tangA[1] - tangD[1]) < T_threshold):
 982             return normpath.normsubpath([normpath.normline_pt(A[0], A[1], D[0], D[1])])
 983
 984         result = normpath.normsubpath(epsilon=epsilon)
 985         # is there enough space on the normals before they intersect?
 986         a, d = intersection(orig_A, orig_D, (-tangA[1], tangA[0]), (-tangD[1], tangD[0]))
 987         # a,d are the lengths to the intersection points:
 988         # for a (and equally for b) we can proceed in one of the following cases:
 989         #   a is None (means parallel normals)
 990         #   a and dist have opposite signs (and the same for b)
 991         #   a has the same sign but is bigger
 992         if ( (a is None or a*dist < 0 or abs(a) > abs(dist) + epsilon) or
 993              (d is None or d*dist < 0 or abs(d) > abs(dist) + epsilon) ):
 994             # the original path is long enough to draw a parallel piece
 995             # this is the generic case. Get the parallel curves
 996             orig_curvA, orig_curvD = normcurve.curvature_pt([startparam, endparam])
 997             # normpath.invalid may not appear here because we have asked
 998             # for this already at the tangents
 999             assert orig_curvA is not normpath.invalid
1000             assert orig_curvD is not normpath.invalid
1001             curvA = orig_curvA / (1.0 - dist*orig_curvA)
1002             curvD = orig_curvD / (1.0 - dist*orig_curvD)
1003
1004             # first try to approximate the normcurve with a single item
1005             controldistpairs = controldists_from_endgeometry_pt(A, D, tangA, tangD, curvA, curvD)
1006
1007             if controldistpairs:
1008                 # TODO: is it good enough to get the first entry here?
1009                 #       from testing: this fails if there are loops in the original curve
1010                 a, d = controldistpairs[0]
1011                 if a >= 0 and d >= 0:
1012                     if a < epsilon and d < epsilon:
1013                         result = normpath.normsubpath([normpath.normline_pt(A[0], A[1], D[0], D[1])], epsilon=epsilon)
1014                     else:
1015                         # we avoid curves with invalid parameterization
1016                         a = max(a, epsilon)
1017                         d = max(d, epsilon)
1018                         result = normpath.normsubpath([normpath.normcurve_pt(
1019                             A[0], A[1],
1020                             A[0] + a * tangA[0], A[1] + a * tangA[1],
1021                             D[0] - d * tangD[0], D[1] - d * tangD[1],
1022                             D[0], D[1])], epsilon=epsilon)
1023
1024             # then try with two items, recursive call
1025             if ((not result.normsubpathitems) or
1026                 (self.checkdistanceparams and result.normsubpathitems
1027                  and not self.distchecked(normcurve, result, epsilon, startparam, endparam))):
1028                 # TODO: does this ever converge?
1029                 # TODO: what if this hits epsilon?
1030                 firstnsp = self.deformnicecurve(normcurve, epsilon, startparam, 0.5*(startparam+endparam))
1031                 secondnsp = self.deformnicecurve(normcurve, epsilon, 0.5*(startparam+endparam), endparam)
1032                 if not (firstnsp.normsubpathitems and secondnsp.normsubpathitems):
1033                     result = normpath.normsubpath(
1034                         [normpath.normline_pt(A[0], A[1], D[0], D[1])], epsilon=epsilon)
1035                 else:
1036                     # we will get problems if the curves are too short:
1037                     result = firstnsp.joined(secondnsp)
1038
1039         return result
1040         # >>>
1041
1042     def distchecked(self, orig_normcurve, parallel_normsubpath, epsilon, tstart, tend): # <<<
1043
1044         """Checks the distances between orig_normcurve and parallel_normsubpath
1045
1046         The checking is done at parameters self.checkdistanceparams of orig_normcurve."""
1047
1048         dist = self.dist_pt
1049         # do not look closer than epsilon:
1050         dist_relerr = mathutils.sign(dist) * max(abs(self.relerr*dist), epsilon)
1051
1052         checkdistanceparams = [tstart + (tend-tstart)*t for t in self.checkdistanceparams]
1053
1054         for param, P, rotation in zip(checkdistanceparams,
1055                                       orig_normcurve.at_pt(checkdistanceparams),
1056                                       orig_normcurve.rotation(checkdistanceparams)):
1057             # check if the distance is really the wanted distance
1058             # measure the distance in the "middle" of the original curve
1059             if rotation is normpath.invalid:
1060                 raise InvalidParamException(param)
1061
1062             normal = rotation.apply_pt(0, 1)
1063
1064             # create a short cutline for intersection only:
1065             cutline = normpath.normsubpath([normpath.normline_pt (
1066               P[0] + (dist - 2*dist_relerr) * normal[0],
1067               P[1] + (dist - 2*dist_relerr) * normal[1],
1068               P[0] + (dist + 2*dist_relerr) * normal[0],
1069               P[1] + (dist + 2*dist_relerr) * normal[1])], epsilon=epsilon)
1070
1071             cutparams = parallel_normsubpath.intersect(cutline)
1072             distances = [math.hypot(P[0] - cutpoint[0], P[1] - cutpoint[1])
1073                          for cutpoint in cutline.at_pt(cutparams[1])]
1074
1075             if (not distances) or (abs(min(distances) - abs(dist)) > abs(dist_relerr)):
1076                 return 0
1077
1078         return 1
1079     # >>>
1080     def distcrossingparameters(self, normcurve, epsilon, tstart=0, tend=1): # <<<
1081
1082         """Returns a list of parameters where the curvature is 1/distance"""
1083
1084         dist = self.dist_pt
1085
1086         # we _need_ to do this with the curvature, not with the radius
1087         # because the curvature is continuous at the straight line and the radius is not:
1088         # when passing from one slightly curved curve to the other with opposite curvature sign,
1089         # via the straight line, then the curvature changes its sign at curv=0, while the
1090         # radius changes its sign at +/-infinity
1091         # this causes instabilities for nearly straight curves
1092
1093         # include tstart and tend
1094         params = [tstart + i * (tend - tstart) * 1.0 / (self.lookforcurvatures - 1)
1095                   for i in range(self.lookforcurvatures)]
1096         curvs = normcurve.curvature_pt(params)
1097
1098         # break everything at invalid curvatures
1099         for param, curv in zip(params, curvs):
1100             if curv is normpath.invalid:
1101                 raise InvalidParamException(param)
1102
1103         parampairs = list(zip(params[:-1], params[1:]))
1104         curvpairs = list(zip(curvs[:-1], curvs[1:]))
1105
1106         crossingparams = []
1107         for parampair, curvpair in zip(parampairs, curvpairs):
1108             begparam, endparam = parampair
1109             begcurv, endcurv = curvpair
1110             if (endcurv*dist - 1)*(begcurv*dist - 1) < 0:
1111                 # the curvature crosses the value 1/dist
1112                 # get the parmeter value by linear interpolation:
1113                 middleparam = (
1114                   (begparam * abs(begcurv*dist - 1) + endparam * abs(endcurv*dist - 1)) /
1115                   (abs(begcurv*dist - 1) + abs(endcurv*dist - 1)))
1116                 middleradius = normcurve.curveradius_pt([middleparam])[0]
1117
1118                 if middleradius is normpath.invalid:
1119                     raise InvalidParamException(middleparam)
1120
1121                 if abs(middleradius - dist) < epsilon:
1122                     # get the parmeter value by linear interpolation:
1123                     crossingparams.append(middleparam)
1124                 else:
1125                     # call recursively:
1126                     cps = self.distcrossingparameters(normcurve, epsilon, tstart=begparam, tend=endparam)
1127                     crossingparams += cps
1128
1129         return crossingparams
1130         # >>>
1131     def valid_near_rotation(self, nspitem, param, otherparam, stepsize, epsilon): # <<<
1132         p = param
1133         rot = nspitem.rotation([p])[0]
1134         # run towards otherparam searching for a valid rotation
1135         while rot is normpath.invalid:
1136             p = (1-stepsize)*p + stepsize*otherparam
1137             rot = nspitem.rotation([p])[0]
1138         # walk back to param until near enough
1139         # but do not go further if an invalid point is hit
1140         end, new = nspitem.at_pt([param, p])
1141         far = math.hypot(end[0]-new[0], end[1]-new[1])
1142         pnew = p
1143         while far > epsilon:
1144             pnew = (1-stepsize)*pnew + stepsize*param
1145             end, new = nspitem.at_pt([param, pnew])
1146             far = math.hypot(end[0]-new[0], end[1]-new[1])
1147             if nspitem.rotation([pnew])[0] is normpath.invalid:
1148                 break
1149             else:
1150                 p = pnew
1151         return p, nspitem.rotation([p])[0]
1152     # >>>
1153     def length_pt(self, path, param1, param2): # <<<
1154         point1, point2 = path.at_pt([param1, param2])
1155         return math.hypot(point1[0] - point2[0], point1[1] - point2[1])
1156     # >>>
1157
1158     def normpath_selfintersections(self, np, epsilon): # <<<
1159
1160         """return all self-intersection points of normpath np.
1161
1162         This does not include the intersections of a single normcurve with itself,
1163         but all intersections of one normpathitem with a different one in the path"""
1164
1165         n = len(np)
1166         linearparams = []
1167         parampairs = []
1168         paramsriap = {}
1169         for nsp_i in range(n):
1170             for nsp_j in range(nsp_i, n):
1171                 for nspitem_i in range(len(np[nsp_i])):
1172                     if nsp_j == nsp_i:
1173                         nspitem_j_range = list(range(nspitem_i+1, len(np[nsp_j])))
1174                     else:
1175                         nspitem_j_range = list(range(len(np[nsp_j])))
1176                     for nspitem_j in nspitem_j_range:
1177                         intsparams = np[nsp_i][nspitem_i].intersect(np[nsp_j][nspitem_j], epsilon)
1178                         if intsparams:
1179                             for intsparam_i, intsparam_j in intsparams:
1180                                 if ( (abs(intsparam_i) < epsilon and abs(1-intsparam_j) < epsilon) or
1181                                      (abs(intsparam_j) < epsilon and abs(1-intsparam_i) < epsilon) ):
1182                                      continue
1183                                 npp_i = normpath.normpathparam(np, nsp_i, float(nspitem_i)+intsparam_i)
1184                                 npp_j = normpath.normpathparam(np, nsp_j, float(nspitem_j)+intsparam_j)
1185                                 linearparams.append(npp_i)
1186                                 linearparams.append(npp_j)
1187                                 paramsriap[id(npp_i)] = len(parampairs)
1188                                 paramsriap[id(npp_j)] = len(parampairs)
1189                                 parampairs.append((npp_i, npp_j))
1190         linearparams.sort()
1191         return linearparams, parampairs, paramsriap
1192
1193     # >>>
1194     def can_continue(self, par_np, param1, param2): # <<<
1195         dist = self.dist_pt
1196
1197         rot1, rot2 = par_np.rotation([param1, param2])
1198         if rot1 is normpath.invalid or rot2 is normpath.invalid:
1199             return 0
1200         curv1, curv2 = par_np.curvature_pt([param1, param2])
1201         tang2 = rot2.apply_pt(1, 0)
1202         norm1 = rot1.apply_pt(0, -1)
1203         norm1 = (dist*norm1[0], dist*norm1[1])
1204
1205         # the self-intersection is valid if the tangents
1206         # point into the correct direction or, for parallel tangents,
1207         # if the curvature is such that the on-going path does not
1208         # enter the region defined by dist
1209         mult12 = norm1[0]*tang2[0] + norm1[1]*tang2[1]
1210         eps = 1.0e-6
1211         if abs(mult12) > eps:
1212             return (mult12 < 0)
1213         else:
1214             # tang1 and tang2 are parallel
1215             if curv2 is normpath.invalid or curv1 is normpath.invalid:
1216                 return 0
1217             if dist > 0:
1218                 return (curv2 <= curv1)
1219             else:
1220                 return (curv2 >= curv1)
1221     # >>>
1222     def rebuild_intersected_normpath(self, par_np, orig_np, epsilon): # <<<
1223
1224         dist = self.dist_pt
1225
1226         # calculate the self-intersections of the par_np
1227         selfintparams, selfintpairs, selfintsriap = self.normpath_selfintersections(par_np, epsilon)
1228         # calculate the intersections of the par_np with the original path
1229         origintparams = par_np.intersect(orig_np)[0]
1230
1231         # visualize the intersection points: # <<<
1232         if self.debug is not None:
1233             for param1, param2 in selfintpairs:
1234                 point1, point2 = par_np.at([param1, param2])
1235                 self.debug.fill(path.circle(point1[0], point1[1], 0.05), [color.rgb.red])
1236                 self.debug.fill(path.circle(point2[0], point2[1], 0.03), [color.rgb.black])
1237             for param in origintparams:
1238                 point = par_np.at([param])[0]
1239                 self.debug.fill(path.circle(point[0], point[1], 0.05), [color.rgb.green])
1240         # >>>
1241
1242         result = normpath.normpath()
1243         if not selfintparams:
1244             if origintparams:
1245                 return result
1246             else:
1247                 return par_np
1248
1249         beginparams = []
1250         endparams = []
1251         for i in range(len(par_np)):
1252             beginparams.append(normpath.normpathparam(par_np, i, 0))
1253             endparams.append(normpath.normpathparam(par_np, i, len(par_np[i])))
1254
1255         allparams = selfintparams + origintparams + beginparams + endparams
1256         allparams.sort()
1257         allparamindices = {}
1258         for i, param in enumerate(allparams):
1259             allparamindices[id(param)] = i
1260
1261         done = {}
1262         for param in allparams:
1263             done[id(param)] = 0
1264
1265         def otherparam(p): # <<<
1266             pair = selfintpairs[selfintsriap[id(p)]]
1267             if (p is pair[0]):
1268                 return pair[1]
1269             else:
1270                 return pair[0]
1271         # >>>
1272         def trial_parampairs(startp): # <<<
1273             tried = {}
1274             for param in allparams:
1275                 tried[id(param)] = done[id(param)]
1276
1277             lastp = startp
1278             currentp = allparams[allparamindices[id(startp)] + 1]
1279             result = []
1280
1281             while 1:
1282                 if currentp is startp:
1283                     result.append((lastp, currentp))
1284                     return result
1285                 if currentp in selfintparams and otherparam(currentp) is startp:
1286                     result.append((lastp, currentp))
1287                     return result
1288                 if currentp in endparams:
1289                     result.append((lastp, currentp))
1290                     return result
1291                 if tried[id(currentp)]:
1292                     return []
1293                 if currentp in origintparams:
1294                     return []
1295                 # follow the crossings on valid startpairs until
1296                 # the normsubpath is closed or the end is reached
1297                 if (currentp in selfintparams and
1298                     self.can_continue(par_np, currentp, otherparam(currentp))):
1299                     # go to the next pair on the curve, seen from currentpair[1]
1300                     result.append((lastp, currentp))
1301                     lastp = otherparam(currentp)
1302                     tried[id(currentp)] = 1
1303                     tried[id(otherparam(currentp))] = 1
1304                     currentp = allparams[allparamindices[id(otherparam(currentp))] + 1]
1305                 else:
1306                     # go to the next pair on the curve, seen from currentpair[0]
1307                     tried[id(currentp)] = 1
1308                     tried[id(otherparam(currentp))] = 1
1309                     currentp = allparams[allparamindices[id(currentp)] + 1]
1310             assert 0
1311         # >>>
1312
1313         # first the paths that start at the beginning of a subnormpath:
1314         for startp in beginparams + selfintparams:
1315             if done[id(startp)]:
1316                 continue
1317
1318             parampairs = trial_parampairs(startp)
1319             if not parampairs:
1320                 continue
1321
1322             # collect all the pieces between parampairs
1323             add_nsp = normpath.normsubpath(epsilon=epsilon)
1324             for begin, end in parampairs:
1325                 # check that trial_parampairs works correctly
1326                 assert begin is not end
1327                 # we do not cross the border of a normsubpath here
1328                 assert begin.normsubpathindex is end.normsubpathindex
1329                 for item in par_np[begin.normsubpathindex].segments(
1330                     [begin.normsubpathparam, end.normsubpathparam])[0].normsubpathitems:
1331                     # TODO: this should be obsolete with an improved intersection algorithm
1332                     #       guaranteeing epsilon
1333                     if add_nsp.normsubpathitems:
1334                         item = item.modifiedbegin_pt(*(add_nsp.atend_pt()))
1335                     add_nsp.append(item)
1336
1337                 if begin in selfintparams:
1338                     done[id(begin)] = 1
1339                     #done[otherparam(begin)] = 1
1340                 if end in selfintparams:
1341                     done[id(end)] = 1
1342                     #done[otherparam(end)] = 1
1343
1344             # eventually close the path
1345             if add_nsp and (parampairs[0][0] is parampairs[-1][-1] or
1346                 (parampairs[0][0] in selfintparams and otherparam(parampairs[0][0]) is parampairs[-1][-1])):
1347                 add_nsp.normsubpathitems[-1] = add_nsp.normsubpathitems[-1].modifiedend_pt(*add_nsp.atbegin_pt())
1348                 add_nsp.close()
1349
1350             result.extend([add_nsp])
1351
1352         return result
1353
1354     # >>>
1355
1356 # >>>
1357
1358 parallel.clear = attr.clearclass(parallel)
1359
1360 class linesmoothed(baseclasses.deformer): # <<<
1361
1362     def __init__(self, tension=1, atleast=False, lcurl=1, rcurl=1):
1363         """Tension and atleast control the tension of the replacement curves.
1364         l/rcurl control the curlynesses at (possible) endpoints. If a curl is
1365         set to None, the angle is taken from the original path."""
1366         if atleast:
1367             self.tension = -abs(tension)
1368         else:
1369             self.tension = abs(tension)
1370         self.lcurl = lcurl
1371         self.rcurl = rcurl
1372
1373     def __call__(self, tension=_marker, atleast=_marker, lcurl=_marker, rcurl=_marker):
1374         if tension is _marker:
1375             tension = self.tension
1376         if atleast is _marker:
1377             atleast = (self.tension < 0)
1378         if lcurl is _marker:
1379             lcurl = self.lcurl
1380         if rcurl is _marker:
1381             rcurl = self.rcurl
1382         return linesmoothed(tension, atleast, lcurl, rcurl)
1383
1384     def deform(self, basepath):
1385         newnp = normpath.normpath()
1386         for nsp in basepath.normpath().normsubpaths:
1387             newnp += self.deformsubpath(nsp)
1388         return newnp
1389
1390     def deformsubpath(self, nsp):
1391         from metapost import path as mppath
1392         """Returns a path/normpath from the points in the given normsubpath"""
1393         # TODO: epsilon ?
1394         knots = []
1395
1396         # first point
1397         x_pt, y_pt = nsp.atbegin_pt()
1398         if nsp.closed:
1399             knots.append(mppath.smoothknot_pt(x_pt, y_pt))
1400         elif self.lcurl is None:
1401             rot = nsp.rotation([0])[0]
1402             dx, dy = rot.apply_pt(1, 0)
1403             angle = math.atan2(dy, dx)
1404             knots.append(mppath.beginknot_pt(x_pt, y_pt, angle=angle))
1405         else:
1406             knots.append(mppath.beginknot_pt(x_pt, y_pt, curl=self.lcurl))
1407
1408         # intermediate points:
1409         for npelem in nsp[:-1]:
1410             knots.append(mppath.tensioncurve(self.tension))
1411             knots.append(mppath.smoothknot_pt(*npelem.atend_pt()))
1412
1413         # last point
1414         knots.append(mppath.tensioncurve(self.tension))
1415         x_pt, y_pt = nsp.atend_pt()
1416         if nsp.closed:
1417             pass
1418         elif self.rcurl is None:
1419             rot = nsp.rotation([len(nsp)])[0]
1420             dx, dy = rot.apply_pt(1, 0)
1421             angle = math.atan2(dy, dx)
1422             knots.append(mppath.endknot_pt(x_pt, y_pt, angle=angle))
1423         else:
1424             knots.append(mppath.endknot_pt(x_pt, y_pt, curl=self.rcurl))
1425
1426         return mppath.path(knots)
1427 # >>>
1428
1429 linesmoothed.clear = attr.clearclass(linesmoothed)
1430
1431
1432 # vim:foldmethod=marker:foldmarker=<<<,>>>