pyx/deformer.py

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