pyx/deformer.py

   1 #!/usr/bin/env python
   2 # -*- coding: ISO-8859-1 -*-
   3 #
   4 #
   5 # Copyright (C) 2003-2006 Michael Schindler <m-schindler@users.sourceforge.net>
   6 # Copyright (C) 2003-2005 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, mathutils, path, normpath, unit, color
  27 from path import degrees, radians
  28
  29 # specific exception for an invalid parameterization point
  30 # used in parallel
  31 class InvalidParamException(Exception):
  32
  33     def __init__(self, param):
  34         self.normsubpathitemparam = param
  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, epsilon, 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, epsilon=epsilon)[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 deformer(attr.attr):
 315
 316     def deform (self, basepath):
 317         return basepath
 318
 319 class cycloid(deformer): # <<<
 320     """Wraps a cycloid around a path.
 321
 322     The outcome looks like a spring with the originalpath as the axis.
 323     radius: radius of the cycloid
 324     halfloops:  number of halfloops
 325     skipfirst/skiplast: undeformed end lines of the original path
 326     curvesperhloop:
 327     sign: start left (1) or right (-1) with the first halfloop
 328     turnangle: angle of perspective on a (3D) spring
 329                turnangle=0 will produce a sinus-like cycloid,
 330                turnangle=90 will procude a row of connected circles
 331
 332     """
 333
 334     def __init__(self, radius=0.5*unit.t_cm, halfloops=10,
 335     skipfirst=1*unit.t_cm, skiplast=1*unit.t_cm, curvesperhloop=3, sign=1, turnangle=45):
 336         self.skipfirst = skipfirst
 337         self.skiplast = skiplast
 338         self.radius = radius
 339         self.halfloops = halfloops
 340         self.curvesperhloop = curvesperhloop
 341         self.sign = sign
 342         self.turnangle = turnangle
 343
 344     def __call__(self, radius=None, halfloops=None,
 345     skipfirst=None, skiplast=None, curvesperhloop=None, sign=None, turnangle=None):
 346         if radius is None:
 347             radius = self.radius
 348         if halfloops is None:
 349             halfloops = self.halfloops
 350         if skipfirst is None:
 351             skipfirst = self.skipfirst
 352         if skiplast is None:
 353             skiplast = self.skiplast
 354         if curvesperhloop is None:
 355             curvesperhloop = self.curvesperhloop
 356         if sign is None:
 357             sign = self.sign
 358         if turnangle is None:
 359             turnangle = self.turnangle
 360
 361         return cycloid(radius=radius, halfloops=halfloops, skipfirst=skipfirst, skiplast=skiplast,
 362                        curvesperhloop=curvesperhloop, sign=sign, turnangle=turnangle)
 363
 364     def deform(self, basepath):
 365         resultnormsubpaths = [self.deformsubpath(nsp) for nsp in basepath.normpath().normsubpaths]
 366         return normpath.normpath(resultnormsubpaths)
 367
 368     def deformsubpath(self, normsubpath):
 369
 370         skipfirst = abs(unit.topt(self.skipfirst))
 371         skiplast = abs(unit.topt(self.skiplast))
 372         radius = abs(unit.topt(self.radius))
 373         turnangle = degrees(self.turnangle)
 374         sign = mathutils.sign(self.sign)
 375
 376         cosTurn = math.cos(turnangle)
 377         sinTurn = math.sin(turnangle)
 378
 379         # make list of the lengths and parameters at points on normsubpath
 380         # where we will add cycloid-points
 381         totlength = normsubpath.arclen_pt()
 382         if totlength <= skipfirst + skiplast + 2*radius*sinTurn:
 383             warnings.warn("normsubpath is too short for deformation with cycloid -- skipping...")
 384             return normsubpath
 385
 386         # parameterization is in rotation-angle around the basepath
 387         # differences in length, angle ... between two basepoints
 388         # and between basepoints and controlpoints
 389         Dphi = math.pi / self.curvesperhloop
 390         phis = [i * Dphi for i in range(self.halfloops * self.curvesperhloop + 1)]
 391         DzDphi = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * math.pi * cosTurn)
 392         # Dz = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * self.curvesperhloop * cosTurn)
 393         # zs = [i * Dz for i in range(self.halfloops * self.curvesperhloop + 1)]
 394         # from path._arctobcurve:
 395         # optimal relative distance along tangent for second and third control point
 396         L = 4 * radius * (1 - math.cos(Dphi/2)) / (3 * math.sin(Dphi/2))
 397
 398         # Now the transformation of z into the turned coordinate system
 399         Zs = [ skipfirst + radius*sinTurn # here the coordinate z starts
 400              - sinTurn*radius*math.cos(phi) + cosTurn*DzDphi*phi # the transformed z-coordinate
 401              for phi in phis]
 402         params = normsubpath._arclentoparam_pt(Zs)[0]
 403
 404         # get the positions of the splitpoints in the cycloid
 405         points = []
 406         for phi, param in zip(phis, params):
 407             # the cycloid is a circle that is stretched along the normsubpath
 408             # here are the points of that circle
 409             basetrafo = normsubpath.trafo([param])[0]
 410
 411             # The point on the cycloid, in the basepath's local coordinate system
 412             baseZ, baseY = 0, radius*math.sin(phi)
 413
 414             # The tangent there, also in local coords
 415             tangentX = -cosTurn*radius*math.sin(phi) + sinTurn*DzDphi
 416             tangentY = radius*math.cos(phi)
 417             tangentZ = sinTurn*radius*math.sin(phi) + DzDphi*cosTurn
 418             norm = math.sqrt(tangentX*tangentX + tangentY*tangentY + tangentZ*tangentZ)
 419             tangentY, tangentZ = tangentY/norm, tangentZ/norm
 420
 421             # Respect the curvature of the basepath for the cycloid's curvature
 422             # XXX this is only a heuristic, not a "true" expression for
 423             #     the curvature in curved coordinate systems
 424             pathradius = normsubpath.curveradius_pt([param])[0]
 425             if pathradius is not normpath.invalid:
 426                 factor = (pathradius - baseY) / pathradius
 427                 factor = abs(factor)
 428             else:
 429                 factor = 1
 430             l = L * factor
 431
 432             # The control points prior and after the point on the cycloid
 433             preeZ, preeY = baseZ - l * tangentZ, baseY - l * tangentY
 434             postZ, postY = baseZ + l * tangentZ, baseY + l * tangentY
 435
 436             # Now put everything at the proper place
 437             points.append(basetrafo.apply_pt(preeZ, sign * preeY) +
 438                           basetrafo.apply_pt(baseZ, sign * baseY) +
 439                           basetrafo.apply_pt(postZ, sign * postY))
 440
 441         if len(points) <= 1:
 442             warnings.warn("normsubpath is too short for deformation with cycloid -- skipping...")
 443             return normsubpath
 444
 445         # Build the path from the pointlist
 446         # containing (control x 2,  base x 2, control x 2)
 447         if skipfirst > normsubpath.epsilon:
 448             normsubpathitems = normsubpath.segments([0, params[0]])[0]
 449             normsubpathitems.append(normpath.normcurve_pt(*(points[0][2:6] + points[1][0:4])))
 450         else:
 451             normsubpathitems = [normpath.normcurve_pt(*(points[0][2:6] + points[1][0:4]))]
 452         for i in range(1, len(points)-1):
 453             normsubpathitems.append(normpath.normcurve_pt(*(points[i][2:6] + points[i+1][0:4])))
 454         if skiplast > normsubpath.epsilon:
 455             for nsp in normsubpath.segments([params[-1], len(normsubpath)]):
 456                 normsubpathitems.extend(nsp.normsubpathitems)
 457
 458         # That's it
 459         return normpath.normsubpath(normsubpathitems, epsilon=normsubpath.epsilon)
 460 # >>>
 461
 462 cycloid.clear = attr.clearclass(cycloid)
 463
 464 class smoothed(deformer): # <<<
 465
 466     """Bends corners in a normpath.
 467
 468     This decorator replaces corners in a normpath with bezier curves. There are two cases:
 469     - If the corner lies between two lines, _two_ bezier curves will be used
 470       that are highly optimized to look good (their curvature is to be zero at the ends
 471       and has to have zero derivative in the middle).
 472       Additionally, it can controlled by the softness-parameter.
 473     - If the corner lies between curves then _one_ bezier is used that is (except in some
 474       special cases) uniquely determined by the tangents and curvatures at its end-points.
 475       In some cases it is necessary to use only the absolute value of the curvature to avoid a
 476       cusp-shaped connection of the new bezier to the old path. In this case the use of
 477       "obeycurv=0" allows the sign of the curvature to switch.
 478     - The radius argument gives the arclength-distance of the corner to the points where the
 479       old path is cut and the beziers are inserted.
 480     - Path elements that are too short (shorter than the radius) are skipped
 481     """
 482
 483     def __init__(self, radius, softness=1, obeycurv=0, relskipthres=0.01):
 484         self.radius = radius
 485         self.softness = softness
 486         self.obeycurv = obeycurv
 487         self.relskipthres = relskipthres
 488
 489     def __call__(self, radius=None, softness=None, obeycurv=None, relskipthres=None):
 490         if radius is None:
 491             radius = self.radius
 492         if softness is None:
 493             softness = self.softness
 494         if obeycurv is None:
 495             obeycurv = self.obeycurv
 496         if relskipthres is None:
 497             relskipthres = self.relskipthres
 498         return smoothed(radius=radius, softness=softness, obeycurv=obeycurv, relskipthres=relskipthres)
 499
 500     def deform(self, basepath):
 501         return normpath.normpath([self.deformsubpath(normsubpath)
 502                               for normsubpath in basepath.normpath().normsubpaths])
 503
 504     def deformsubpath(self, normsubpath):
 505         radius_pt = unit.topt(self.radius)
 506         epsilon = normsubpath.epsilon
 507
 508         # remove too short normsubpath items (shorter than self.relskipthres*radius_pt or epsilon)
 509         pertinentepsilon = max(epsilon, self.relskipthres*radius_pt)
 510         pertinentnormsubpath = normpath.normsubpath(normsubpath.normsubpathitems,
 511                                                 epsilon=pertinentepsilon)
 512         pertinentnormsubpath.flushskippedline()
 513         pertinentnormsubpathitems = pertinentnormsubpath.normsubpathitems
 514
 515         # calculate the splitting parameters for the pertinentnormsubpathitems
 516         arclens_pt = []
 517         params = []
 518         for pertinentnormsubpathitem in pertinentnormsubpathitems:
 519             arclen_pt = pertinentnormsubpathitem.arclen_pt(epsilon)
 520             arclens_pt.append(arclen_pt)
 521             l1_pt = min(radius_pt, 0.5*arclen_pt)
 522             l2_pt = max(0.5*arclen_pt, arclen_pt - radius_pt)
 523             params.append(pertinentnormsubpathitem.arclentoparam_pt([l1_pt, l2_pt], epsilon))
 524
 525         # handle the first and last pertinentnormsubpathitems for a non-closed normsubpath
 526         if not normsubpath.closed:
 527             l1_pt = 0
 528             l2_pt = max(0, arclens_pt[0] - radius_pt)
 529             params[0] = pertinentnormsubpathitems[0].arclentoparam_pt([l1_pt, l2_pt], epsilon)
 530             l1_pt = min(radius_pt, arclens_pt[-1])
 531             l2_pt = arclens_pt[-1]
 532             params[-1] = pertinentnormsubpathitems[-1].arclentoparam_pt([l1_pt, l2_pt], epsilon)
 533
 534         newnormsubpath = normpath.normsubpath(epsilon=normsubpath.epsilon)
 535         for i in range(len(pertinentnormsubpathitems)):
 536             this = i
 537             next = (i+1) % len(pertinentnormsubpathitems)
 538             thisparams = params[this]
 539             nextparams = params[next]
 540             thisnormsubpathitem = pertinentnormsubpathitems[this]
 541             nextnormsubpathitem = pertinentnormsubpathitems[next]
 542             thisarclen_pt = arclens_pt[this]
 543             nextarclen_pt = arclens_pt[next]
 544
 545             # insert the middle segment
 546             newnormsubpath.append(thisnormsubpathitem.segments(thisparams)[0])
 547
 548             # insert replacement curves for the corners
 549             if next or normsubpath.closed:
 550
 551                 t1 = thisnormsubpathitem.rotation([thisparams[1]])[0].apply_pt(1, 0)
 552                 t2 = nextnormsubpathitem.rotation([nextparams[0]])[0].apply_pt(1, 0)
 553                 # TODO: normpath.invalid
 554
 555                 if (isinstance(thisnormsubpathitem, normpath.normline_pt) and
 556                     isinstance(nextnormsubpathitem, normpath.normline_pt)):
 557
 558                     # case of two lines -> replace by two curves
 559                     d1, g1, f1, e, f2, g2, d2 = curvescontrols_from_endlines_pt(
 560                         thisnormsubpathitem.atend_pt(), t1, t2,
 561                         thisarclen_pt*(1-thisparams[1]), nextarclen_pt*(nextparams[0]), softness=self.softness)
 562
 563                     p1 = thisnormsubpathitem.at_pt([thisparams[1]])[0]
 564                     p2 = nextnormsubpathitem.at_pt([nextparams[0]])[0]
 565
 566                     newnormsubpath.append(normpath.normcurve_pt(*(d1 + g1 + f1 + e)))
 567                     newnormsubpath.append(normpath.normcurve_pt(*(e + f2 + g2 + d2)))
 568
 569                 else:
 570
 571                     # generic case -> replace by a single curve with prescribed tangents and curvatures
 572                     p1 = thisnormsubpathitem.at_pt([thisparams[1]])[0]
 573                     p2 = nextnormsubpathitem.at_pt([nextparams[0]])[0]
 574                     c1 = thisnormsubpathitem.curvature_pt([thisparams[1]])[0]
 575                     c2 = nextnormsubpathitem.curvature_pt([nextparams[0]])[0]
 576                     # TODO: normpath.invalid
 577
 578                     # TODO: more intelligent fallbacks:
 579                     #   circle -> line
 580                     #   circle -> circle
 581
 582                     if not self.obeycurv:
 583                         # do not obey the sign of the curvature but
 584                         # make the sign such that the curve smoothly passes to the next point
 585                         # this results in a discontinuous curvature
 586                         # (but the absolute value is still continuous)
 587                         s1 = +mathutils.sign(t1[0] * (p2[1]-p1[1]) - t1[1] * (p2[0]-p1[0]))
 588                         s2 = -mathutils.sign(t2[0] * (p2[1]-p1[1]) - t2[1] * (p2[0]-p1[0]))
 589                         c1 = s1 * abs(c1)
 590                         c2 = s2 * abs(c2)
 591
 592                     # get the length of the control "arms"
 593                     controldists = controldists_from_endgeometry_pt(p1, p2, t1, t2, c1, c2, epsilon=epsilon)
 594
 595                     if controldists and (controldists[0][0] >= 0 and controldists[0][1] >= 0):
 596                         # use the first entry in the controldists
 597                         # this should be the "smallest" pair
 598                         a, d = controldists[0]
 599                         # avoid curves with invalid parameterization
 600                         a = max(a, epsilon)
 601                         d = max(d, epsilon)
 602
 603                         # avoid overshooting at the corners:
 604                         # this changes not only the sign of the curvature
 605                         # but also the magnitude
 606                         if not self.obeycurv:
 607                             t, s = intersection(p1, p2, t1, t2)
 608                             if (t is not None and s is not None and
 609                                 t > 0 and s < 0):
 610                                 a = min(a, abs(t))
 611                                 d = min(d, abs(s))
 612
 613                     else:
 614                         # use a fallback
 615                         t, s = intersection(p1, p2, t1, t2)
 616                         if t is not None and s is not None:
 617                             a = 0.65 * abs(t)
 618                             d = 0.65 * abs(s)
 619                         else:
 620                             # if there is no useful result:
 621                             # take an arbitrary smoothing curve that does not obey
 622                             # the curvature constraints
 623                             dist = math.hypot(p1[0] - p2[0], p1[1] - p2[1])
 624                             a = dist / (3.0 * math.hypot(*t1))
 625                             d = dist / (3.0 * math.hypot(*t2))
 626
 627                     # calculate the two missing control points
 628                     q1 = p1[0] + a * t1[0], p1[1] + a * t1[1]
 629                     q2 = p2[0] - d * t2[0], p2[1] - d * t2[1]
 630
 631                     newnormsubpath.append(normpath.normcurve_pt(*(p1 + q1 + q2 + p2)))
 632
 633         if normsubpath.closed:
 634             newnormsubpath.close()
 635         return newnormsubpath
 636
 637 # >>>
 638
 639 smoothed.clear = attr.clearclass(smoothed)
 640
 641 class parallel(deformer): # <<<
 642
 643     """creates a parallel normpath with constant distance to the original normpath
 644
 645     A positive 'distance' results in a curve left of the original one -- and a
 646     negative 'distance' in a curve at the right. Left/Right are understood in
 647     terms of the parameterization of the original curve. For each path element
 648     a parallel curve/line is constructed. At corners, either a circular arc is
 649     drawn around the corner, or, if possible, the parallel curve is cut in
 650     order to also exhibit a corner.
 651
 652     distance:            the distance of the parallel normpath
 653     relerr:              distance*relerr is the maximal allowed error in the parallel distance
 654     sharpoutercorners:   make the outer corners not round but sharp.
 655                          The inner corners (corners after inflection points) will stay round
 656     dointersection:      boolean for doing the intersection step (default: 1).
 657                          Set this value to 0 if you want the whole parallel path
 658     checkdistanceparams: a list of parameter values in the interval (0,1) where the
 659                          parallel distance is checked on each normpathitem
 660     lookforcurvatures:   number of points per normpathitem where is looked for
 661                          a critical value of the curvature
 662     """
 663
 664     # TODO:
 665     # * do testing for curv=0, T=0, D=0, E=0 cases
 666     # * do testing for several random curves
 667     #   -- does the recursive deformnicecurve converge?
 668
 669
 670     def __init__(self, distance, relerr=0.05, sharpoutercorners=0, dointersection=1,
 671                        checkdistanceparams=[0.5], lookforcurvatures=11, debug=None):
 672         self.distance = distance
 673         self.relerr = relerr
 674         self.sharpoutercorners = sharpoutercorners
 675         self.checkdistanceparams = checkdistanceparams
 676         self.lookforcurvatures = lookforcurvatures
 677         self.dointersection = dointersection
 678         self.debug = debug
 679
 680     def __call__(self, distance=None, relerr=None, sharpoutercorners=None, dointersection=None,
 681                        checkdistanceparams=None, lookforcurvatures=None, debug=None):
 682         # returns a copy of the deformer with different parameters
 683         if distance is None:
 684             distance = self.distance
 685         if relerr is None:
 686             relerr = self.relerr
 687         if sharpoutercorners is None:
 688             sharpoutercorners = self.sharpoutercorners
 689         if dointersection is None:
 690             dointersection = self.dointersection
 691         if checkdistanceparams is None:
 692             checkdistanceparams = self.checkdistanceparams
 693         if lookforcurvatures is None:
 694             lookforcurvatures = self.lookforcurvatures
 695         if debug is None:
 696             debug = self.debug
 697
 698         return parallel(distance=distance, relerr=relerr,
 699                         sharpoutercorners=sharpoutercorners,
 700                         dointersection=dointersection,
 701                         checkdistanceparams=checkdistanceparams,
 702                         lookforcurvatures=lookforcurvatures,
 703                         debug=debug)
 704
 705     def deform(self, basepath):
 706         self.dist_pt = unit.topt(self.distance)
 707         resultnormsubpaths = []
 708         for nsp in basepath.normpath().normsubpaths:
 709             parallel_normpath = self.deformsubpath(nsp)
 710             resultnormsubpaths += parallel_normpath.normsubpaths
 711         result = normpath.normpath(resultnormsubpaths)
 712         return result
 713
 714     def deformsubpath(self, orig_nsp): # <<<
 715
 716         """returns a list of normsubpaths building the parallel curve"""
 717
 718         dist = self.dist_pt
 719         epsilon = orig_nsp.epsilon
 720
 721         # avoid too small dists: we would run into instabilities
 722         if abs(dist) < abs(epsilon):
 723             return orig_nsp
 724
 725         result = normpath.normpath()
 726
 727         # iterate over the normsubpath in the following manner:
 728         # * for each item first append the additional arc / intersect
 729         #   and then add the next parallel piece
 730         # * for the first item only add the parallel piece
 731         #   (because this is done for next_orig_nspitem, we need to start with next=0)
 732         for i in range(len(orig_nsp.normsubpathitems)):
 733             prev = i-1
 734             next = i
 735             prev_orig_nspitem = orig_nsp.normsubpathitems[prev]
 736             next_orig_nspitem = orig_nsp.normsubpathitems[next]
 737
 738             stepsize = 0.01
 739             prev_param, prev_rotation = self.valid_near_rotation(prev_orig_nspitem, 1, 0, stepsize, 0.5*epsilon)
 740             next_param, next_rotation = self.valid_near_rotation(next_orig_nspitem, 0, 1, stepsize, 0.5*epsilon)
 741             # TODO: eventually shorten next_orig_nspitem
 742
 743             prev_tangent = prev_rotation.apply_pt(1, 0)
 744             next_tangent = next_rotation.apply_pt(1, 0)
 745
 746             # get the next parallel piece for the normpath
 747             try:
 748                 next_parallel_normpath = self.deformsubpathitem(next_orig_nspitem, epsilon)
 749             except InvalidParamException, e:
 750                 invalid_nspitem_param = e.normsubpathitemparam
 751                 # split the nspitem apart and continue with pieces that do not contain
 752                 # the invalid point anymore. At the end, simply take one piece, otherwise two.
 753                 stepsize = 0.01
 754                 if self.length_pt(next_orig_nspitem, invalid_nspitem_param, 0) > epsilon:
 755                     if self.length_pt(next_orig_nspitem, invalid_nspitem_param, 1) > epsilon:
 756                         p1, foo = self.valid_near_rotation(next_orig_nspitem, invalid_nspitem_param, 0, stepsize, 0.5*epsilon)
 757                         p2, foo = self.valid_near_rotation(next_orig_nspitem, invalid_nspitem_param, 1, stepsize, 0.5*epsilon)
 758                         segments = next_orig_nspitem.segments([0, p1, p2, 1])[0:3:2]
 759                         segments[1] = segments[1].modifiedbegin_pt(*(segments[0].atend_pt()))
 760                     else:
 761                         p1, foo = self.valid_near_rotation(next_orig_nspitem, invalid_nspitem_param, 0, stepsize, 0.5*epsilon)
 762                         segments = next_orig_nspitem.segments([0, p1])
 763                 else:
 764                     p2, foo = self.valid_near_rotation(next_orig_nspitem, invalid_nspitem_param, 1, stepsize, 0.5*epsilon)
 765                     segments = next_orig_nspitem.segments([p2, 1])
 766
 767                 next_parallel_normpath = self.deformsubpath(normpath.normsubpath(segments, epsilon=epsilon))
 768
 769             if not (next_parallel_normpath.normsubpaths and next_parallel_normpath[0].normsubpathitems):
 770                 continue
 771
 772             # this starts the whole normpath
 773             if not result.normsubpaths:
 774                 result = next_parallel_normpath
 775                 continue
 776
 777             # sinus of the angle between the tangents
 778             # sinangle > 0 for a left-turning nexttangent
 779             # sinangle < 0 for a right-turning nexttangent
 780             sinangle = prev_tangent[0]*next_tangent[1] - prev_tangent[1]*next_tangent[0]
 781             cosangle = prev_tangent[0]*next_tangent[0] + prev_tangent[1]*next_tangent[1]
 782             if cosangle < 0 or abs(dist*math.asin(sinangle)) >= epsilon:
 783                 if self.sharpoutercorners and dist*sinangle < 0:
 784                     A1, A2 = result.atend_pt(), next_parallel_normpath.atbegin_pt()
 785                     t1, t2 = intersection(A1, A2, prev_tangent, next_tangent)
 786                     B = A1[0] + t1 * prev_tangent[0], A1[1] + t1 * prev_tangent[1]
 787                     arc_normpath = normpath.normpath([normpath.normsubpath([
 788                         normpath.normline_pt(A1[0], A1[1], B[0], B[1]),
 789                         normpath.normline_pt(B[0], B[1], A2[0], A2[1])
 790                         ])])
 791                 else:
 792                     # We must append an arc around the corner
 793                     arccenter = next_orig_nspitem.atbegin_pt()
 794                     arcbeg = result.atend_pt()
 795                     arcend = next_parallel_normpath.atbegin_pt()
 796                     angle1 = math.atan2(arcbeg[1] - arccenter[1], arcbeg[0] - arccenter[0])
 797                     angle2 = math.atan2(arcend[1] - arccenter[1], arcend[0] - arccenter[0])
 798
 799                     # depending on the direction we have to use arc or arcn
 800                     if dist > 0:
 801                         arcclass = path.arcn_pt
 802                     else:
 803                         arcclass = path.arc_pt
 804                     arc_normpath = path.path(arcclass(
 805                       arccenter[0], arccenter[1], abs(dist),
 806                       degrees(angle1), degrees(angle2))).normpath(epsilon=epsilon)
 807
 808                 # append the arc to the parallel path
 809                 result.join(arc_normpath)
 810                 # append the next parallel piece to the path
 811                 result.join(next_parallel_normpath)
 812             else:
 813                 # The path is quite straight between prev and next item:
 814                 # normpath.normpath.join adds a straight line if necessary
 815                 result.join(next_parallel_normpath)
 816
 817
 818         # end here if nothing has been found so far
 819         if not (result.normsubpaths and result[-1].normsubpathitems):
 820             return result
 821
 822         # the curve around the closing corner may still be missing
 823         if orig_nsp.closed:
 824             # TODO: normpath.invalid
 825             stepsize = 0.01
 826             prev_param, prev_rotation = self.valid_near_rotation(result[-1][-1], 1, 0, stepsize, 0.5*epsilon)
 827             next_param, next_rotation = self.valid_near_rotation(result[0][0], 0, 1, stepsize, 0.5*epsilon)
 828             # TODO: eventually shorten next_orig_nspitem
 829
 830             prev_tangent = prev_rotation.apply_pt(1, 0)
 831             next_tangent = next_rotation.apply_pt(1, 0)
 832             sinangle = prev_tangent[0]*next_tangent[1] - prev_tangent[1]*next_tangent[0]
 833             cosangle = prev_tangent[0]*next_tangent[0] + prev_tangent[1]*next_tangent[1]
 834
 835             if cosangle < 0 or abs(dist*math.asin(sinangle)) >= epsilon:
 836                 # We must append an arc around the corner
 837                 # TODO: avoid the code dublication
 838                 if self.sharpoutercorners and dist*sinangle < 0:
 839                     A1, A2 = result.atend_pt(), result.atbegin_pt()
 840                     t1, t2 = intersection(A1, A2, prev_tangent, next_tangent)
 841                     B = A1[0] + t1 * prev_tangent[0], A1[1] + t1 * prev_tangent[1]
 842                     arc_normpath = normpath.normpath([normpath.normsubpath([
 843                         normpath.normline_pt(A1[0], A1[1], B[0], B[1]),
 844                         normpath.normline_pt(B[0], B[1], A2[0], A2[1])
 845                         ])])
 846                 else:
 847                     arccenter = orig_nsp.atend_pt()
 848                     arcbeg = result.atend_pt()
 849                     arcend = result.atbegin_pt()
 850                     angle1 = math.atan2(arcbeg[1] - arccenter[1], arcbeg[0] - arccenter[0])
 851                     angle2 = math.atan2(arcend[1] - arccenter[1], arcend[0] - arccenter[0])
 852
 853                     # depending on the direction we have to use arc or arcn
 854                     if dist > 0:
 855                         arcclass = path.arcn_pt
 856                     else:
 857                         arcclass = path.arc_pt
 858                     arc_normpath = path.path(arcclass(
 859                         arccenter[0], arccenter[1], abs(dist),
 860                         degrees(angle1), degrees(angle2))).normpath(epsilon=epsilon)
 861
 862                 # append the arc to the parallel path
 863                 if (result.normsubpaths and result[-1].normsubpathitems and
 864                     arc_normpath.normsubpaths and arc_normpath[-1].normsubpathitems):
 865                     result.join(arc_normpath)
 866
 867             if len(result) == 1:
 868                 result[0].close()
 869             else:
 870                 # if the parallel normpath is split into several subpaths anyway,
 871                 # then use the natural beginning and ending
 872                 # closing is not possible anymore
 873                 for nspitem in result[0]:
 874                     result[-1].append(nspitem)
 875                 result.normsubpaths = result.normsubpaths[1:]
 876
 877         if self.dointersection:
 878             result = self.rebuild_intersected_normpath(result, normpath.normpath([orig_nsp]), epsilon)
 879
 880         return result
 881         # >>>
 882     def deformsubpathitem(self, nspitem, epsilon): # <<<
 883
 884         """Returns a parallel normpath for a single normsubpathitem
 885
 886         Analyzes the curvature of a normsubpathitem and returns a normpath with
 887         the appropriate number of normsubpaths. This must be a normpath because
 888         a normcurve can be strongly curved, such that the parallel path must
 889         contain a hole"""
 890
 891         dist = self.dist_pt
 892
 893         # for a simple line we return immediately
 894         if isinstance(nspitem, normpath.normline_pt):
 895             normal = nspitem.rotation([0])[0].apply_pt(0, 1)
 896             start = nspitem.atbegin_pt()
 897             end = nspitem.atend_pt()
 898             return path.line_pt(
 899                 start[0] + dist * normal[0], start[1] + dist * normal[1],
 900                 end[0] + dist * normal[0], end[1] + dist * normal[1]).normpath(epsilon=epsilon)
 901
 902         # for a curve we have to check if the curvatures
 903         # cross the singular value 1/dist
 904         crossings = self.distcrossingparameters(nspitem, epsilon)
 905
 906         # depending on the number of crossings we must consider
 907         # three different cases:
 908         if crossings:
 909             # The curvature crosses the borderline 1/dist
 910             # the parallel curve contains points with infinite curvature!
 911             result = normpath.normpath()
 912
 913             # we need the endpoints of the nspitem
 914             if self.length_pt(nspitem, crossings[0], 0) > epsilon:
 915                 crossings.insert(0, 0)
 916             if self.length_pt(nspitem, crossings[-1], 1) > epsilon:
 917                 crossings.append(1)
 918
 919             for i in range(len(crossings) - 1):
 920                 middleparam = 0.5*(crossings[i] + crossings[i+1])
 921                 middlecurv = nspitem.curvature_pt([middleparam])[0]
 922                 if middlecurv is normpath.invalid:
 923                     raise InvalidParamException(middleparam)
 924                 # the radius is good if
 925                 #  - middlecurv and dist have opposite signs or
 926                 #  - middlecurv is "smaller" than 1/dist
 927                 if middlecurv*dist < 0 or abs(dist*middlecurv) < 1:
 928                     parallel_nsp = self.deformnicecurve(nspitem.segments(crossings[i:i+2])[0], epsilon)
 929                     # never append empty normsubpaths
 930                     if parallel_nsp.normsubpathitems:
 931                         result.append(parallel_nsp)
 932
 933             return result
 934
 935         else:
 936             # the curvature is either bigger or smaller than 1/dist
 937             middlecurv = nspitem.curvature_pt([0.5])[0]
 938             if dist*middlecurv < 0 or abs(dist*middlecurv) < 1:
 939                 # The curve is everywhere less curved than 1/dist
 940                 # We can proceed finding the parallel curve for the whole piece
 941                 parallel_nsp = self.deformnicecurve(nspitem, epsilon)
 942                 # never append empty normsubpaths
 943                 if parallel_nsp.normsubpathitems:
 944                     return normpath.normpath([parallel_nsp])
 945                 else:
 946                     return normpath.normpath()
 947             else:
 948                 # the curve is everywhere stronger curved than 1/dist
 949                 # There is nothing to be returned.
 950                 return normpath.normpath()
 951
 952         # >>>
 953     def deformnicecurve(self, normcurve, epsilon, startparam=0.0, endparam=1.0): # <<<
 954
 955         """Returns a parallel normsubpath for the normcurve.
 956
 957         This routine assumes that the normcurve is everywhere
 958         'less' curved than 1/dist and contains no point with an
 959         invalid parameterization
 960         """
 961         dist = self.dist_pt
 962         T_threshold = 1.0e-5
 963
 964         # normalized tangent directions
 965         tangA, tangD = normcurve.rotation([startparam, endparam])
 966         # if we find an unexpected normpath.invalid we have to
 967         # parallelise this normcurve on the level of split normsubpaths
 968         if tangA is normpath.invalid:
 969             raise InvalidParamException(startparam)
 970         if tangD is normpath.invalid:
 971             raise InvalidParamException(endparam)
 972         tangA = tangA.apply_pt(1, 0)
 973         tangD = tangD.apply_pt(1, 0)
 974
 975         # the new starting points
 976         orig_A, orig_D = normcurve.at_pt([startparam, endparam])
 977         A = orig_A[0] - dist * tangA[1], orig_A[1] + dist * tangA[0]
 978         D = orig_D[0] - dist * tangD[1], orig_D[1] + dist * tangD[0]
 979
 980         # we need to end this _before_ we will run into epsilon-problems
 981         # when creating curves we do not want to calculate the length of
 982         # or even split it for recursive calls
 983         if (math.hypot(A[0] - D[0], A[1] - D[1]) < epsilon and
 984             math.hypot(tangA[0] - tangD[0], tangA[1] - tangD[1]) < T_threshold):
 985             return normpath.normsubpath([normpath.normline_pt(A[0], A[1], D[0], D[1])], epsilon=epsilon)
 986
 987         result = normpath.normsubpath(epsilon=epsilon)
 988         # is there enough space on the normals before they intersect?
 989         a, d = intersection(orig_A, orig_D, (-tangA[1], tangA[0]), (-tangD[1], tangD[0]))
 990         # a,d are the lengths to the intersection points:
 991         # for a (and equally for b) we can proceed in one of the following cases:
 992         #   a is None (means parallel normals)
 993         #   a and dist have opposite signs (and the same for b)
 994         #   a has the same sign but is bigger
 995         if ( (a is None or a*dist < 0 or abs(a) > abs(dist) + epsilon) or
 996              (d is None or d*dist < 0 or abs(d) > abs(dist) + epsilon) ):
 997             # the original path is long enough to draw a parallel piece
 998             # this is the generic case. Get the parallel curves
 999             orig_curvA, orig_curvD = normcurve.curvature_pt([startparam, endparam])
1000             # normpath.invalid may not appear here because we have asked
1001             # for this already at the tangents
1002             assert orig_curvA is not normpath.invalid
1003             assert orig_curvD is not normpath.invalid
1004             curvA = orig_curvA / (1.0 - dist*orig_curvA)
1005             curvD = orig_curvD / (1.0 - dist*orig_curvD)
1006
1007             # first try to approximate the normcurve with a single item
1008             controldistpairs = controldists_from_endgeometry_pt(A, D, tangA, tangD, curvA, curvD, epsilon=epsilon)
1009
1010             if controldistpairs:
1011                 # TODO: is it good enough to get the first entry here?
1012                 #       from testing: this fails if there are loops in the original curve
1013                 a, d = controldistpairs[0]
1014                 if a >= 0 and d >= 0:
1015                     if a < epsilon and d < epsilon:
1016                         result = normpath.normsubpath([normpath.normline_pt(A[0], A[1], D[0], D[1])], epsilon=epsilon)
1017                     else:
1018                         # we avoid curves with invalid parameterization
1019                         a = max(a, epsilon)
1020                         d = max(d, epsilon)
1021                         result = normpath.normsubpath([normpath.normcurve_pt(
1022                             A[0], A[1],
1023                             A[0] + a * tangA[0], A[1] + a * tangA[1],
1024                             D[0] - d * tangD[0], D[1] - d * tangD[1],
1025                             D[0], D[1])], epsilon=epsilon)
1026
1027             # then try with two items, recursive call
1028             if ((not result.normsubpathitems) or
1029                 (self.checkdistanceparams and result.normsubpathitems
1030                  and not self.distchecked(normcurve, result, epsilon, startparam, endparam))):
1031                 # TODO: does this ever converge?
1032                 # TODO: what if this hits epsilon?
1033                 firstnsp = self.deformnicecurve(normcurve, epsilon, startparam, 0.5*(startparam+endparam))
1034                 secondnsp = self.deformnicecurve(normcurve, epsilon, 0.5*(startparam+endparam), endparam)
1035                 if not (firstnsp.normsubpathitems and secondnsp.normsubpathitems):
1036                     result = normpath.normsubpath(
1037                         [normpath.normline_pt(A[0], A[1], D[0], D[1])], epsilon=epsilon)
1038                 else:
1039                     # we will get problems if the curves are too short:
1040                     result = firstnsp.joined(secondnsp)
1041
1042         return result
1043         # >>>
1044
1045     def distchecked(self, orig_normcurve, parallel_normsubpath, epsilon, tstart, tend): # <<<
1046
1047         """Checks the distances between orig_normcurve and parallel_normsubpath
1048
1049         The checking is done at parameters self.checkdistanceparams of orig_normcurve."""
1050
1051         dist = self.dist_pt
1052         # do not look closer than epsilon:
1053         dist_relerr = mathutils.sign(dist) * max(abs(self.relerr*dist), epsilon)
1054
1055         checkdistanceparams = [tstart + (tend-tstart)*t for t in self.checkdistanceparams]
1056
1057         for param, P, rotation in zip(checkdistanceparams,
1058                                       orig_normcurve.at_pt(checkdistanceparams),
1059                                       orig_normcurve.rotation(checkdistanceparams)):
1060             # check if the distance is really the wanted distance
1061             # measure the distance in the "middle" of the original curve
1062             if rotation is normpath.invalid:
1063                 raise InvalidParamException(param)
1064
1065             normal = rotation.apply_pt(0, 1)
1066
1067             # create a short cutline for intersection only:
1068             cutline = normpath.normsubpath([normpath.normline_pt (
1069               P[0] + (dist - 2*dist_relerr) * normal[0],
1070               P[1] + (dist - 2*dist_relerr) * normal[1],
1071               P[0] + (dist + 2*dist_relerr) * normal[0],
1072               P[1] + (dist + 2*dist_relerr) * normal[1])], epsilon=epsilon)
1073
1074             cutparams = parallel_normsubpath.intersect(cutline)
1075             distances = [math.hypot(P[0] - cutpoint[0], P[1] - cutpoint[1])
1076                          for cutpoint in cutline.at_pt(cutparams[1])]
1077
1078             if (not distances) or (abs(min(distances) - abs(dist)) > abs(dist_relerr)):
1079                 return 0
1080
1081         return 1
1082     # >>>
1083     def distcrossingparameters(self, normcurve, epsilon, tstart=0, tend=1): # <<<
1084
1085         """Returns a list of parameters where the curvature is 1/distance"""
1086
1087         dist = self.dist_pt
1088
1089         # we _need_ to do this with the curvature, not with the radius
1090         # because the curvature is continuous at the straight line and the radius is not:
1091         # when passing from one slightly curved curve to the other with opposite curvature sign,
1092         # via the straight line, then the curvature changes its sign at curv=0, while the
1093         # radius changes its sign at +/-infinity
1094         # this causes instabilities for nearly straight curves
1095
1096         # include tstart and tend
1097         params = [tstart + i * (tend - tstart) * 1.0 / (self.lookforcurvatures - 1)
1098                   for i in range(self.lookforcurvatures)]
1099         curvs = normcurve.curvature_pt(params)
1100
1101         # break everything at invalid curvatures
1102         for param, curv in zip(params, curvs):
1103             if curv is normpath.invalid:
1104                 raise InvalidParamException(param)
1105
1106         parampairs = zip(params[:-1], params[1:])
1107         curvpairs = zip(curvs[:-1], curvs[1:])
1108
1109         crossingparams = []
1110         for parampair, curvpair in zip(parampairs, curvpairs):
1111             begparam, endparam = parampair
1112             begcurv, endcurv = curvpair
1113             if (endcurv*dist - 1)*(begcurv*dist - 1) < 0:
1114                 # the curvature crosses the value 1/dist
1115                 # get the parmeter value by linear interpolation:
1116                 middleparam = (
1117                   (begparam * abs(begcurv*dist - 1) + endparam * abs(endcurv*dist - 1)) /
1118                   (abs(begcurv*dist - 1) + abs(endcurv*dist - 1)))
1119                 middleradius = normcurve.curveradius_pt([middleparam])[0]
1120
1121                 if middleradius is normpath.invalid:
1122                     raise InvalidParamException(middleparam)
1123
1124                 if abs(middleradius - dist) < epsilon:
1125                     # get the parmeter value by linear interpolation:
1126                     crossingparams.append(middleparam)
1127                 else:
1128                     # call recursively:
1129                     cps = self.distcrossingparameters(normcurve, epsilon, tstart=begparam, tend=endparam)
1130                     crossingparams += cps
1131
1132         return crossingparams
1133         # >>>
1134     def valid_near_rotation(self, nspitem, param, otherparam, stepsize, epsilon): # <<<
1135         p = param
1136         rot = nspitem.rotation([p])[0]
1137         # run towards otherparam searching for a valid rotation
1138         while rot is normpath.invalid:
1139             p = (1-stepsize)*p + stepsize*otherparam
1140             rot = nspitem.rotation([p])[0]
1141         # walk back to param until near enough
1142         # but do not go further if an invalid point is hit
1143         end, new = nspitem.at_pt([param, p])
1144         far = math.hypot(end[0]-new[0], end[1]-new[1])
1145         pnew = p
1146         while far > epsilon:
1147             pnew = (1-stepsize)*pnew + stepsize*param
1148             end, new = nspitem.at_pt([param, pnew])
1149             far = math.hypot(end[0]-new[0], end[1]-new[1])
1150             if nspitem.rotation([pnew])[0] is normpath.invalid:
1151                 break
1152             else:
1153                 p = pnew
1154         return p, nspitem.rotation([p])[0]
1155     # >>>
1156     def length_pt(self, path, param1, param2): # <<<
1157         point1, point2 = path.at_pt([param1, param2])
1158         return math.hypot(point1[0] - point2[0], point1[1] - point2[1])
1159     # >>>
1160
1161     def normpath_selfintersections(self, np, epsilon): # <<<
1162
1163         """return all self-intersection points of normpath np.
1164
1165         This does not include the intersections of a single normcurve with itself,
1166         but all intersections of one normpathitem with a different one in the path"""
1167
1168         n = len(np)
1169         linearparams = []
1170         parampairs = []
1171         paramsriap = {}
1172         for nsp_i in range(n):
1173             for nsp_j in range(nsp_i, n):
1174                 for nspitem_i in range(len(np[nsp_i])):
1175                     if nsp_j == nsp_i:
1176                         nspitem_j_range = range(nspitem_i+1, len(np[nsp_j]))
1177                     else:
1178                         nspitem_j_range = range(len(np[nsp_j]))
1179                     for nspitem_j in nspitem_j_range:
1180                         intsparams = np[nsp_i][nspitem_i].intersect(np[nsp_j][nspitem_j], epsilon)
1181                         if intsparams:
1182                             for intsparam_i, intsparam_j in intsparams:
1183                                 if ( (abs(intsparam_i) < epsilon and abs(1-intsparam_j) < epsilon) or
1184                                      (abs(intsparam_j) < epsilon and abs(1-intsparam_i) < epsilon) ):
1185                                      continue
1186                                 npp_i = normpath.normpathparam(np, nsp_i, float(nspitem_i)+intsparam_i)
1187                                 npp_j = normpath.normpathparam(np, nsp_j, float(nspitem_j)+intsparam_j)
1188                                 linearparams.append(npp_i)
1189                                 linearparams.append(npp_j)
1190                                 paramsriap[npp_i] = len(parampairs)
1191                                 paramsriap[npp_j] = len(parampairs)
1192                                 parampairs.append((npp_i, npp_j))
1193         linearparams.sort()
1194         return linearparams, parampairs, paramsriap
1195
1196     # >>>
1197     def can_continue(self, par_np, param1, param2): # <<<
1198         dist = self.dist_pt
1199
1200         rot1, rot2 = par_np.rotation([param1, param2])
1201         if rot1 is normpath.invalid or rot2 is normpath.invalid:
1202             return 0
1203         curv1, curv2 = par_np.curvature_pt([param1, param2])
1204         tang2 = rot2.apply_pt(1, 0)
1205         norm1 = rot1.apply_pt(0, -1)
1206         norm1 = (dist*norm1[0], dist*norm1[1])
1207
1208         # the self-intersection is valid if the tangents
1209         # point into the correct direction or, for parallel tangents,
1210         # if the curvature is such that the on-going path does not
1211         # enter the region defined by dist
1212         mult12 = norm1[0]*tang2[0] + norm1[1]*tang2[1]
1213         eps = 1.0e-6
1214         if abs(mult12) > eps:
1215             return (mult12 < 0)
1216         else:
1217             # tang1 and tang2 are parallel
1218             if curv2 is normpath.invalid or curv1 is normpath.invalid:
1219                 return 0
1220             if dist > 0:
1221                 return (curv2 <= curv1)
1222             else:
1223                 return (curv2 >= curv1)
1224     # >>>
1225     def rebuild_intersected_normpath(self, par_np, orig_np, epsilon): # <<<
1226
1227         dist = self.dist_pt
1228
1229         # calculate the self-intersections of the par_np
1230         selfintparams, selfintpairs, selfintsriap = self.normpath_selfintersections(par_np, epsilon)
1231         # calculate the intersections of the par_np with the original path
1232         origintparams = par_np.intersect(orig_np)[0]
1233
1234         # visualize the intersection points: # <<<
1235         if self.debug is not None:
1236             for param1, param2 in selfintpairs:
1237                 point1, point2 = par_np.at([param1, param2])
1238                 self.debug.fill(path.circle(point1[0], point1[1], 0.05), [color.rgb.red])
1239                 self.debug.fill(path.circle(point2[0], point2[1], 0.03), [color.rgb.black])
1240             for param in origintparams:
1241                 point = par_np.at([param])[0]
1242                 self.debug.fill(path.circle(point[0], point[1], 0.05), [color.rgb.green])
1243         # >>>
1244
1245         result = normpath.normpath()
1246         if not selfintparams:
1247             if origintparams:
1248                 return result
1249             else:
1250                 return par_np
1251
1252         beginparams = []
1253         endparams = []
1254         for i in range(len(par_np)):
1255             beginparams.append(normpath.normpathparam(par_np, i, 0))
1256             endparams.append(normpath.normpathparam(par_np, i, len(par_np[i])))
1257
1258         allparams = selfintparams + origintparams + beginparams + endparams
1259         allparams.sort()
1260         allparamindices = {}
1261         for i, param in enumerate(allparams):
1262             allparamindices[param] = i
1263
1264         done = {}
1265         for param in allparams:
1266             done[param] = 0
1267
1268         def otherparam(p): # <<<
1269             pair = selfintpairs[selfintsriap[p]]
1270             if (p is pair[0]):
1271                 return pair[1]
1272             else:
1273                 return pair[0]
1274         # >>>
1275         def trial_parampairs(startp): # <<<
1276             tried = {}
1277             for param in allparams:
1278                 tried[param] = done[param]
1279
1280             lastp = startp
1281             currentp = allparams[allparamindices[startp] + 1]
1282             result = []
1283
1284             while 1:
1285                 if currentp is startp:
1286                     result.append((lastp, currentp))
1287                     return result
1288                 if currentp in selfintparams and otherparam(currentp) is startp:
1289                     result.append((lastp, currentp))
1290                     return result
1291                 if currentp in endparams:
1292                     result.append((lastp, currentp))
1293                     return result
1294                 if tried[currentp]:
1295                     return []
1296                 if currentp in origintparams:
1297                     return []
1298                 # follow the crossings on valid startpairs until
1299                 # the normsubpath is closed or the end is reached
1300                 if (currentp in selfintparams and
1301                     self.can_continue(par_np, currentp, otherparam(currentp))):
1302                     # go to the next pair on the curve, seen from currentpair[1]
1303                     result.append((lastp, currentp))
1304                     lastp = otherparam(currentp)
1305                     tried[currentp] = 1
1306                     tried[otherparam(currentp)] = 1
1307                     currentp = allparams[allparamindices[otherparam(currentp)] + 1]
1308                 else:
1309                     # go to the next pair on the curve, seen from currentpair[0]
1310                     tried[currentp] = 1
1311                     tried[otherparam(currentp)] = 1
1312                     currentp = allparams[allparamindices[currentp] + 1]
1313             assert 0
1314         # >>>
1315
1316         # first the paths that start at the beginning of a subnormpath:
1317         for startp in beginparams + selfintparams:
1318             if done[startp]:
1319                 continue
1320
1321             parampairs = trial_parampairs(startp)
1322             if not parampairs:
1323                 continue
1324
1325             # collect all the pieces between parampairs
1326             add_nsp = normpath.normsubpath(epsilon=epsilon)
1327             for begin, end in parampairs:
1328                 # check that trial_parampairs works correctly
1329                 assert begin is not end
1330                 # we do not cross the border of a normsubpath here
1331                 assert begin.normsubpathindex is end.normsubpathindex
1332                 for item in par_np[begin.normsubpathindex].segments(
1333                     [begin.normsubpathparam, end.normsubpathparam])[0].normsubpathitems:
1334                     # TODO: this should be obsolete with an improved intersection algorithm
1335                     #       guaranteeing epsilon
1336                     if add_nsp.normsubpathitems:
1337                         item = item.modifiedbegin_pt(*(add_nsp.atend_pt()))
1338                     add_nsp.append(item)
1339
1340                 if begin in selfintparams:
1341                     done[begin] = 1
1342                     #done[otherparam(begin)] = 1
1343                 if end in selfintparams:
1344                     done[end] = 1
1345                     #done[otherparam(end)] = 1
1346
1347             # eventually close the path
1348             if add_nsp and (parampairs[0][0] is parampairs[-1][-1] or
1349                 (parampairs[0][0] in selfintparams and otherparam(parampairs[0][0]) is parampairs[-1][-1])):
1350                 add_nsp.normsubpathitems[-1] = add_nsp.normsubpathitems[-1].modifiedend_pt(*add_nsp.atbegin_pt())
1351                 add_nsp.close()
1352
1353             result.extend([add_nsp])
1354
1355         return result
1356
1357     # >>>
1358
1359 # >>>
1360
1361 parallel.clear = attr.clearclass(parallel)
1362
1363 # vim:foldmethod=marker:foldmarker=<<<,>>>