pyx/deformer.py

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