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