pyx/deformer.py

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