pyx/deformer.py

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