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