pyx/deformer.py

   1 # -*- encoding: utf-8 -*-
   2 #
   3 #
   4 # Copyright (C) 2003-2013 Michael Schindler <m-schindler@users.sourceforge.net>
   5 # Copyright (C) 2003-2005 André Wobst <wobsta@users.sourceforge.net>
   6 #
   7 # This file is part of PyX (http://pyx.sourceforge.net/).
   8 #
   9 # PyX is free software; you can redistribute it and/or modify
  10 # it under the terms of the GNU General Public License as published by
  11 # the Free Software Foundation; either version 2 of the License, or
  12 # (at your option) any later version.
  13 #
  14 # PyX is distributed in the hope that it will be useful,
  15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  17 # GNU General Public License for more details.
  18 #
  19 # You should have received a copy of the GNU General Public License
  20 # along with PyX; if not, write to the Free Software
  21 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  22
  23 import functools, logging, math
  24 from . import attr, baseclasses, mathutils, path, normpath, unit, color
  25
  26 normpath.invalid = 175e175 # Just a very crude workaround to get the code running again. normpath.invalid does not exist anymore.
  27
  28 logger = logging.getLogger("pyx")
  29
  30 # specific exception for an invalid parameterization point
  31 # used in parallel
  32 class InvalidParamException(Exception):
  33
  34     def __init__(self, param):
  35         self.normsubpathitemparam = param
  36
  37 # error raised in parallel if we are trying to get badly defined intersections
  38 class IntersectionError(Exception): pass
  39
  40 # None has a meaning in linesmoothed
  41 class _marker: pass
  42
  43 class inf_curvature: pass
  44
  45 def curvescontrols_from_endlines_pt(B, tangent1, tangent2, r1, r2, softness): # <<<
  46     # calculates the parameters for two bezier curves connecting two lines (curvature=0)
  47     # starting at B - r1*tangent1
  48     # ending at   B + r2*tangent2
  49     #
  50     # Takes the corner B
  51     # and two tangent vectors heading to and from B
  52     # and two radii r1 and r2:
  53     # All arguments must be in Points
  54     # Returns the seven control points of the two bezier curves:
  55     #  - start d1
  56     #  - control points g1 and f1
  57     #  - midpoint e
  58     #  - control points f2 and g2
  59     #  - endpoint d2
  60
  61     # make direction vectors d1: from B to A
  62     #                        d2: from B to C
  63     d1 = -tangent1[0] / math.hypot(*tangent1), -tangent1[1] / math.hypot(*tangent1)
  64     d2 =  tangent2[0] / math.hypot(*tangent2),  tangent2[1] / math.hypot(*tangent2)
  65
  66     # 0.3192 has turned out to be the maximum softness available
  67     # for straight lines ;-)
  68     f = 0.3192 * softness
  69     g = (15.0 * f + math.sqrt(-15.0*f*f + 24.0*f))/12.0
  70
  71     # make the control points of the two bezier curves
  72     f1 = B[0] + f * r1 * d1[0], B[1] + f * r1 * d1[1]
  73     f2 = B[0] + f * r2 * d2[0], B[1] + f * r2 * d2[1]
  74     g1 = B[0] + g * r1 * d1[0], B[1] + g * r1 * d1[1]
  75     g2 = B[0] + g * r2 * d2[0], B[1] + g * r2 * d2[1]
  76     d1 = B[0] +     r1 * d1[0], B[1] +     r1 * d1[1]
  77     d2 = B[0] +     r2 * d2[0], B[1] +     r2 * d2[1]
  78     e  = 0.5 * (f1[0] + f2[0]), 0.5 * (f1[1] + f2[1])
  79
  80     return (d1, g1, f1, e, f2, g2, d2)
  81 # >>>
  82
  83 def controldists_from_endgeometry_pt(A, B, tangA, tangB, curvA, curvB, allownegative=False, curv_epsilon=1.0e-8): # <<<
  84
  85     """For a curve with given tangents and curvatures at the endpoints this gives the distances between the controlpoints
  86
  87     This helper routine returns a list of two distances between the endpoints and the
  88     corresponding control points of a (cubic) bezier curve that has
  89     prescribed tangents tangentA, tangentB and curvatures curvA, curvB at the
  90     end points.
  91
  92     Note: The returned distances are not always positive.
  93           But only positive values are geometrically correct, so please check!
  94           The outcome is sorted so that the first entry is expected to be the
  95           most reasonable one
  96     """
  97     debug = 0
  98
  99     def test_divisions(T, D, E, AB, curvA, curvB, debug):# <<<
 100         small = AB * 1.0e-4 # at infinite curvature, avoid setting control points exactly on the startpoint
 101                             # TODO: is this consistent with the avoiding of curv=inf in normpath?
 102         arbitrary = AB * 0.33 # at zero curvature, we know nothing about a or b
 103
 104         def is_zero(x):
 105             return abs(x) < curv_epsilon
 106             # the following gave different results for forward/reverse paths
 107             # in test/functional/test_deformer parallel G
 108             #try:
 109             #    1.0 / x
 110             #except ZeroDivisionError:
 111             #    return True
 112             #return False
 113
 114
 115         if is_zero(T):
 116
 117             if curvA is inf_curvature:
 118                a = small
 119                if curvB is inf_curvature:
 120                    b = small
 121                elif is_zero(curvB):
 122                     assert abs(E) < 1.0e-10
 123                     b = arbitrary
 124                else:
 125                     b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
 126             elif is_zero(curvA):
 127                 assert abs(D) < 1.0e-10
 128                 a = arbitrary
 129                 if curvB is inf_curvature:
 130                     b = small
 131                 elif is_zero(curvB):
 132                     assert abs(E) < 1.0e-10
 133                     b = arbitrary
 134                 else:
 135                     b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
 136             else:
 137                 a = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA)
 138                 if curvB is inf_curvature:
 139                     b = small
 140                 elif is_zero(curvB):
 141                     assert abs(E) < 1.0e-10
 142                     b = arbitrary
 143                 else:
 144                     b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
 145
 146         else:
 147             if curvA is inf_curvature:
 148                a = small
 149                if curvB is inf_curvature:
 150                    b = small
 151                elif is_zero(curvB):
 152                    b = arbitrary
 153                else:
 154                    b1 = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
 155                    b2 = D / T
 156                    if abs(b1) < abs(b2):
 157                        b = b1
 158                    else:
 159                        b = b2
 160             elif curvB is inf_curvature:
 161                b = small
 162                if is_zero(curvA):
 163                    a = arbitrary
 164                else:
 165                    a1 = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA)
 166                    a2 = E / T
 167                    if abs(a1) < abs(a2):
 168                        a = a1
 169                    else:
 170                        a = a2
 171             elif is_zero(curvA):
 172                 b = D / T
 173                 a = (E - 1.5*curvB*b*abs(b)) / T
 174             elif is_zero(curvB):
 175                 a = E / T
 176                 b = (D - 1.5*curvA*a*abs(a)) / T
 177             else:
 178                 return []
 179
 180         if debug:
 181             print("fallback with exact zero value")
 182         return [(a, b)]
 183     # >>>
 184     def fallback_smallT(T, D, E, AB, curvA, curvB, threshold, debug):# <<<
 185         a = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA)
 186         b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
 187         q1 = min(abs(1.5*a*a*curvA), abs(D))
 188         q2 = min(abs(1.5*b*b*curvB), abs(E))
 189         if (a >= 0 and b >= 0 and
 190             abs(b*T) < threshold * q1 and abs(1.5*a*abs(a)*curvA - D) < threshold * q1 and
 191             abs(a*T) < threshold * q2 and abs(1.5*b*abs(b)*curvB - E) < threshold * q2):
 192             if debug:
 193                 print("fallback with T approx 0")
 194             return [(a, b)]
 195         return []
 196     # >>>
 197     def fallback_smallcurv(T, D, E, AB, curvA, curvB, threshold, debug):# <<<
 198         result = []
 199
 200         # is curvB approx zero?
 201         a = E / T
 202         b = (D - 1.5*curvA*a*abs(a)) / T
 203         if (a >= 0 and b >= 0 and
 204             abs(1.5*b*b*curvB) < threshold * min(abs(a*T), abs(E)) and
 205             abs(a*T - E) < threshold * min(abs(a*T), abs(E))):
 206             if debug:
 207                 print("fallback with curvB approx 0")
 208             result.append((a, b))
 209
 210         # is curvA approx zero?
 211         b = D / T
 212         a = (E - 1.5*curvB*b*abs(b)) / T
 213         if (a >= 0 and b >= 0 and
 214             abs(1.5*a*a*curvA) < threshold * min(abs(b*T), abs(D)) and
 215             abs(b*T - D) < threshold * min(abs(b*T), abs(D))):
 216             if debug:
 217                 print("fallback with curvA approx 0")
 218             result.append((a, b))
 219
 220         return result
 221     # >>>
 222     def findnearest(x, ys): # <<<
 223         I = 0
 224         Y = ys[I]
 225         mindist = abs(x - Y)
 226
 227         # find the value in ys which is nearest to x
 228         for i, y in enumerate(ys[1:]):
 229             dist = abs(x - y)
 230             if dist < mindist:
 231                 I, Y, mindist = i, y, dist
 232
 233         return I, Y
 234     # >>>
 235
 236     # some shortcuts
 237     T = tangA[0] * tangB[1] - tangA[1] * tangB[0]
 238     D = tangA[0] * (B[1]-A[1]) - tangA[1] * (B[0]-A[0])
 239     E = tangB[0] * (A[1]-B[1]) - tangB[1] * (A[0]-B[0])
 240     AB = math.hypot(A[0] - B[0], A[1] - B[1])
 241
 242     # try if one of the prefactors is exactly zero
 243     testsols = test_divisions(T, D, E, AB, curvA, curvB, debug)
 244     if testsols:
 245         return testsols
 246
 247     # The general case:
 248     # we try to find all the zeros of the decoupled 4th order problem
 249     # for the combined problem:
 250     # The control points of a cubic Bezier curve are given by a, b:
 251     #     A, A + a*tangA, B - b*tangB, B
 252     # for the derivation see /design/beziers.tex
 253     #     0 = 1.5 a |a| curvA + b * T - D
 254     #     0 = 1.5 b |b| curvB + a * T - E
 255     # because of the absolute values we get several possibilities for the signs
 256     # in the equation. We test all signs, also the invalid ones!
 257     if allownegative:
 258         signs = [(+1, +1), (-1, +1), (+1, -1), (-1, -1)]
 259     else:
 260         signs = [(+1, +1)]
 261
 262     candidates_a = []
 263     candidates_b = []
 264     for sign_a, sign_b in signs:
 265         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)
 266         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)
 267         candidates_a += [root for root in mathutils.realpolyroots(*coeffs_a) if sign_a*root >= 0]
 268         candidates_b += [root for root in mathutils.realpolyroots(*coeffs_b) if sign_b*root >= 0]
 269     solutions = []
 270     if candidates_a and candidates_b:
 271         for a in candidates_a:
 272             i, b = findnearest((D - 1.5*curvA*a*abs(a))/T, candidates_b)
 273             solutions.append((a, b))
 274
 275     # try if there is an approximate solution
 276     for thr in [1.0e-2, 1.0e-1]:
 277         if not solutions:
 278             solutions = fallback_smallT(T, D, E, AB, curvA, curvB, thr, debug)
 279         if not solutions:
 280             solutions = fallback_smallcurv(T, D, E, AB, curvA, curvB, thr, debug)
 281
 282     # sort the solutions: the more reasonable values at the beginning
 283     def mycmp(x,y): # <<<
 284         # first the pairs that are purely positive, then all the pairs with some negative signs
 285         # inside the two sets: sort by magnitude
 286         sx = (x[0] > 0 and x[1] > 0)
 287         sy = (y[0] > 0 and y[1] > 0)
 288
 289         # experimental stuff:
 290         # what criterion should be used for sorting ?
 291         #
 292         #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)
 293         #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)
 294         # # For each equation, a value like
 295         # #   abs(1.5*curvA*y[0]*abs(y[0]) + y[1]*T - D) / abs(curvA*(D - y[1]*T))
 296         # # indicates how good the solution is. In order to avoid the division,
 297         # # we here multiply with all four denominators:
 298         # 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)) ),
 299         #            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)) ))
 300         # 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)) ),
 301         #            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)) ))
 302         #errx = (abs(curvA*x[0]) - 1.0)**2 + (abs(curvB*x[1]) - 1.0)**2
 303         #erry = (abs(curvA*y[0]) - 1.0)**2 + (abs(curvB*y[1]) - 1.0)**2
 304
 305         errx = x[0]**2 + x[1]**2
 306         erry = y[0]**2 + y[1]**2
 307
 308         if sx == 1 and sy == 1:
 309             # try to use longer solutions if there are any crossings in the control-arms
 310             # the following combination yielded fewest sorting errors in test_bezier.py
 311             t, s = intersection(A, B, tangA, tangB)
 312             t, s = abs(t), abs(s)
 313             if (t > 0 and t < x[0] and s > 0 and s < x[1]):
 314                 if (t > 0 and t < y[0] and s > 0 and s < y[1]):
 315                     # use the shorter one
 316                     return (errx > erry) - (errx < erry)
 317                 else:
 318                     # use the longer one
 319                     return -1
 320             else:
 321                 if (t > 0 and t < y[0] and s > 0 and s < y[1]):
 322                     # use the longer one
 323                     return 1
 324                 else:
 325                     # use the shorter one
 326                     return (errx > erry) - (errx < erry)
 327             #return cmp(x[0]**2 + x[1]**2, y[0]**2 + y[1]**2)
 328         else:
 329             return (sy > sx) - (sy < sx)
 330     # >>>
 331     solutions.sort(key=functools.cmp_to_key(mycmp))
 332
 333     return solutions
 334 # >>>
 335
 336 def normcurve_from_endgeometry_pt(A, B, tangA, tangB, curvA, curvB): # <<<
 337     a, b = controldists_from_endgeometry_pt(A, B, tangA, tangB, curvA, curvB)[0]
 338     return normpath.normcurve_pt(A[0], A[1],
 339         A[0] + a * tangA[0], A[1] + a * tangA[1],
 340         B[0] - b * tangB[0], B[1] - b * tangB[1], B[0], B[1])
 341     # >>>
 342
 343 def intersection(A, D, tangA, tangD): # <<<
 344
 345     """returns the intersection parameters of two evens
 346
 347     they are defined by:
 348       x(t) = A + t * tangA
 349       x(s) = D + s * tangD
 350     """
 351     det = -tangA[0] * tangD[1] + tangA[1] * tangD[0]
 352     try:
 353         1.0 / det
 354     except ArithmeticError:
 355         return None, None
 356
 357     DA = D[0] - A[0], D[1] - A[1]
 358
 359     t = (-tangD[1]*DA[0] + tangD[0]*DA[1]) / det
 360     s = (-tangA[1]*DA[0] + tangA[0]*DA[1]) / det
 361
 362     return t, s
 363 # >>>
 364
 365 class cycloid(baseclasses.deformer): # <<<
 366     """Wraps a cycloid around a path.
 367
 368     The outcome looks like a spring with the originalpath as the axis.
 369     radius: radius of the cycloid
 370     halfloops:  number of halfloops
 371     skipfirst/skiplast: undeformed end lines of the original path
 372     curvesperhloop:
 373     sign: start left (1) or right (-1) with the first halfloop
 374     turnangle: angle of perspective on a (3D) spring
 375                turnangle=0 will produce a sinus-like cycloid,
 376                turnangle=90 will procude a row of connected circles
 377
 378     """
 379
 380     def __init__(self, radius=0.5*unit.t_cm, halfloops=10,
 381     skipfirst=1*unit.t_cm, skiplast=1*unit.t_cm, curvesperhloop=3, sign=1, turnangle=45):
 382         self.skipfirst = skipfirst
 383         self.skiplast = skiplast
 384         self.radius = radius
 385         self.halfloops = halfloops
 386         self.curvesperhloop = curvesperhloop
 387         self.sign = sign
 388         self.turnangle = turnangle
 389
 390     def __call__(self, radius=None, halfloops=None,
 391     skipfirst=None, skiplast=None, curvesperhloop=None, sign=None, turnangle=None):
 392         if radius is None:
 393             radius = self.radius
 394         if halfloops is None:
 395             halfloops = self.halfloops
 396         if skipfirst is None:
 397             skipfirst = self.skipfirst
 398         if skiplast is None:
 399             skiplast = self.skiplast
 400         if curvesperhloop is None:
 401             curvesperhloop = self.curvesperhloop
 402         if sign is None:
 403             sign = self.sign
 404         if turnangle is None:
 405             turnangle = self.turnangle
 406
 407         return cycloid(radius=radius, halfloops=halfloops, skipfirst=skipfirst, skiplast=skiplast,
 408                        curvesperhloop=curvesperhloop, sign=sign, turnangle=turnangle)
 409
 410     def deform(self, basepath):
 411         resultnormsubpaths = [self.deformsubpath(nsp) for nsp in basepath.normpath().normsubpaths]
 412         return normpath.normpath(resultnormsubpaths)
 413
 414     def deformsubpath(self, normsubpath):
 415
 416         skipfirst = abs(unit.topt(self.skipfirst))
 417         skiplast = abs(unit.topt(self.skiplast))
 418         radius = abs(unit.topt(self.radius))
 419         turnangle = math.radians(self.turnangle)
 420         sign = mathutils.sign(self.sign)
 421
 422         cosTurn = math.cos(turnangle)
 423         sinTurn = math.sin(turnangle)
 424
 425         # make list of the lengths and parameters at points on normsubpath
 426         # where we will add cycloid-points
 427         totlength = normsubpath.arclen_pt()
 428         if totlength <= skipfirst + skiplast + 2*radius*sinTurn:
 429             logger.warning("normsubpath is too short for deformation with cycloid -- skipping...")
 430             return normsubpath
 431
 432         # parameterization is in rotation-angle around the basepath
 433         # differences in length, angle ... between two basepoints
 434         # and between basepoints and controlpoints
 435         Dphi = math.pi / self.curvesperhloop
 436         phis = [i * Dphi for i in range(self.halfloops * self.curvesperhloop + 1)]
 437         DzDphi = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * math.pi * cosTurn)
 438         # Dz = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * self.curvesperhloop * cosTurn)
 439         # zs = [i * Dz for i in range(self.halfloops * self.curvesperhloop + 1)]
 440         # from path._arctobcurve:
 441         # optimal relative distance along tangent for second and third control point
 442         L = 4 * radius * (1 - math.cos(Dphi/2)) / (3 * math.sin(Dphi/2))
 443
 444         # Now the transformation of z into the turned coordinate system
 445         Zs = [ skipfirst + radius*sinTurn # here the coordinate z starts
 446              - sinTurn*radius*math.cos(phi) + cosTurn*DzDphi*phi # the transformed z-coordinate
 447              for phi in phis]
 448         params = normsubpath._arclentoparam_pt(Zs)[0]
 449
 450         # get the positions of the splitpoints in the cycloid
 451         points = []
 452         for phi, param in zip(phis, params):
 453             # the cycloid is a circle that is stretched along the normsubpath
 454             # here are the points of that circle
 455             basetrafo = normsubpath.trafo([param])[0]
 456
 457             # The point on the cycloid, in the basepath's local coordinate system
 458             baseZ, baseY = 0, radius*math.sin(phi)
 459
 460             # The tangent there, also in local coords
 461             tangentX = -cosTurn*radius*math.sin(phi) + sinTurn*DzDphi
 462             tangentY = radius*math.cos(phi)
 463             tangentZ = sinTurn*radius*math.sin(phi) + DzDphi*cosTurn
 464             norm = math.sqrt(tangentX*tangentX + tangentY*tangentY + tangentZ*tangentZ)
 465             tangentY, tangentZ = tangentY/norm, tangentZ/norm
 466
 467             # Respect the curvature of the basepath for the cycloid's curvature
 468             # XXX this is only a heuristic, not a "true" expression for
 469             #     the curvature in curved coordinate systems
 470             try:
 471                 pathradius = 1/normsubpath.curvature_pt([param])[0]
 472             except ArithmeticError:
 473                 factor = 1
 474             else:
 475                 factor = (pathradius - baseY) / pathradius
 476                 factor = abs(factor)
 477             l = L * factor
 478
 479             # The control points prior and after the point on the cycloid
 480             preeZ, preeY = baseZ - l * tangentZ, baseY - l * tangentY
 481             postZ, postY = baseZ + l * tangentZ, baseY + l * tangentY
 482
 483             # Now put everything at the proper place
 484             points.append(basetrafo.apply_pt(preeZ, sign * preeY) +
 485                           basetrafo.apply_pt(baseZ, sign * baseY) +
 486                           basetrafo.apply_pt(postZ, sign * postY))
 487
 488         if len(points) <= 1:
 489             logger.warning("normsubpath is too short for deformation with cycloid -- skipping...")
 490             return normsubpath
 491
 492         # Build the path from the pointlist
 493         # containing (control x 2,  base x 2, control x 2)
 494         if skipfirst > normsubpath.epsilon:
 495             normsubpathitems = normsubpath.segments([0, params[0]])[0]
 496             normsubpathitems.append(normpath.normcurve_pt(*(points[0][2:6] + points[1][0:4])))
 497         else:
 498             normsubpathitems = [normpath.normcurve_pt(*(points[0][2:6] + points[1][0:4]))]
 499         for i in range(1, len(points)-1):
 500             normsubpathitems.append(normpath.normcurve_pt(*(points[i][2:6] + points[i+1][0:4])))
 501         if skiplast > normsubpath.epsilon:
 502             for nsp in normsubpath.segments([params[-1], len(normsubpath)]):
 503                 normsubpathitems.extend(nsp.normsubpathitems)
 504
 505         # That's it
 506         return normpath.normsubpath(normsubpathitems, epsilon=normsubpath.epsilon)
 507 # >>>
 508
 509 cycloid.clear = attr.clearclass(cycloid)
 510
 511 class cornersmoothed(baseclasses.deformer): # <<<
 512
 513     """Bends corners in a normpath.
 514
 515     This decorator replaces corners in a normpath with bezier curves. There are two cases:
 516     - If the corner lies between two lines, _two_ bezier curves will be used
 517       that are highly optimized to look good (their curvature is to be zero at the ends
 518       and has to have zero derivative in the middle).
 519       Additionally, it can controlled by the softness-parameter.
 520     - If the corner lies between curves then _one_ bezier is used that is (except in some
 521       special cases) uniquely determined by the tangents and curvatures at its end-points.
 522       In some cases it is necessary to use only the absolute value of the curvature to avoid a
 523       cusp-shaped connection of the new bezier to the old path. In this case the use of
 524       "obeycurv=0" allows the sign of the curvature to switch.
 525     - The radius argument gives the arclength-distance of the corner to the points where the
 526       old path is cut and the beziers are inserted.
 527     - Path elements that are too short (shorter than the radius) are skipped
 528     """
 529
 530     def __init__(self, radius, softness=1, obeycurv=0, relskipthres=0.01):
 531         self.radius = radius
 532         self.softness = softness
 533         self.obeycurv = obeycurv
 534         self.relskipthres = relskipthres
 535
 536     def __call__(self, radius=None, softness=None, obeycurv=None, relskipthres=None):
 537         if radius is None:
 538             radius = self.radius
 539         if softness is None:
 540             softness = self.softness
 541         if obeycurv is None:
 542             obeycurv = self.obeycurv
 543         if relskipthres is None:
 544             relskipthres = self.relskipthres
 545         return cornersmoothed(radius=radius, softness=softness, obeycurv=obeycurv, relskipthres=relskipthres)
 546
 547     def deform(self, basepath):
 548         return normpath.normpath([self.deformsubpath(normsubpath)
 549                               for normsubpath in basepath.normpath().normsubpaths])
 550
 551     def deformsubpath(self, normsubpath):
 552         radius_pt = unit.topt(self.radius)
 553         epsilon = normsubpath.epsilon
 554
 555         # remove too short normsubpath items (shorter than self.relskipthres*radius_pt or epsilon)
 556         pertinentepsilon = max(epsilon, self.relskipthres*radius_pt)
 557         pertinentnormsubpath = normpath.normsubpath(normsubpath.normsubpathitems,
 558                                                 epsilon=pertinentepsilon)
 559         pertinentnormsubpath.flushskippedline()
 560         pertinentnormsubpathitems = pertinentnormsubpath.normsubpathitems
 561
 562         # calculate the splitting parameters for the pertinentnormsubpathitems
 563         arclens_pt = []
 564         params = []
 565         for pertinentnormsubpathitem in pertinentnormsubpathitems:
 566             arclen_pt = pertinentnormsubpathitem.arclen_pt(epsilon)
 567             arclens_pt.append(arclen_pt)
 568             l1_pt = min(radius_pt, 0.5*arclen_pt)
 569             l2_pt = max(0.5*arclen_pt, arclen_pt - radius_pt)
 570             params.append(pertinentnormsubpathitem.arclentoparam_pt([l1_pt, l2_pt], epsilon))
 571
 572         # handle the first and last pertinentnormsubpathitems for a non-closed normsubpath
 573         if not normsubpath.closed:
 574             l1_pt = 0
 575             l2_pt = max(0, arclens_pt[0] - radius_pt)
 576             params[0] = pertinentnormsubpathitems[0].arclentoparam_pt([l1_pt, l2_pt], epsilon)
 577             l1_pt = min(radius_pt, arclens_pt[-1])
 578             l2_pt = arclens_pt[-1]
 579             params[-1] = pertinentnormsubpathitems[-1].arclentoparam_pt([l1_pt, l2_pt], epsilon)
 580
 581         newnormsubpath = normpath.normsubpath(epsilon=normsubpath.epsilon)
 582         for i in range(len(pertinentnormsubpathitems)):
 583             this = i
 584             next = (i+1) % len(pertinentnormsubpathitems)
 585             thisparams = params[this]
 586             nextparams = params[next]
 587             thisnormsubpathitem = pertinentnormsubpathitems[this]
 588             nextnormsubpathitem = pertinentnormsubpathitems[next]
 589             thisarclen_pt = arclens_pt[this]
 590             nextarclen_pt = arclens_pt[next]
 591
 592             # insert the middle segment
 593             newnormsubpath.append(thisnormsubpathitem.segments(thisparams)[0])
 594
 595             # insert replacement curves for the corners
 596             if next or normsubpath.closed:
 597
 598                 t1 = thisnormsubpathitem.rotation([thisparams[1]])[0].apply_pt(1, 0)
 599                 t2 = nextnormsubpathitem.rotation([nextparams[0]])[0].apply_pt(1, 0)
 600                 # TODO: normpath.invalid
 601
 602                 if (isinstance(thisnormsubpathitem, normpath.normline_pt) and
 603                     isinstance(nextnormsubpathitem, normpath.normline_pt)):
 604
 605                     # case of two lines -> replace by two curves
 606                     d1, g1, f1, e, f2, g2, d2 = curvescontrols_from_endlines_pt(
 607                         thisnormsubpathitem.atend_pt(), t1, t2,
 608                         thisarclen_pt*(1-thisparams[1]), nextarclen_pt*(nextparams[0]), softness=self.softness)
 609
 610                     p1 = thisnormsubpathitem.at_pt([thisparams[1]])[0]
 611                     p2 = nextnormsubpathitem.at_pt([nextparams[0]])[0]
 612
 613                     newnormsubpath.append(normpath.normcurve_pt(*(d1 + g1 + f1 + e)))
 614                     newnormsubpath.append(normpath.normcurve_pt(*(e + f2 + g2 + d2)))
 615
 616                 else:
 617
 618                     # generic case -> replace by a single curve with prescribed tangents and curvatures
 619                     p1 = thisnormsubpathitem.at_pt([thisparams[1]])[0]
 620                     p2 = nextnormsubpathitem.at_pt([nextparams[0]])[0]
 621                     c1 = thisnormsubpathitem.curvature_pt([thisparams[1]])[0]
 622                     c2 = nextnormsubpathitem.curvature_pt([nextparams[0]])[0]
 623                     # TODO: normpath.invalid
 624
 625                     # TODO: more intelligent fallbacks:
 626                     #   circle -> line
 627                     #   circle -> circle
 628
 629                     if not self.obeycurv:
 630                         # do not obey the sign of the curvature but
 631                         # make the sign such that the curve smoothly passes to the next point
 632                         # this results in a discontinuous curvature
 633                         # (but the absolute value is still continuous)
 634                         s1 = +mathutils.sign(t1[0] * (p2[1]-p1[1]) - t1[1] * (p2[0]-p1[0]))
 635                         s2 = -mathutils.sign(t2[0] * (p2[1]-p1[1]) - t2[1] * (p2[0]-p1[0]))
 636                         c1 = s1 * abs(c1)
 637                         c2 = s2 * abs(c2)
 638
 639                     # get the length of the control "arms"
 640                     controldists = controldists_from_endgeometry_pt(p1, p2, t1, t2, c1, c2)
 641
 642                     if controldists and (controldists[0][0] >= 0 and controldists[0][1] >= 0):
 643                         # use the first entry in the controldists
 644                         # this should be the "smallest" pair
 645                         a, d = controldists[0]
 646                         # avoid curves with invalid parameterization
 647                         a = max(a, epsilon)
 648                         d = max(d, epsilon)
 649
 650                         # avoid overshooting at the corners:
 651                         # this changes not only the sign of the curvature
 652                         # but also the magnitude
 653                         if not self.obeycurv:
 654                             t, s = intersection(p1, p2, t1, t2)
 655                             if (t is not None and s is not None and
 656                                 t > 0 and s < 0):
 657                                 a = min(a, abs(t))
 658                                 d = min(d, abs(s))
 659
 660                     else:
 661                         # use a fallback
 662                         t, s = intersection(p1, p2, t1, t2)
 663                         if t is not None and s is not None:
 664                             a = 0.65 * abs(t)
 665                             d = 0.65 * abs(s)
 666                         else:
 667                             # if there is no useful result:
 668                             # take an arbitrary smoothing curve that does not obey
 669                             # the curvature constraints
 670                             dist = math.hypot(p1[0] - p2[0], p1[1] - p2[1])
 671                             a = dist / (3.0 * math.hypot(*t1))
 672                             d = dist / (3.0 * math.hypot(*t2))
 673
 674                     # calculate the two missing control points
 675                     q1 = p1[0] + a * t1[0], p1[1] + a * t1[1]
 676                     q2 = p2[0] - d * t2[0], p2[1] - d * t2[1]
 677
 678                     newnormsubpath.append(normpath.normcurve_pt(*(p1 + q1 + q2 + p2)))
 679
 680         if normsubpath.closed:
 681             newnormsubpath.close()
 682         return newnormsubpath
 683
 684 # >>>
 685
 686 cornersmoothed.clear = attr.clearclass(cornersmoothed)
 687 smoothed = cornersmoothed
 688 smoothed.clear = attr.clearclass(smoothed)
 689
 690 class mynormpathparam(normpath.normpathparam): # <<<
 691     """In the parallel deformer we use properties such as the curvature, which
 692     are not continuous on a path (at points between normsubpathitems). We
 693     therefore require a better parameter class which really resolves the
 694     nsp-item """
 695
 696     # TODO: find reasonable values for these eps:
 697     rounding_eps = 1.0e-8
 698
 699     def __init__(self, np, normsubpathindex, normsubpathitemindex, normsubpathitemparam):
 700         normpath.normpathparam.__init__(self, np, normsubpathindex, normsubpathitemindex + normsubpathitemparam)
 701         self.normsubpathitemindex = normsubpathitemindex
 702         self.normsubpathitemparam = normsubpathitemparam
 703         self.beg_nspitem_known = False
 704         self.end_nspitem_known = False
 705
 706         # guarantee that normpath.normpathparam always gets the correct item:
 707         if int(self.normsubpathparam) != self.normsubpathitemindex:
 708             if int(self.normsubpathparam) == self.normsubpathitemindex - 1:
 709                 self.normsubpathparam = self.normsubpathitemindex + self.rounding_eps
 710                 self.beg_nspitem_known = True
 711             elif int(self.normsubpathparam) == self.normsubpathitemindex + 1:
 712                 self.normsubpathparam = (self.normsubpathitemindex + 1) - self.rounding_eps
 713                 self.end_nspitem_known = True
 714             else:
 715                 assert False
 716         assert int(self.normsubpathparam) == self.normsubpathitemindex
 717         #assert 0 <= self.normsubpathparam - self.normsubpathitemindex
 718         #assert 1 >= self.normsubpathparam - self.normsubpathitemindex
 719
 720     def __str__(self):
 721         return "npp(%d, %d, %.3f)" % (self.normsubpathindex, self.normsubpathitemindex, self.normsubpathitemparam)
 722
 723     def __eq__(self, other):
 724         if isinstance(other, mynormpathparam):
 725             assert self.normpath is other.normpath, "normpathparams have to belong to the same normpath"
 726             return (self.normsubpathindex, self.normsubpathitemindex, self.normsubpathitemparam) == (other.normsubpathindex, other.normsubpathitemindex, other.normsubpathitemparam)
 727         else:
 728             normpath.normpathparam.__eq__(self, other)
 729
 730     def __lt__(self, other):
 731         if isinstance(other, mynormpathparam):
 732             assert self.normpath is other.normpath, "normpathparams have to belong to the same normpath"
 733             return (self.normsubpathindex, self.normsubpathitemindex, self.normsubpathitemparam) < (other.normsubpathindex, other.normsubpathitemindex, other.normsubpathitemparam)
 734         else:
 735             normpath.normpathparam.__eq__(self, other)
 736
 737     def __hash__(self):
 738         return id(self)
 739
 740     def smaller_equiv(self, epsilon=None):
 741         """Returns smaller equivalent parameter, if self is between two nsp-items"""
 742         if not self.is_beg_of_nspitem(epsilon):
 743             return self
 744         nsp = self.normpath[self.normsubpathindex]
 745         nspi_index = self.normsubpathitemindex - 1
 746         nspi_param = 1
 747         if nsp.closed:
 748             nspi_index = nspi_index % len(nsp)
 749         elif nspi_index < 0:
 750             nspi_index = 0
 751             nspi_param = 0
 752         other = mynormpathparam(self.normpath, self.normsubpathindex, nspi_index, nspi_param)
 753         if self.is_equiv(other, epsilon):
 754             return other
 755         return self
 756
 757     def larger_equiv(self, epsilon=None):
 758         """Returns smaller equivalent parameter, if self is between two nsp-items"""
 759         if not self.is_end_of_nspitem(epsilon):
 760             return self
 761         nsp = self.normpath[self.normsubpathindex]
 762         nspi_index = self.normsubpathitemindex + 1
 763         nspi_param = 0
 764         if nsp.closed:
 765             nspi_index = nspi_index % len(nsp)
 766         elif nspi_index >= len(nsp):
 767             nspi_index = len(nsp) - 1
 768             nspi_param = 1
 769         other = mynormpathparam(self.normpath, self.normsubpathindex, nspi_index, nspi_param)
 770         if self.is_equiv(other, epsilon):
 771             return other
 772         return self
 773
 774     def is_equiv(self, other, epsilon=None):
 775         """Test whether the two params yield essentially the same point"""
 776         assert self.normpath is other.normpath, "normpathparams have to belong to the same normpath"
 777         if self.normsubpathindex != other.normsubpathindex:
 778             return False
 779         nsp = self.normpath[self.normsubpathindex]
 780         if epsilon is None:
 781             epsilon = nsp.epsilon
 782         A, B = self.normpath.at_pt([self, other])
 783         return math.hypot(A[0]-B[0], A[1]-B[1]) < epsilon
 784
 785     def is_beg_of_nspitem(self, epsilon=None):
 786         if self.beg_nspitem_known:
 787             return True
 788         return self.is_equiv(mynormpathparam(self.normpath, self.normsubpathindex, self.normsubpathitemindex, 0), epsilon)
 789
 790     def is_end_of_nspitem(self, epsilon=None):
 791         if self.end_nspitem_known:
 792             return True
 793         return self.is_equiv(mynormpathparam(self.normpath, self.normsubpathindex, self.normsubpathitemindex, 1), epsilon)
 794
 795     def is_beg_of_nsp(self, epsilon=None):
 796         if self.normsubpathitemindex > 0:
 797             return False
 798         return self.is_equiv(mynormpathparam(self.normpath, self.normsubpathindex, 0, 0), epsilon)
 799
 800     def is_end_of_nsp(self, epsilon=None):
 801         n = len(self.normpath[self.normsubpathindex]) - 1
 802         if self.normsubpathitemindex < n:
 803             return False
 804         return self.is_equiv(mynormpathparam(self.normpath, self.normsubpathindex, n, 1), epsilon)
 805
 806 # >>>
 807 def _length_pt(path, param1, param2): # <<<
 808     point1, point2 = path.at_pt([param1, param2])
 809     return math.hypot(point1[0] - point2[0], point1[1] - point2[1])
 810 # >>>
 811 class parallel(baseclasses.deformer): # <<<
 812
 813     """creates a parallel normpath with constant distance to the original normpath
 814
 815     A positive 'distance' results in a curve left of the original one -- and a
 816     negative 'distance' in a curve at the right. Left/right are understood in
 817     terms of the parameterization of the original curve. The construction of
 818     the paralel path is done in two steps: First, an "extended" parallel path
 819     is built. For each path element a parallel curve/line is constructed, which
 820     can be too long or too short, depending on the presence of corners. At
 821     corners, either a circular arc is drawn around the corner, or, if possible,
 822     the parallel curve is cut in order to also exhibit a corner. In a second
 823     step all self-intersection points are determined and unnecessary parts of
 824     the path are cut away.
 825
 826     distance:            the distance of the parallel normpath
 827     relerr:              distance*relerr is the maximal allowed error in the parallel distance
 828     sharpoutercorners:   make the outer corners not round but sharp.
 829                          The inner corners (corners after inflection points) will stay round
 830     dointersection:      boolean for doing the intersection step (default: 1).
 831                          Set this value to 0 if you want the whole parallel path
 832     checkdistanceparams: a list of parameter values in the interval (0,1) where the
 833                          parallel distance is checked on each normpathitem
 834     lookforcurvatures:   number of points per normpathitem where is looked for
 835                          a critical value of the curvature
 836     """
 837
 838     # TODO:
 839     # * DECIDE MEANING of arcs around corners (see case L in test/functional/test_deformer.py)
 840     # * eliminate double, triple, ... pairs
 841     # * implement self-intersection of normcurve_pt
 842     # * implement _between_paths also for normcurve_pt
 843
 844
 845     def __init__(self, distance, relerr=0.05, sharpoutercorners=False, dointersection=True,
 846                        checkdistanceparams=[0.5], lookforcurvatures=11, searchstep=0.01, debug=None):
 847         self.distance = distance
 848         self.relerr = relerr
 849         self.sharpoutercorners = sharpoutercorners
 850         self.checkdistanceparams = checkdistanceparams
 851         self.lookforcurvatures = lookforcurvatures
 852         self.dointersection = dointersection
 853         self.searchstep = searchstep
 854         self.debug = debug
 855
 856     def __call__(self, distance=None, relerr=None, sharpoutercorners=None, dointersection=None,
 857                        checkdistanceparams=None, lookforcurvatures=None, searchstep=None, debug=None):
 858         # returns a copy of the deformer with different parameters
 859         if distance is None:
 860             distance = self.distance
 861         if relerr is None:
 862             relerr = self.relerr
 863         if sharpoutercorners is None:
 864             sharpoutercorners = self.sharpoutercorners
 865         if dointersection is None:
 866             dointersection = self.dointersection
 867         if checkdistanceparams is None:
 868             checkdistanceparams = self.checkdistanceparams
 869         if lookforcurvatures is None:
 870             lookforcurvatures = self.lookforcurvatures
 871         if searchstep is None:
 872             searchstep = self.searchstep
 873         if debug is None:
 874             debug = self.debug
 875
 876         return parallel(distance=distance, relerr=relerr,
 877                         sharpoutercorners=sharpoutercorners,
 878                         dointersection=dointersection,
 879                         checkdistanceparams=checkdistanceparams,
 880                         lookforcurvatures=lookforcurvatures,
 881                         searchstep=searchstep,
 882                         debug=debug)
 883
 884     def deform(self, basepath):
 885         basepath = basepath.normpath()
 886         self.dist_pt = unit.topt(self.distance)
 887         resultnormsubpaths = []
 888         par_to_orig = {}
 889         for nsp in basepath.normsubpaths:
 890             parallel_normpath, tmp1, tmp2, par2orig = self.deformsubpath(nsp)
 891             resultnormsubpaths += parallel_normpath.normsubpaths
 892             for key in par2orig:
 893                 par_to_orig[key] = par2orig[key]
 894         result = normpath.normpath(resultnormsubpaths)
 895
 896         if self.dointersection:
 897             result = self.rebuild_intersected_normpath(result, basepath, par_to_orig)
 898
 899         return result
 900
 901     def deformsubpath(self, orig_nsp): # <<<
 902
 903         """Performs the first step of the deformation: Returns a list of
 904         normsubpaths building the parallel to the given normsubpath.
 905         Then calls the intersection routine to do the second step."""
 906         # the default case is not to join the result.
 907
 908         dist = self.dist_pt
 909         epsilon = orig_nsp.epsilon
 910         assert len(orig_nsp.normsubpathitems) != 0
 911
 912         # avoid too small dists: we would run into instabilities
 913         if abs(dist) < abs(epsilon):
 914             assert orig_nsp.normsubpathitems
 915             par_to_orig = {}
 916             for nspitem in orig_nsp:
 917                 par_to_orig[nspitem] = nspitem
 918             return normpath.normpath([orig_nsp]), None, None, par_to_orig
 919
 920         result = None
 921         par_to_orig = None
 922         join_begin = None
 923         prev_joinend = None
 924
 925         # iterate over the normsubpath in the following way:
 926         # * for each item first append the additional arc
 927         #   and then add the next parallel piece
 928         # * for the first item only add the parallel piece
 929         #   (because this is done for curr_orig_nspitem, we need to start with i=0)
 930         for i in range(len(orig_nsp.normsubpathitems)):
 931             prev_orig_nspitem = orig_nsp.normsubpathitems[i-1]
 932             curr_orig_nspitem = orig_nsp.normsubpathitems[i]
 933
 934             # get the next parallel piece for the normpath
 935             next_parallel_normpath, joinbeg, joinend, par2orig = self.deformsubpathitem(curr_orig_nspitem, epsilon)
 936             if result is None:
 937                 if join_begin is None:
 938                     join_begin = joinbeg
 939                 else:
 940                     join_begin = (join_begin and joinbeg)
 941
 942             if not (next_parallel_normpath.normsubpaths and next_parallel_normpath[0].normsubpathitems):
 943                 if prev_joinend is None:
 944                     prev_joinend = joinend
 945                 else:
 946                     prev_joinend = (prev_joinend and joinend)
 947                 continue
 948
 949             # this starts the whole normpath
 950             if result is None:
 951                 result = next_parallel_normpath
 952                 par_to_orig = {}
 953                 for key in par2orig:
 954                     par_to_orig[key] = par2orig[key]
 955                 prev_joinend = joinend
 956                 continue # there is nothing to join
 957
 958             prev_tangent, next_tangent, is_straight, is_concave = self._get_angles(prev_orig_nspitem, curr_orig_nspitem, epsilon)
 959             if not (joinbeg and prev_joinend): # split due to loo large curvature
 960                 result += next_parallel_normpath
 961             elif is_straight:
 962                 # The path is quite straight between prev and next item:
 963                 # normpath.normpath.join adds a straight line if necessary
 964                 result.join(next_parallel_normpath)
 965             else:
 966                 # The parallel path can be extended continuously.
 967                 # We must add a corner or an arc around the corner:
 968                 cornerpath = self._path_around_corner(curr_orig_nspitem.atbegin_pt(), result.atend_pt(), next_parallel_normpath.atbegin_pt(),
 969                                                       prev_tangent, next_tangent, is_concave, epsilon)
 970                 result.join(cornerpath)
 971                 assert len(cornerpath) == 1
 972                 corner = curr_orig_nspitem.atbegin_pt()
 973                 for cp_nspitem in cornerpath[0]:
 974                     par_to_orig[cp_nspitem] = corner
 975                 # append the next parallel piece to the path
 976                 result.join(next_parallel_normpath)
 977             for key in par2orig:
 978                 par_to_orig[key] = par2orig[key]
 979             prev_joinend = joinend
 980
 981         # end here if nothing has been found so far
 982         if result is None:
 983             return normpath.normpath(), False, False, {}
 984
 985         # the curve around the closing corner may still be missing
 986         if orig_nsp.closed:
 987             prev_tangent, next_tangent, is_straight, is_concave = self._get_angles(orig_nsp.normsubpathitems[-1], orig_nsp.normsubpathitems[0], epsilon)
 988             if not (joinend and join_begin): # do not close because of infinite curvature
 989                 do_close = False
 990             elif is_straight:
 991                 # The path is quite straight at end and beginning
 992                 do_close = True
 993             else:
 994                 # The parallel path can be extended continuously.
 995                 do_close = True
 996                 # We must add a corner or an arc around the corner:
 997                 cornerpath = self._path_around_corner(orig_nsp.atend_pt(), result.atend_pt(), result.atbegin_pt(),
 998                                                       prev_tangent, next_tangent, is_concave, epsilon)
 999                 result.join(cornerpath)
1000                 corner = orig_nsp.atend_pt()
1001                 assert len(cornerpath) == 1
1002                 for cp_nspitem in cornerpath[0]:
1003                     par_to_orig[cp_nspitem] = corner
1004
1005             if do_close:
1006                 if len(result) == 1:
1007                     result[0].close()
1008                 else:
1009                     # if the parallel normpath is split into several subpaths anyway,
1010                     # then use the natural beginning and ending
1011                     # closing is not possible anymore
1012                     for nspitem in result[0]:
1013                         result[-1].append(nspitem)
1014                     result.normsubpaths = result.normsubpaths[1:]
1015                 join_begin, joinend = False, False
1016
1017         return result, join_begin, joinend, par_to_orig
1018         # >>>
1019     def deformsubpathitem(self, nspitem, epsilon): # <<<
1020
1021         """Returns a parallel normpath for a single normsubpathitem
1022
1023         Analyzes the curvature of a normsubpathitem and returns a normpath with
1024         the appropriate number of normsubpaths. This must be a normpath because
1025         a normcurve can be strongly curved, such that the parallel path must
1026         contain a hole"""
1027         # the default case is to join the result. Only if there was an infinite
1028         # curvature at beginning/end, we return info not to join it.
1029
1030         dist = self.dist_pt
1031
1032         # for a simple line we return immediately
1033         if isinstance(nspitem, normpath.normline_pt):
1034             normal = nspitem.rotation([0])[0].apply_pt(0, 1)
1035             start = nspitem.atbegin_pt()
1036             end = nspitem.atend_pt()
1037             result = path.line_pt(start[0] + dist * normal[0], start[1] + dist * normal[1],
1038                                   end[0] + dist * normal[0], end[1] + dist * normal[1]).normpath(epsilon=epsilon)
1039             assert len(result) == 1 and len(result[0]) == 1
1040             return result, True, True, {result[0][0]:nspitem}
1041
1042         # for a curve we have to check if the curvatures
1043         # cross the singular value 1/dist
1044         crossings = list(self._distcrossingparameters(nspitem, epsilon))
1045         crossings.sort()
1046
1047         # depending on the number of crossings we must consider
1048         # three different cases:
1049         if crossings:
1050             # The curvature crosses the borderline 1/dist
1051             # the parallel curve contains points with infinite curvature!
1052             parallcurvs = [inf_curvature]*len(crossings)
1053
1054             result = normpath.normpath()
1055             join_begin, join_end = False, False
1056             par_to_orig = {}
1057
1058             # we need the endpoints of the nspitem
1059             if _length_pt(nspitem, crossings[0], 0) > epsilon:
1060                 crossings.insert(0, 0)
1061                 parallcurvs.insert(0, None)
1062             if _length_pt(nspitem, crossings[-1], 1) > epsilon:
1063                 crossings.append(1)
1064                 parallcurvs.append(None)
1065
1066             for i in range(len(crossings) - 1):
1067                 middleparam = 0.5*(crossings[i] + crossings[i+1])
1068                 middlecurv = nspitem.curvature_pt([middleparam])[0]
1069                 # the radius is good if
1070                 #  - middlecurv and dist have opposite signs : distance vector points "out" of the original curve
1071                 #  - middlecurv is "smaller" than 1/dist : original curve is less curved than +-1/dist
1072                 if dist*middlecurv < 0 or abs(dist*middlecurv) < 1:
1073                     if i == 0:
1074                         join_begin = True
1075                     elif i == len(crossings) - 2:
1076                         join_end = True
1077                     parallel_nsp, par2orig = self.deformnicecurve(nspitem.segments(crossings[i:i+2])[0], epsilon, curvA=parallcurvs[i], curvD=parallcurvs[i+1])
1078                     # never append empty normsubpaths
1079                     if parallel_nsp.normsubpathitems:
1080                         # infinite curvatures interrupt the path and start a new nsp
1081                         result.append(parallel_nsp)
1082                         for key in par2orig:
1083                             par_to_orig[key] = par2orig[key]
1084             if not (result.normsubpaths and result[0].normsubpathitems):
1085                 return normpath.normpath(), True, True, {}
1086             return result, join_begin, join_end, par_to_orig
1087
1088         # the curvature is either bigger or smaller than 1/dist
1089         middlecurv = nspitem.curvature_pt([0.5])[0]
1090         if dist*middlecurv < 0 or abs(dist*middlecurv) < 1:
1091             # The curve is everywhere less curved than 1/dist
1092             # We can proceed finding the parallel curve for the whole piece
1093             parallel_nsp, par2orig = self.deformnicecurve(nspitem, epsilon)
1094             # never append empty normsubpaths
1095             if parallel_nsp.normsubpathitems:
1096                 par_to_orig = {}
1097                 for key in par2orig:
1098                     par_to_orig[key] = par2orig[key]
1099                 return normpath.normpath([parallel_nsp]), True, True, par_to_orig
1100
1101         # the curve is everywhere stronger curved than 1/dist
1102         # There is nothing to be returned.
1103         return normpath.normpath(), False, False, {}
1104         # >>>
1105     def deformnicecurve(self, normcurve, epsilon, startparam=0.0, endparam=1.0, curvA=None, curvD=None): # <<<
1106         """Returns a parallel normsubpath for the normcurve.
1107
1108         This routine assumes that the normcurve is everywhere
1109         'less' curved than 1/dist. Only at the ends, the curvature
1110         can be exactly +-1/dist, which is marked by curvA and/or curvD.
1111         """
1112         dist = self.dist_pt
1113
1114         # normalized tangent directions
1115         tangA, tangD = normcurve.rotation([startparam, endparam])
1116         tangA = tangA.apply_pt(1, 0)
1117         tangD = tangD.apply_pt(1, 0)
1118
1119         # the new starting points
1120         orig_A, orig_D = normcurve.at_pt([startparam, endparam])
1121         A = orig_A[0] - dist * tangA[1], orig_A[1] + dist * tangA[0]
1122         D = orig_D[0] - dist * tangD[1], orig_D[1] + dist * tangD[0]
1123
1124         # we need to end this _before_ we will run into epsilon-problems
1125         # when creating curves we do not want to calculate the length of
1126         # or even split it for recursive calls
1127         if (math.hypot(A[0] - D[0], A[1] - D[1]) < epsilon and
1128             abs(dist)*(tangA[0]*tangD[1] - tangA[1]*tangD[0]) < epsilon):
1129             nl = normpath.normline_pt(A[0], A[1], D[0], D[1])
1130             return normpath.normsubpath([nl]), {nl:normcurve}
1131
1132         result = normpath.normsubpath(epsilon=epsilon)
1133         # is there enough space on the normals before they intersect?
1134         a, d = intersection(orig_A, orig_D, (-tangA[1], tangA[0]), (-tangD[1], tangD[0]))
1135         # a,d are the lengths to the intersection points:
1136         # for a (and equally for b) we can proceed in one of the following cases:
1137         #   a is None (means parallel normals)
1138         #   a and dist have opposite signs (and the same for b)
1139         #   a has the same sign but is bigger
1140         if ( (a is None or a*dist < 0 or abs(a) > abs(dist) + epsilon) or
1141              (d is None or d*dist < 0 or abs(d) > abs(dist) + epsilon) ):
1142             # the original path is long enough to draw a parallel piece
1143             # this is the generic case. Get the parallel curves
1144             orig_curvA, orig_curvD = normcurve.curvature_pt([startparam, endparam])
1145             if curvA is None:
1146                 curvA = orig_curvA / (1.0 - dist*orig_curvA)
1147             if curvD is None:
1148                 curvD = orig_curvD / (1.0 - dist*orig_curvD)
1149
1150             # first try to approximate the normcurve with a single item
1151             controldistpairs = controldists_from_endgeometry_pt(A, D, tangA, tangD, curvA, curvD)
1152
1153             if controldistpairs:
1154                 # TODO: is it good enough to get the first entry here?
1155                 #       from testing: this fails if there are loops in the original curve
1156                 a, d = controldistpairs[0]
1157                 if a >= 0 and d >= 0:
1158                     # we avoid to create curves with invalid parameterization
1159                     if a < epsilon and d < epsilon:
1160                         result = normpath.normsubpath([normpath.normline_pt(A[0], A[1], D[0], D[1])], epsilon=epsilon)
1161                     else:
1162                         a = max(a, epsilon)
1163                         d = max(d, epsilon)
1164                         result = normpath.normsubpath([normpath.normcurve_pt(
1165                             A[0], A[1],
1166                             A[0] + a * tangA[0], A[1] + a * tangA[1],
1167                             D[0] - d * tangD[0], D[1] - d * tangD[1],
1168                             D[0], D[1])], epsilon=epsilon)
1169
1170             # then try with two items, recursive call
1171             if ((not result.normsubpathitems) or
1172                 (self.checkdistanceparams and result.normsubpathitems
1173                  and not self._distchecked(normcurve, result, epsilon, startparam, endparam))):
1174                 # TODO: does this ever converge?
1175                 # TODO: what if this hits epsilon?
1176                 middleparam = 0.5*(startparam + endparam)
1177                 firstnsp, first_par2orig = self.deformnicecurve(normcurve, epsilon, startparam, middleparam, curvA, None)
1178                 secondnsp, second_par2orig = self.deformnicecurve(normcurve, epsilon, middleparam, endparam, None, curvD)
1179                 if not (firstnsp.normsubpathitems and secondnsp.normsubpathitems):
1180                     result = normpath.normsubpath(
1181                         [normpath.normline_pt(A[0], A[1], D[0], D[1])], epsilon=epsilon)
1182                 else:
1183                     result = firstnsp.joined(secondnsp)
1184
1185         par_to_orig = {}
1186         for key in result:
1187             par_to_orig[key] = normcurve
1188         return result, par_to_orig
1189         # >>>
1190
1191     def _path_around_corner(self, corner_pt, beg_pt, end_pt, beg_tangent, end_tangent, is_concave, epsilon): # <<<
1192         """Helper routine for parallel.deformsubpath: Draws an arc around a convex corner"""
1193         if self.sharpoutercorners and not is_concave:
1194             # straight lines:
1195             t1, t2 = intersection(beg_pt, end_pt, beg_tangent, end_tangent)
1196             B = beg_pt[0] + t1 * beg_tangent[0], beg_pt[1] + t1 * beg_tangent[1]
1197             return normpath.normpath([normpath.normsubpath([
1198                 normpath.normline_pt(beg_pt[0], beg_pt[1], B[0], B[1]),
1199                 normpath.normline_pt(B[0], B[1], end_pt[0], end_pt[1])
1200                 ])])
1201
1202         # We append an arc around the corner
1203         # these asserts fail in test case "E"
1204         #assert abs(math.hypot(beg_pt[1] - corner_pt[1], beg_pt[0] - corner_pt[0]) - abs(self.dist_pt)) < epsilon
1205         #assert abs(math.hypot(end_pt[1] - corner_pt[1], end_pt[0] - corner_pt[0]) - abs(self.dist_pt)) < epsilon
1206         angle1 = math.atan2(beg_pt[1] - corner_pt[1], beg_pt[0] - corner_pt[0])
1207         angle2 = math.atan2(end_pt[1] - corner_pt[1], end_pt[0] - corner_pt[0])
1208
1209         # depending on the direction we have to use arc or arcn
1210         sinangle = beg_tangent[0]*end_tangent[1] - beg_tangent[1]*end_tangent[0] # >0 for left-turning, <0 for right-turning
1211         if self.dist_pt > 0:
1212             arcclass = path.arcn_pt
1213         else:
1214             arcclass = path.arc_pt
1215         return path.path(arcclass(
1216           corner_pt[0], corner_pt[1], abs(self.dist_pt),
1217           math.degrees(angle1), math.degrees(angle2))).normpath(epsilon=epsilon)
1218     # >>>
1219     def _distchecked(self, orig_normcurve, parallel_normsubpath, epsilon, tstart, tend): # <<<
1220         """Helper routine for parallel.deformnicecurve: Checks the distances between orig_normcurve and parallel_normsubpath.
1221
1222         The checking is done at parameters self.checkdistanceparams of orig_normcurve."""
1223
1224         dist = self.dist_pt
1225         # do not look closer than epsilon:
1226         dist_err = mathutils.sign(dist) * max(abs(self.relerr*dist), epsilon)
1227
1228         checkdistanceparams = [tstart + (tend-tstart)*t for t in self.checkdistanceparams]
1229
1230         for param, P, rotation in zip(checkdistanceparams,
1231                                       orig_normcurve.at_pt(checkdistanceparams),
1232                                       orig_normcurve.rotation(checkdistanceparams)):
1233             normal = rotation.apply_pt(0, 1)
1234
1235             # create a short cutline for intersection only:
1236             cutline = normpath.normsubpath([normpath.normline_pt(
1237               P[0] + (dist - 2*dist_err) * normal[0], P[1] + (dist - 2*dist_err) * normal[1],
1238               P[0] + (dist + 2*dist_err) * normal[0], P[1] + (dist + 2*dist_err) * normal[1])], epsilon=epsilon)
1239
1240             cutparams = parallel_normsubpath.intersect(cutline)
1241             distances = [math.hypot(P[0] - cutpoint[0], P[1] - cutpoint[1])
1242                          for cutpoint in cutline.at_pt(cutparams[1])]
1243
1244             if (not distances) or (abs(min(distances) - abs(dist)) > abs(dist_err)):
1245                 return False
1246
1247         return True
1248     # >>>
1249     def _distcrossingparameters(self, normcurve, epsilon, tstart=0, tend=1): # <<<
1250         """Helper routine for parallel.deformsubpathitem: Returns a list of parameters where the curvature of normcurve is 1/distance"""
1251
1252         assert tstart < tend
1253         dist = self.dist_pt
1254
1255         # we _need_ to do this with the curvature, not with the radius
1256         # because the curvature is continuous at the straight line and the radius is not:
1257         # when passing from one slightly curved curve to the other with opposite curvature sign,
1258         # via the straight line, then the curvature changes its sign at curv=0, while the
1259         # radius changes its sign at +/-infinity
1260         # this causes instabilities for nearly straight curves
1261
1262         # include tstart and tend
1263         params = [tstart + i * (tend - tstart) / (self.lookforcurvatures - 1.0)
1264                   for i in range(self.lookforcurvatures)]
1265         curvs = normcurve.curvature_pt(params)
1266
1267         parampairs = list(zip(params[:-1], params[1:]))
1268         curvpairs = list(zip(curvs[:-1], curvs[1:]))
1269
1270         crossingparams = set()
1271         for parampair, curvpair in zip(parampairs, curvpairs):
1272             begparam, endparam = parampair
1273             begcurv, endcurv = curvpair
1274             begchange = begcurv*dist - 1
1275             endchange = endcurv*dist - 1
1276             if begchange*endchange < 0:
1277                 # the curvature crosses the value 1/dist
1278                 # get the parmeter value by linear interpolation:
1279                 middleparam = (
1280                   (begparam * abs(begchange) + endparam * abs(endchange)) /
1281                   (abs(begchange) + abs(endchange)))
1282                 try:
1283                     middleradius = 1/normcurve.curvature_pt([middleparam])[0]
1284                 except ArithmeticError:
1285                     raise InvalidParamException(middleparam)
1286
1287                 if abs(middleradius - dist) < epsilon or endparam-begparam < 1.0e-14:
1288                     # get the parmeter value by linear interpolation:
1289                     crossingparams.add(middleparam)
1290                 else:
1291                     # call recursively:
1292                     for x in self._distcrossingparameters(normcurve, epsilon, tstart=begparam, tend=endparam):
1293                         crossingparams.add(x)
1294             else:
1295                 if begchange == 0:
1296                     crossingparams.add(begparam)
1297                 if endchange == 0:
1298                     crossingparams.add(endparam)
1299
1300         return crossingparams
1301         # >>>
1302     def _get_angles(self, prev_nspitem, next_nspitem, epsilon): # <<<
1303         prev_rotation = prev_nspitem.rotation([1])[0]
1304         next_rotation = next_nspitem.rotation([0])[0]
1305         prev_tangent = prev_rotation.apply_pt(1, 0)
1306         prev_orthogo = prev_rotation.apply_pt(0, self.dist_pt) # points towards parallel path (prev_nspitem is on original path)
1307         next_tangent = next_rotation.apply_pt(1, 0)
1308         #sinangle = prev_tangent[0]*next_tangent[1] - prev_tangent[1]*next_tangent[0] # >0 for left-turning, <0 for right-turning
1309         cosangle = prev_tangent[0]*next_tangent[0] + prev_tangent[1]*next_tangent[1]
1310         proj = prev_orthogo[0]*next_tangent[0] + prev_orthogo[1]*next_tangent[1]
1311         is_straight = (cosangle > 0 and abs(proj) < epsilon)
1312         is_concave = (proj > 0)
1313         return prev_tangent, next_tangent, is_straight, is_concave
1314     # >>>
1315
1316     def rebuild_intersected_normpath(self, par_np, orig_np, par2orig, epsilon=None): # <<<
1317
1318         dist = self.dist_pt
1319         if epsilon is None:
1320             epsilon = orig_np.normsubpaths[0].epsilon
1321         eps_comparepairs = 10*epsilon
1322
1323         # calculate the self-intersections of the par_np
1324         forwardpairs, backwardpairs = self.normpath_selfintersections(par_np, epsilon, eps_comparepairs)
1325         # calculate the intersections of the par_np with the original path
1326         origintparams, orig_origintparams = self.normpath_origintersections(orig_np, par_np, epsilon)
1327         if not forwardpairs:
1328             if origintparams:
1329                 return normpath.normpath()
1330             else:
1331                 return par_np
1332
1333         # parameters at begin and end of subnormpaths:
1334         # omit those which start/end on the original path
1335         beginparams = []
1336         endparams = []
1337         testparams = origintparams + list(forwardpairs.keys()) + list(forwardpairs.values())
1338         for i, nsp in enumerate(par_np):
1339             beginparam = mynormpathparam(par_np, i, 0, 0)
1340             is_new = True
1341             for param in testparams:
1342                 if beginparam.is_equiv(param):
1343                     is_new = False
1344                     break
1345             if is_new:
1346                 beginparams.append(beginparam)
1347
1348             endparam = mynormpathparam(par_np, i, len(nsp)-1, 1)
1349             is_new = True
1350             for param in testparams:
1351                 if endparam.is_equiv(param):
1352                     is_new = False
1353                     break
1354             if is_new:
1355                 endparams.append(endparam)
1356         beginparams.sort()
1357         endparams.sort()
1358
1359         # we need a way to get the "next" param on the normpath
1360         # XXX why + beginparams + endparams ?
1361         allparams = list(forwardpairs.keys()) + list(backwardpairs.keys()) + origintparams + beginparams + endparams
1362         allparams.sort()
1363         done = {}
1364         for param in allparams:
1365             done[param] = False
1366         nextp = {}
1367         for i, param in enumerate(allparams[:-1]):
1368             nextp[param] = allparams[i+1]
1369         for endparam in endparams:
1370             if par_np[endparam.normsubpathindex].closed:
1371                 begparam = [p for p in allparams if p.normsubpathindex == endparam.normsubpathindex][0]
1372                 assert begparam.normsubpathitemindex == 0
1373                 assert begparam.normsubpathitemparam == 0
1374                 nextp[endparam] = begparam
1375             else:
1376                 nextp[endparam] = None
1377
1378         # exclude all intersections that are between the original and the parallel path:
1379         # See for example test/functional/test_deformer (parallel Z): There can
1380         # remain a little piece of the path (triangle) that lies between a lot
1381         # of intersection points. Simple intersection rules such as thoe in
1382         # trial_parampairs cannot exclude this piece.
1383         for param in forwardpairs:
1384             if done[param] or done[forwardpairs[param]]:
1385                 done[param] = done[forwardpairs[param]] = True
1386             elif self._between_paths(par_np.at_pt(param), par2orig, 4*epsilon):
1387                 done[param] = done[forwardpairs[param]] = True
1388         for param in beginparams + endparams:
1389             if self._between_paths(par_np.at_pt(param), par2orig, 4*epsilon):
1390                 done[param] = True
1391
1392         # visualize the intersection points: # <<<
1393         if self.debug is not None:
1394             for param1, param2 in forwardpairs.items():
1395                 point1, point2 = par_np.at([param1, param2])
1396                 if not done[param1]:
1397                     self.debug.fill(path.circle(point1[0], point1[1], 0.05), [color.rgb.red])
1398                 if not done[param2]:
1399                     self.debug.fill(path.circle(point2[0], point2[1], 0.03), [color.rgb.black])
1400             for param in origintparams:
1401                 #assert done[param]
1402                 point = par_np.at([param])[0]
1403                 self.debug.fill(path.circle(point[0], point[1], 0.05), [color.rgb.green])
1404             for i, nsp in enumerate(par_np):
1405               for j, nspi in enumerate(nsp):
1406                 x, y = nspi.at_pt([0.5])[0]
1407                 self.debug.text_pt(x, y, "{}/{}".format(i,j))#, [text.halign.center, text.vshift.mathaxis])
1408             print("aborted path intersection due to debug")
1409             return par_np
1410         # >>>
1411
1412         def ptype(param): # <<<
1413             if param in forwardpairs : return "fw with partner %s" % (forwardpairs[param])
1414             if param in backwardpairs : return "bw with partner %s" % (backwardpairs[param])
1415             if param in origintparams: return "orig"
1416             if param in beginparams: return "begin"
1417             if param in endparams: return "end"
1418         # >>>
1419         def trial_parampairs(startp): # <<<
1420             """Starting at startp, try to find a valid series of intersection parameters"""
1421             tried = {} # a local copy of done
1422             for param in allparams:
1423                 tried[param] = done[param]
1424
1425             previousp = startp
1426             currentp = nextp[previousp]
1427             result = []
1428
1429             while True:
1430                 # successful and unsuccessful termination conditions:
1431                 if tried[currentp]:
1432                     # we reached a branch that has already been treated
1433                     # ==> this is not a valid parallel path
1434                     return []
1435                 if currentp in origintparams:
1436                     # we cross the original path
1437                     # ==> this is not a valid parallel path
1438                     return []
1439                 if currentp in backwardpairs:
1440                     # we reached a branch that should be followed from another part
1441                     # ==> this is not a valid parallel path
1442                     return []
1443                 if currentp is startp:
1444                     # we have reached again the starting point on a closed subpath.
1445                     assert startp in beginparams
1446                     assert previousp in endparams
1447                     return result
1448                 if currentp in forwardpairs:
1449                     result.append((previousp, currentp))
1450                     if forwardpairs[currentp] is startp:
1451                         # we have found the same point as the startp (its pair partner)
1452                         return result
1453                     previousp = forwardpairs[currentp]
1454                 if currentp in endparams:
1455                     result.append((previousp, currentp))
1456                     if nextp[currentp] is None: # open subpath
1457                         # we have found the end of a non-closed subpath
1458                         return result
1459                     previousp = currentp # closed subpath
1460                 # follow the crossings on valid startpairs
1461                 tried[currentp] = True
1462                 tried[previousp] = True
1463                 currentp = nextp[previousp]
1464             assert False # never reach this point
1465         # >>>
1466
1467         # first the paths that start at the beginning of a subnormpath:
1468         result = normpath.normpath()
1469         # paths can start on subnormpaths or crossings where we can "get away":
1470         bwkeys = list(backwardpairs.keys())
1471         bwkeys.sort()
1472         for startp in beginparams + bwkeys:
1473             if done[startp]:
1474                 continue
1475
1476             # try to find a valid series of intersection points:
1477             parampairs = trial_parampairs(startp)
1478             if not parampairs:
1479                 continue
1480
1481             # collect all the pieces between parampairs:
1482             add_nsp = normpath.normsubpath(epsilon=epsilon)
1483             for begin, end in parampairs:
1484                 # check that trial_parampairs works correctly
1485                 assert begin is not end
1486                 for item in par_np[begin.normsubpathindex].segments(
1487                     [begin.normsubpathparam, end.normsubpathparam])[0].normsubpathitems:
1488                     # TODO: this should be obsolete with an improved intersection algorithm
1489                     #       guaranteeing epsilon
1490                     if add_nsp.normsubpathitems:
1491                         item = item.modifiedbegin_pt(*(add_nsp.atend_pt()))
1492                     add_nsp.append(item)
1493
1494                 done[begin] = True
1495                 done[end] = True
1496
1497             # close the path if necessary
1498             if add_nsp:
1499                 if ((parampairs[-1][-1] in forwardpairs and forwardpairs[parampairs[-1][-1]] is parampairs[0][0]) or
1500                     (parampairs[-1][-1] in endparams and parampairs[0][0] in beginparams and parampairs[0][0] is nextp[parampairs[-1][-1]])):
1501                     add_nsp.normsubpathitems[-1] = add_nsp.normsubpathitems[-1].modifiedend_pt(*add_nsp.atbegin_pt())
1502                     add_nsp.close()
1503
1504             result.extend([add_nsp])
1505
1506         return result
1507     # >>>
1508     def normpath_selfintersections(self, np, epsilon, eps_comparepairs): # <<<
1509
1510         """Returns all self-intersection points of normpath np.
1511
1512         This does not include the intersections of a single normcurve with itself,
1513         but all intersections of one normpathitem with a different one in the path.
1514         The intersection pairs are such that the parallel path can be continued
1515         from the first to the second parameter, but not vice-versa."""
1516
1517         dist = self.dist_pt
1518
1519         n = len(np)
1520         forwardpairs = {}
1521         for nsp_i in range(n):
1522             for nsp_j in range(nsp_i, n):
1523                 for nspitem_i in range(len(np[nsp_i])):
1524                     if nsp_j == nsp_i:
1525                         nspitem_j_range = list(range(nspitem_i+1, len(np[nsp_j])))
1526                     else:
1527                         nspitem_j_range = list(range(len(np[nsp_j])))
1528                     for nspitem_j in nspitem_j_range:
1529                         intsparams = np[nsp_i][nspitem_i].intersect(np[nsp_j][nspitem_j], epsilon)
1530                         if intsparams:
1531                             for intsparam_i, intsparam_j in intsparams:
1532                                 npp_i = mynormpathparam(np, nsp_i, nspitem_i, intsparam_i)
1533                                 npp_j = mynormpathparam(np, nsp_j, nspitem_j, intsparam_j)
1534
1535                                 # skip successive nsp-items
1536                                 if nsp_i == nsp_j:
1537                                     if nspitem_j == nspitem_i+1 and (npp_i.is_end_of_nspitem(epsilon) or npp_j.is_beg_of_nspitem(epsilon)):
1538                                         continue
1539                                     if np[nsp_i].closed and ((npp_i.is_beg_of_nsp(epsilon) and npp_j.is_end_of_nsp(epsilon)) or
1540                                                              (npp_j.is_beg_of_nsp(epsilon) and npp_i.is_end_of_nsp(epsilon))):
1541                                         continue
1542
1543                                 # correct the order of the pair, such that we can use it to continue on the path
1544                                 if not self._can_continue(npp_i, npp_j, epsilon):
1545                                     assert self._can_continue(npp_j, npp_i, epsilon)
1546                                     npp_i, npp_j = npp_j, npp_i
1547
1548                                 # if the intersection is between two nsp-items, take the smallest -> largest
1549                                 npp_i = npp_i.smaller_equiv(5*epsilon)
1550                                 npp_j = npp_j.larger_equiv(5*epsilon)
1551
1552                                 # because of the above change of npp_ij, and because there may be intersections between nsp-items,
1553                                 # it may happen that we try to insert two times the same pair
1554                                 if self._skip_intersection_doublet(npp_i, npp_j, forwardpairs, eps_comparepairs):
1555                                     continue
1556                                 forwardpairs[npp_i] = npp_j
1557
1558         # this is partially done in _skip_intersection_doublet
1559         #forwardpairs = self._elim_intersection_doublets(forwardpairs, eps_comparepairs)
1560         # create the reverse mapping
1561         backwardpairs = {}
1562         for p, q in forwardpairs.items():
1563             backwardpairs[q] = p
1564         return forwardpairs, backwardpairs
1565
1566     # >>>
1567     def normpath_origintersections(self, orig_np, par_np, epsilon): # <<<
1568         """return all intersection points of the original path and the parallel path"""
1569
1570         # this code became necessary with introduction of mynormpathparam
1571         params = []
1572         oparams = []
1573         for nsp_i in range(len(orig_np)):
1574             for nsp_j in range(len(par_np)):
1575                 for nspitem_i in range(len(orig_np[nsp_i])):
1576                     for nspitem_j in range(len(par_np[nsp_j])):
1577                         intsparams = orig_np[nsp_i][nspitem_i].intersect(par_np[nsp_j][nspitem_j], epsilon)
1578                         if intsparams:
1579                             for intsparam_i, intsparam_j in intsparams:
1580                                 npp_i = mynormpathparam(orig_np, nsp_i, nspitem_i, intsparam_i)
1581                                 npp_j = mynormpathparam(par_np, nsp_j, nspitem_j, intsparam_j)
1582
1583                                 oparams.append(npp_i)
1584                                 params.append(npp_j)
1585         return params, oparams
1586     # >>>
1587     def _can_continue(self, param1, param2, epsilon=None): # <<<
1588         """Test whether the parallel path can be continued at the param-pair (param1, param2)"""
1589         par_np = param1.normpath
1590         if epsilon is None:
1591             epsilon = par_np[0].epsilon
1592
1593         rot1, rot2 = par_np.rotation([param1, param2])
1594         orth1 = rot1.apply_pt(0, self.dist_pt) # directs away from original path (as seen from parallel path)
1595         tang2 = rot2.apply_pt(1, 0)
1596
1597         # the self-intersection is valid if the tangents
1598         # point into the correct direction or, for parallel tangents,
1599         # if the curvature is such that the on-going path does not
1600         # enter the region defined by dist
1601         proj = orth1[0]*tang2[0] + orth1[1]*tang2[1]
1602         if abs(proj) > epsilon: # the curves are not parallel
1603             # tang2 must go away from the original path
1604             return (proj > 0)
1605
1606         # tang1 and tang2 are parallel.
1607         curv1, curv2 = par_np.curvature_pt([param1, param2])
1608
1609         # We need to treat also cases where the params are nspitem-endpoints.
1610         # There, we know that the tangents are continuous, but the curvature is
1611         # not necessarily continuous. We have to test whether the curve *after*
1612         # param2 has curvature such that it enters the forbidden side of the
1613         # curve after param1
1614         if param1.is_end_of_nspitem(epsilon):
1615             curv1 = par_np.curvature_pt([param1.larger_equiv(epsilon)])[0]
1616         if param2.is_end_of_nspitem(epsilon):
1617             curv2 = par_np.curvature_pt([param2.larger_equiv(epsilon)])[0]
1618
1619         tang1 = rot1.apply_pt(1, 0)
1620         running_back = (tang1[0]*tang2[0] + tang1[1]*tang2[1] < 0)
1621         if running_back:
1622             # the second curve is running "back" -- the curvature sign appears to be switched
1623             curv2 = -curv2
1624             # endpoints of normsubpaths must be treated differently:
1625
1626         if (not running_back) and param1.is_end_of_nsp(epsilon):
1627             return True
1628
1629         if curv1 == curv2:
1630             raise IntersectionError("Cannot determine whether curves intersect (parallel and equally curved)")
1631
1632         if self.dist_pt > 0:
1633             return (curv2 > curv1)
1634         else:
1635             return (curv2 < curv1)
1636     # >>>
1637     def _skip_intersection_doublet(self, npp_i, npp_j, parampairs, epsilon): # <<<
1638         # An intersection point that lies exactly between two nsp-items can occur twice or more
1639         # times if we calculate all mutual intersections. We should take only
1640         # one such parameter pair, namely the one with smallest first and
1641         # largest last param.
1642         result = False
1643         delete_keys = []
1644         delete_values = []
1645         # TODO: improve complexity?
1646         for pi, pj in parampairs.items():
1647             if npp_i.is_equiv(pi, epsilon) and npp_j.is_equiv(pj, epsilon):
1648                 #print("double pair: ", npp_i, npp_j, pi, pj)
1649                 #print("... replacing ", pi, parampairs[pi], "by", min(npp_i, pi), max(npp_j, pj))
1650                 delete_keys.append(pi)
1651                 delete_values.append(pj)
1652                 result = True # we have already added this one
1653         newkey = min([npp_i] + delete_keys)
1654         newval = max([npp_j] + delete_values)
1655         for pi in delete_keys:
1656             del parampairs[pi]
1657         parampairs[newkey] = newval
1658         return result
1659     # >>>
1660     def _elim_intersection_doublets(self, parampairs, epsilon): # <<<
1661         # It may always happen that three intersections coincide. (It will
1662         # occur often with degenerate distances for technical designs such as
1663         # those used in microfluidics). We then have two equivalent pairs in our
1664         # forward list, and we must throw away one of them.
1665         # One of them is indeed forbidden by the _can_continue of the other.
1666
1667         # TODO implement this
1668
1669         keys = list(parampairs.keys())
1670         n = len(keys)
1671         for i in range(n):
1672             start = "equivalent pairs\n"
1673             for j in range(i+1, n):
1674                 key1, key2 = keys[i], keys[j]
1675                 npp1 = parampairs[key1]
1676                 npp2 = parampairs[key2]
1677                 #assert key1.is_equiv(npp1, epsilon)
1678                 #if not key2.is_equiv(npp2, epsilon):
1679                 #    np = key2.normpath
1680                 #    print(np.at_pt(key2), np.at_pt(npp2), _length_pt(np, key2, npp2)/epsilon)
1681                 #assert key2.is_equiv(npp2, epsilon)
1682                 if ((key1.is_equiv(key2, epsilon) and npp1.is_equiv(npp2, epsilon)) or
1683                     (key1.is_equiv(npp2, epsilon) and npp1.is_equiv(key2, epsilon))):
1684                     print(start,"pair: ", key1, npp1, " and ", key2, npp2)
1685                     start = ""
1686             if not start:
1687                 print()
1688         return parampairs
1689     # >>>
1690     def _between_paths(self, pos, par2orig, epsilon): # <<<
1691         """Tests whether the given point (pos) is found in the forbidden zone between an original and a parallel nsp-item (these are in par2orig)
1692
1693         The test uses epsilon close to the original/parallel path, and sharp comparison at their ends."""
1694         dist = self.dist_pt
1695         for par_nspitem in par2orig:
1696             origobj = par2orig[par_nspitem]
1697             if isinstance(origobj, normpath.normline_pt):
1698                 rot = origobj.rotation([0])[0]
1699                 t, s = intersection(pos, origobj.atbegin_pt(), rot.apply_pt(0, mathutils.sign(dist)), rot.apply_pt(origobj.arclen_pt(epsilon), 0))
1700                 if 0 <= s <= 1 and -abs(dist)+epsilon < t < -epsilon:
1701                     return True
1702             elif isinstance(origobj, normpath.normcurve_pt):
1703                 # TODO: implement this
1704                 # TODO: pre-sort par2orig as a list to fasten up this code
1705                 pass
1706             else:
1707                 cx, cy = origobj
1708                 if math.hypot(pos[0]-cx, pos[1]-cy) < abs(dist) - epsilon:
1709                     if self.dist_pt > 0: # running around (cx,cy) in the negative sense (see _path_around_corner)
1710                         x0, y0 = par_nspitem.atend_pt()
1711                         x1, y1 = par_nspitem.atbegin_pt()
1712                     else: # running around (cx,cy) in the positive sense
1713                         x0, y0 = par_nspitem.atbegin_pt()
1714                         x1, y1 = par_nspitem.atend_pt()
1715                     t0, s0 = intersection(pos, (cx, cy), (-y0+cy, x0-cx), (x0-cx, y0-cy))
1716                     t1, s1 = intersection(pos, (cx, cy), ( y1-cy,-x1+cx), (x1-cx, y1-cy))
1717                     if t0 <= 0 and s0 >= 0 and t1 <= 0 and s1 >= 0:
1718                         return True
1719         return False
1720     # >>>
1721
1722 # >>>
1723
1724 parallel.clear = attr.clearclass(parallel)
1725
1726 class linesmoothed(baseclasses.deformer): # <<<
1727
1728     def __init__(self, tension=1, atleast=False, lcurl=1, rcurl=1):
1729         """Tension and atleast control the tension of the replacement curves.
1730         l/rcurl control the curlynesses at (possible) endpoints. If a curl is
1731         set to None, the angle is taken from the original path."""
1732         if atleast:
1733             self.tension = -abs(tension)
1734         else:
1735             self.tension = abs(tension)
1736         self.lcurl = lcurl
1737         self.rcurl = rcurl
1738
1739     def __call__(self, tension=_marker, atleast=_marker, lcurl=_marker, rcurl=_marker):
1740         if tension is _marker:
1741             tension = self.tension
1742         if atleast is _marker:
1743             atleast = (self.tension < 0)
1744         if lcurl is _marker:
1745             lcurl = self.lcurl
1746         if rcurl is _marker:
1747             rcurl = self.rcurl
1748         return linesmoothed(tension, atleast, lcurl, rcurl)
1749
1750     def deform(self, basepath):
1751         newnp = normpath.normpath()
1752         for nsp in basepath.normpath().normsubpaths:
1753             newnp += self.deformsubpath(nsp)
1754         return newnp
1755
1756     def deformsubpath(self, nsp):
1757         from .metapost import path as mppath
1758         """Returns a path/normpath from the points in the given normsubpath"""
1759         # TODO: epsilon ?
1760         knots = []
1761
1762         # first point
1763         x_pt, y_pt = nsp.atbegin_pt()
1764         if nsp.closed:
1765             knots.append(mppath.smoothknot_pt(x_pt, y_pt))
1766         elif self.lcurl is None:
1767             rot = nsp.rotation([0])[0]
1768             dx, dy = rot.apply_pt(1, 0)
1769             angle = math.atan2(dy, dx)
1770             knots.append(mppath.beginknot_pt(x_pt, y_pt, angle=angle))
1771         else:
1772             knots.append(mppath.beginknot_pt(x_pt, y_pt, curl=self.lcurl))
1773
1774         # intermediate points:
1775         for npelem in nsp[:-1]:
1776             knots.append(mppath.tensioncurve(self.tension))
1777             knots.append(mppath.smoothknot_pt(*npelem.atend_pt()))
1778
1779         # last point
1780         knots.append(mppath.tensioncurve(self.tension))
1781         x_pt, y_pt = nsp.atend_pt()
1782         if nsp.closed:
1783             pass
1784         elif self.rcurl is None:
1785             rot = nsp.rotation([len(nsp)])[0]
1786             dx, dy = rot.apply_pt(1, 0)
1787             angle = math.atan2(dy, dx)
1788             knots.append(mppath.endknot_pt(x_pt, y_pt, angle=angle))
1789         else:
1790             knots.append(mppath.endknot_pt(x_pt, y_pt, curl=self.rcurl))
1791
1792         return mppath.path(knots)
1793 # >>>
1794
1795 linesmoothed.clear = attr.clearclass(linesmoothed)
1796
1797
1798 # vim:foldmethod=marker:foldmarker=<<<,>>>