pyx/deformer.py

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