1# -*- encoding: utf-8 -*-
2#
3#
4# Copyright (C) 2003-2013 Michael Schindler <m-schindler@users.sourceforge.net>
5# Copyright (C) 2003-2005 André Wobst <wobsta@pyx-project.org>
6#
7# This file is part of PyX (https://pyx-project.org/).
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
23import functools, logging, math
24from . import attr, baseclasses, mathutils, path, normpath, unit, color
25
26normpath.invalid = 175e175 # Just a very crude workaround to get the code running again. normpath.invalid does not exist anymore.
27
28logger = logging.getLogger("pyx")
29
30# specific exception for an invalid parameterization point
31# used in parallel
32class InvalidParamException(Exception):
33
34    def __init__(self, param):
35        self.normsubpathitemparam = param
36
37# error raised in parallel if we are trying to get badly defined intersections
38class IntersectionError(Exception): pass
39
40# None has a meaning in linesmoothed
41class _marker: pass
42
43class inf_curvature: pass
44
45def curvescontrols_from_endlines_pt(B, tangent1, tangent2, r1, r2, softness): # <<<
46    # calculates the parameters for two bezier curves connecting two lines (curvature=0)
47    # starting at B - r1*tangent1
48    # ending at   B + r2*tangent2
49    #
50    # Takes the corner B
51    # and two tangent vectors heading to and from B
52    # and two radii r1 and r2:
53    # All arguments must be in Points
54    # Returns the seven control points of the two bezier curves:
55    #  - start d1
56    #  - control points g1 and f1
57    #  - midpoint e
58    #  - control points f2 and g2
59    #  - endpoint d2
60
61    # make direction vectors d1: from B to A
62    #                        d2: from B to C
63    d1 = -tangent1[0] / math.hypot(*tangent1), -tangent1[1] / math.hypot(*tangent1)
64    d2 =  tangent2[0] / math.hypot(*tangent2),  tangent2[1] / math.hypot(*tangent2)
65
66    # 0.3192 has turned out to be the maximum softness available
67    # for straight lines ;-)
68    f = 0.3192 * softness
69    g = (15.0 * f + math.sqrt(-15.0*f*f + 24.0*f))/12.0
70
71    # make the control points of the two bezier curves
72    f1 = B[0] + f * r1 * d1[0], B[1] + f * r1 * d1[1]
73    f2 = B[0] + f * r2 * d2[0], B[1] + f * r2 * d2[1]
74    g1 = B[0] + g * r1 * d1[0], B[1] + g * r1 * d1[1]
75    g2 = B[0] + g * r2 * d2[0], B[1] + g * r2 * d2[1]
76    d1 = B[0] +     r1 * d1[0], B[1] +     r1 * d1[1]
77    d2 = B[0] +     r2 * d2[0], B[1] +     r2 * d2[1]
78    e  = 0.5 * (f1[0] + f2[0]), 0.5 * (f1[1] + f2[1])
79
80    return (d1, g1, f1, e, f2, g2, d2)
81# >>>
82
83def controldists_from_endgeometry_pt(A, B, tangA, tangB, curvA, curvB, allownegative=False, curv_epsilon=1.0e-8): # <<<
84
85    """For a curve with given tangents and curvatures at the endpoints this gives the distances between the controlpoints
86
87    This helper routine returns a list of two distances between the endpoints and the
88    corresponding control points of a (cubic) bezier curve that has
89    prescribed tangents tangentA, tangentB and curvatures curvA, curvB at the
90    end points.
91
92    Note: The returned distances are not always positive.
93          But only positive values are geometrically correct, so please check!
94          The outcome is sorted so that the first entry is expected to be the
95          most reasonable one
96    """
97    debug = 0
98
99    def test_divisions(T, D, E, AB, curvA, curvB, debug):# <<<
100        small = AB * 1.0e-4 # at infinite curvature, avoid setting control points exactly on the startpoint
101                            # TODO: is this consistent with the avoiding of curv=inf in normpath?
102        arbitrary = AB * 0.33 # at zero curvature, we know nothing about a or b
103
104        def is_zero(x):
105            return abs(x) < curv_epsilon
106            # the following gave different results for forward/reverse paths
107            # in test/functional/test_deformer parallel G
108            #try:
109            #    1.0 / x
110            #except ZeroDivisionError:
111            #    return True
112            #return False
113
114
115        if is_zero(T):
116
117            if curvA is inf_curvature:
118               a = small
119               if curvB is inf_curvature:
120                   b = small
121               elif is_zero(curvB):
122                    assert abs(E) < 1.0e-10
123                    b = arbitrary
124               else:
125                    b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
126            elif is_zero(curvA):
127                assert abs(D) < 1.0e-10
128                a = arbitrary
129                if curvB is inf_curvature:
130                    b = small
131                elif is_zero(curvB):
132                    assert abs(E) < 1.0e-10
133                    b = arbitrary
134                else:
135                    b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
136            else:
137                a = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA)
138                if curvB is inf_curvature:
139                    b = small
140                elif is_zero(curvB):
141                    assert abs(E) < 1.0e-10
142                    b = arbitrary
143                else:
144                    b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
145
146        else:
147            if curvA is inf_curvature:
148               a = small
149               if curvB is inf_curvature:
150                   b = small
151               elif is_zero(curvB):
152                   b = arbitrary
153               else:
154                   b1 = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
155                   b2 = D / T
156                   if abs(b1) < abs(b2):
157                       b = b1
158                   else:
159                       b = b2
160            elif curvB is inf_curvature:
161               b = small
162               if is_zero(curvA):
163                   a = arbitrary
164               else:
165                   a1 = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA)
166                   a2 = E / T
167                   if abs(a1) < abs(a2):
168                       a = a1
169                   else:
170                       a = a2
171            elif is_zero(curvA):
172                b = D / T
173                a = (E - 1.5*curvB*b*abs(b)) / T
174            elif is_zero(curvB):
175                a = E / T
176                b = (D - 1.5*curvA*a*abs(a)) / T
177            else:
178                return []
179
180        if debug:
181            print("fallback with exact zero value")
182        return [(a, b)]
183    # >>>
184    def fallback_smallT(T, D, E, AB, curvA, curvB, threshold, debug):# <<<
185        a = math.sqrt(abs(D / (1.5 * curvA))) * mathutils.sign(D*curvA)
186        b = math.sqrt(abs(E / (1.5 * curvB))) * mathutils.sign(E*curvB)
187        q1 = min(abs(1.5*a*a*curvA), abs(D))
188        q2 = min(abs(1.5*b*b*curvB), abs(E))
189        if (a >= 0 and b >= 0 and
190            abs(b*T) < threshold * q1 and abs(1.5*a*abs(a)*curvA - D) < threshold * q1 and
191            abs(a*T) < threshold * q2 and abs(1.5*b*abs(b)*curvB - E) < threshold * q2):
192            if debug:
193                print("fallback with T approx 0")
194            return [(a, b)]
195        return []
196    # >>>
197    def fallback_smallcurv(T, D, E, AB, curvA, curvB, threshold, debug):# <<<
198        result = []
199
200        # is curvB approx zero?
201        a = E / T
202        b = (D - 1.5*curvA*a*abs(a)) / T
203        if (a >= 0 and b >= 0 and
204            abs(1.5*b*b*curvB) < threshold * min(abs(a*T), abs(E)) and
205            abs(a*T - E) < threshold * min(abs(a*T), abs(E))):
206            if debug:
207                print("fallback with curvB approx 0")
208            result.append((a, b))
209
210        # is curvA approx zero?
211        b = D / T
212        a = (E - 1.5*curvB*b*abs(b)) / T
213        if (a >= 0 and b >= 0 and
214            abs(1.5*a*a*curvA) < threshold * min(abs(b*T), abs(D)) and
215            abs(b*T - D) < threshold * min(abs(b*T), abs(D))):
216            if debug:
217                print("fallback with curvA approx 0")
218            result.append((a, b))
219
220        return result
221    # >>>
222    def findnearest(x, ys): # <<<
223        I = 0
224        Y = ys[I]
225        mindist = abs(x - Y)
226
227        # find the value in ys which is nearest to x
228        for i, y in enumerate(ys[1:]):
229            dist = abs(x - y)
230            if dist < mindist:
231                I, Y, mindist = i, y, dist
232
233        return I, Y
234    # >>>
235
236    # some shortcuts
237    T = tangA[0] * tangB[1] - tangA[1] * tangB[0]
238    D = tangA[0] * (B[1]-A[1]) - tangA[1] * (B[0]-A[0])
239    E = tangB[0] * (A[1]-B[1]) - tangB[1] * (A[0]-B[0])
240    AB = math.hypot(A[0] - B[0], A[1] - B[1])
241
242    # try if one of the prefactors is exactly zero
243    testsols = test_divisions(T, D, E, AB, curvA, curvB, debug)
244    if testsols:
245        return testsols
246
247    # The general case:
248    # we try to find all the zeros of the decoupled 4th order problem
249    # for the combined problem:
250    # The control points of a cubic Bezier curve are given by a, b:
251    #     A, A + a*tangA, B - b*tangB, B
252    # for the derivation see /design/beziers.tex
253    #     0 = 1.5 a |a| curvA + b * T - D
254    #     0 = 1.5 b |b| curvB + a * T - E
255    # because of the absolute values we get several possibilities for the signs
256    # in the equation. We test all signs, also the invalid ones!
257    if allownegative:
258        signs = [(+1, +1), (-1, +1), (+1, -1), (-1, -1)]
259    else:
260        signs = [(+1, +1)]
261
262    candidates_a = []
263    candidates_b = []
264    for sign_a, sign_b in signs:
265        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)
266        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)
267        candidates_a += [root for root in mathutils.realpolyroots(*coeffs_a) if sign_a*root >= 0]
268        candidates_b += [root for root in mathutils.realpolyroots(*coeffs_b) if sign_b*root >= 0]
269    solutions = []
270    if candidates_a and candidates_b:
271        for a in candidates_a:
272            i, b = findnearest((D - 1.5*curvA*a*abs(a))/T, candidates_b)
273            solutions.append((a, b))
274
275    # try if there is an approximate solution
276    for thr in [1.0e-2, 1.0e-1]:
277        if not solutions:
278            solutions = fallback_smallT(T, D, E, AB, curvA, curvB, thr, debug)
279        if not solutions:
280            solutions = fallback_smallcurv(T, D, E, AB, curvA, curvB, thr, debug)
281
282    # sort the solutions: the more reasonable values at the beginning
283    def mycmp(x,y): # <<<
284        # first the pairs that are purely positive, then all the pairs with some negative signs
285        # inside the two sets: sort by magnitude
286        sx = (x[0] > 0 and x[1] > 0)
287        sy = (y[0] > 0 and y[1] > 0)
288
289        # experimental stuff:
290        # what criterion should be used for sorting ?
291        #
292        #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)
293        #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)
294        # # For each equation, a value like
295        # #   abs(1.5*curvA*y[0]*abs(y[0]) + y[1]*T - D) / abs(curvA*(D - y[1]*T))
296        # # indicates how good the solution is. In order to avoid the division,
297        # # we here multiply with all four denominators:
298        # 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)) ),
299        #            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)) ))
300        # 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)) ),
301        #            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)) ))
302        #errx = (abs(curvA*x[0]) - 1.0)**2 + (abs(curvB*x[1]) - 1.0)**2
303        #erry = (abs(curvA*y[0]) - 1.0)**2 + (abs(curvB*y[1]) - 1.0)**2
304
305        errx = x[0]**2 + x[1]**2
306        erry = y[0]**2 + y[1]**2
307
308        if sx == 1 and sy == 1:
309            # try to use longer solutions if there are any crossings in the control-arms
310            # the following combination yielded fewest sorting errors in test_bezier.py
311            t, s = intersection(A, B, tangA, tangB)
312            t, s = abs(t), abs(s)
313            if (t > 0 and t < x[0] and s > 0 and s < x[1]):
314                if (t > 0 and t < y[0] and s > 0 and s < y[1]):
315                    # use the shorter one
316                    return int(errx > erry) - int(errx < erry)
317                else:
318                    # use the longer one
319                    return -1
320            else:
321                if (t > 0 and t < y[0] and s > 0 and s < y[1]):
322                    # use the longer one
323                    return 1
324                else:
325                    # use the shorter one
326                    return int(errx > erry) - int(errx < erry)
327            #return cmp(x[0]**2 + x[1]**2, y[0]**2 + y[1]**2)
328        else:
329            return int(sy > sx) - int(sy < sx)
330    # >>>
331    solutions.sort(key=functools.cmp_to_key(mycmp))
332
333    return solutions
334# >>>
335
336def normcurve_from_endgeometry_pt(A, B, tangA, tangB, curvA, curvB): # <<<
337    a, b = controldists_from_endgeometry_pt(A, B, tangA, tangB, curvA, curvB)[0]
338    return normpath.normcurve_pt(A[0], A[1],
339        A[0] + a * tangA[0], A[1] + a * tangA[1],
340        B[0] - b * tangB[0], B[1] - b * tangB[1], B[0], B[1])
341    # >>>
342
343def intersection(A, D, tangA, tangD): # <<<
344
345    """returns the intersection parameters of two evens
346
347    they are defined by:
348      x(t) = A + t * tangA
349      x(s) = D + s * tangD
350    """
351    det = -tangA[0] * tangD[1] + tangA[1] * tangD[0]
352    try:
353        1.0 / det
354    except ArithmeticError:
355        return None, None
356
357    DA = D[0] - A[0], D[1] - A[1]
358
359    t = (-tangD[1]*DA[0] + tangD[0]*DA[1]) / det
360    s = (-tangA[1]*DA[0] + tangA[0]*DA[1]) / det
361
362    return t, s
363# >>>
364
365class cycloid(baseclasses.deformer): # <<<
366    """Wraps a cycloid around a path.
367
368    The outcome looks like a spring with the originalpath as the axis.
369    radius: radius of the cycloid
370    halfloops:  number of halfloops
371    skipfirst/skiplast: undeformed end lines of the original path
372    curvesperhloop:
373    sign: start left (1) or right (-1) with the first halfloop
374    turnangle: angle of perspective on a (3D) spring
375               turnangle=0 will produce a sinus-like cycloid,
376               turnangle=90 will procude a row of connected circles
377
378    """
379
380    def __init__(self, radius=0.5*unit.t_cm, halfloops=10,
381    skipfirst=1*unit.t_cm, skiplast=1*unit.t_cm, curvesperhloop=3, sign=1, turnangle=45):
382        self.skipfirst = skipfirst
383        self.skiplast = skiplast
384        self.radius = radius
385        self.halfloops = halfloops
386        self.curvesperhloop = curvesperhloop
387        self.sign = sign
388        self.turnangle = turnangle
389
390    def __call__(self, radius=None, halfloops=None,
391    skipfirst=None, skiplast=None, curvesperhloop=None, sign=None, turnangle=None):
392        if radius is None:
393            radius = self.radius
394        if halfloops is None:
395            halfloops = self.halfloops
396        if skipfirst is None:
397            skipfirst = self.skipfirst
398        if skiplast is None:
399            skiplast = self.skiplast
400        if curvesperhloop is None:
401            curvesperhloop = self.curvesperhloop
402        if sign is None:
403            sign = self.sign
404        if turnangle is None:
405            turnangle = self.turnangle
406
407        return cycloid(radius=radius, halfloops=halfloops, skipfirst=skipfirst, skiplast=skiplast,
408                       curvesperhloop=curvesperhloop, sign=sign, turnangle=turnangle)
409
410    def deform(self, basepath):
411        resultnormsubpaths = [self.deformsubpath(nsp) for nsp in basepath.normpath().normsubpaths]
412        return normpath.normpath(resultnormsubpaths)
413
414    def deformsubpath(self, normsubpath):
415
416        skipfirst = abs(unit.topt(self.skipfirst))
417        skiplast = abs(unit.topt(self.skiplast))
418        radius = abs(unit.topt(self.radius))
419        turnangle = math.radians(self.turnangle)
420        sign = mathutils.sign(self.sign)
421
422        cosTurn = math.cos(turnangle)
423        sinTurn = math.sin(turnangle)
424
425        # make list of the lengths and parameters at points on normsubpath
426        # where we will add cycloid-points
427        totlength = normsubpath.arclen_pt()
428        if totlength <= skipfirst + skiplast + 2*radius*sinTurn:
429            logger.warning("normsubpath is too short for deformation with cycloid -- skipping...")
430            return normsubpath
431
432        # parameterization is in rotation-angle around the basepath
433        # differences in length, angle ... between two basepoints
434        # and between basepoints and controlpoints
435        Dphi = math.pi / self.curvesperhloop
436        phis = [i * Dphi for i in range(self.halfloops * self.curvesperhloop + 1)]
437        DzDphi = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * math.pi * cosTurn)
438        # Dz = (totlength - skipfirst - skiplast - 2*radius*sinTurn) * 1.0 / (self.halfloops * self.curvesperhloop * cosTurn)
439        # zs = [i * Dz for i in range(self.halfloops * self.curvesperhloop + 1)]
440        # from path._arctobcurve:
441        # optimal relative distance along tangent for second and third control point
442        L = 4 * radius * (1 - math.cos(Dphi/2)) / (3 * math.sin(Dphi/2))
443
444        # Now the transformation of z into the turned coordinate system
445        Zs = [ skipfirst + radius*sinTurn # here the coordinate z starts
446             - sinTurn*radius*math.cos(phi) + cosTurn*DzDphi*phi # the transformed z-coordinate
447             for phi in phis]
448        params = normsubpath._arclentoparam_pt(Zs)[0]
449
450        # get the positions of the splitpoints in the cycloid
451        points = []
452        for phi, param in zip(phis, params):
453            # the cycloid is a circle that is stretched along the normsubpath
454            # here are the points of that circle
455            basetrafo = normsubpath.trafo([param])[0]
456
457            # The point on the cycloid, in the basepath's local coordinate system
458            baseZ, baseY = 0, radius*math.sin(phi)
459
460            # The tangent there, also in local coords
461            tangentX = -cosTurn*radius*math.sin(phi) + sinTurn*DzDphi
462            tangentY = radius*math.cos(phi)
463            tangentZ = sinTurn*radius*math.sin(phi) + DzDphi*cosTurn
464            norm = math.sqrt(tangentX*tangentX + tangentY*tangentY + tangentZ*tangentZ)
465            tangentY, tangentZ = tangentY/norm, tangentZ/norm
466
467            # Respect the curvature of the basepath for the cycloid's curvature
468            # XXX this is only a heuristic, not a "true" expression for
469            #     the curvature in curved coordinate systems
470            try:
471                pathradius = 1/normsubpath.curvature_pt([param])[0]
472            except ArithmeticError:
473                factor = 1
474            else:
475                factor = (pathradius - baseY) / pathradius
476                factor = abs(factor)
477            l = L * factor
478
479            # The control points prior and after the point on the cycloid
480            preeZ, preeY = baseZ - l * tangentZ, baseY - l * tangentY
481            postZ, postY = baseZ + l * tangentZ, baseY + l * tangentY
482
483            # Now put everything at the proper place
484            points.append(basetrafo.apply_pt(preeZ, sign * preeY) +
485                          basetrafo.apply_pt(baseZ, sign * baseY) +
486                          basetrafo.apply_pt(postZ, sign * postY))
487
488        if len(points) <= 1:
489            logger.warning("normsubpath is too short for deformation with cycloid -- skipping...")
490            return normsubpath
491
492        # Build the path from the pointlist
493        # containing (control x 2,  base x 2, control x 2)
494        if skipfirst > normsubpath.epsilon:
495            normsubpathitems = normsubpath.segments([0, params[0]])[0]
496            normsubpathitems.append(normpath.normcurve_pt(*(points[0][2:6] + points[1][0:4])))
497        else:
498            normsubpathitems = [normpath.normcurve_pt(*(points[0][2:6] + points[1][0:4]))]
499        for i in range(1, len(points)-1):
500            normsubpathitems.append(normpath.normcurve_pt(*(points[i][2:6] + points[i+1][0:4])))
501        if skiplast > normsubpath.epsilon:
502            for nsp in normsubpath.segments([params[-1], len(normsubpath)]):
503                normsubpathitems.extend(nsp.normsubpathitems)
504
505        # That's it
506        return normpath.normsubpath(normsubpathitems, epsilon=normsubpath.epsilon)
507# >>>
508
509cycloid.clear = attr.clearclass(cycloid)
510
511class cornersmoothed(baseclasses.deformer): # <<<
512
513    """Bends corners in a normpath.
514
515    This decorator replaces corners in a normpath with bezier curves. There are two cases:
516    - If the corner lies between two lines, _two_ bezier curves will be used
517      that are highly optimized to look good (their curvature is to be zero at the ends
518      and has to have zero derivative in the middle).
519      Additionally, it can controlled by the softness-parameter.
520    - If the corner lies between curves then _one_ bezier is used that is (except in some
521      special cases) uniquely determined by the tangents and curvatures at its end-points.
522      In some cases it is necessary to use only the absolute value of the curvature to avoid a
523      cusp-shaped connection of the new bezier to the old path. In this case the use of
524      "obeycurv=0" allows the sign of the curvature to switch.
525    - The radius argument gives the arclength-distance of the corner to the points where the
526      old path is cut and the beziers are inserted.
527    - Path elements that are too short (shorter than the radius) are skipped
528    """
529
530    def __init__(self, radius, softness=1, obeycurv=0, relskipthres=0.01):
531        self.radius = radius
532        self.softness = softness
533        self.obeycurv = obeycurv
534        self.relskipthres = relskipthres
535
536    def __call__(self, radius=None, softness=None, obeycurv=None, relskipthres=None):
537        if radius is None:
538            radius = self.radius
539        if softness is None:
540            softness = self.softness
541        if obeycurv is None:
542            obeycurv = self.obeycurv
543        if relskipthres is None:
544            relskipthres = self.relskipthres
545        return cornersmoothed(radius=radius, softness=softness, obeycurv=obeycurv, relskipthres=relskipthres)
546
547    def deform(self, basepath):
548        return normpath.normpath([self.deformsubpath(normsubpath)
549                              for normsubpath in basepath.normpath().normsubpaths])
550
551    def deformsubpath(self, normsubpath):
552        radius_pt = unit.topt(self.radius)
553        epsilon = normsubpath.epsilon
554
555        # remove too short normsubpath items (shorter than self.relskipthres*radius_pt or epsilon)
556        pertinentepsilon = max(epsilon, self.relskipthres*radius_pt)
557        pertinentnormsubpath = normpath.normsubpath(normsubpath.normsubpathitems,
558                                                epsilon=pertinentepsilon)
559        pertinentnormsubpath.flushskippedline()
560        pertinentnormsubpathitems = pertinentnormsubpath.normsubpathitems
561
562        # calculate the splitting parameters for the pertinentnormsubpathitems
563        arclens_pt = []
564        params = []
565        for pertinentnormsubpathitem in pertinentnormsubpathitems:
566            arclen_pt = pertinentnormsubpathitem.arclen_pt(epsilon)
567            arclens_pt.append(arclen_pt)
568            l1_pt = min(radius_pt, 0.5*arclen_pt)
569            l2_pt = max(0.5*arclen_pt, arclen_pt - radius_pt)
570            params.append(pertinentnormsubpathitem.arclentoparam_pt([l1_pt, l2_pt], epsilon))
571
572        # handle the first and last pertinentnormsubpathitems for a non-closed normsubpath
573        if not normsubpath.closed:
574            l1_pt = 0
575            l2_pt = max(0, arclens_pt[0] - radius_pt)
576            params[0] = pertinentnormsubpathitems[0].arclentoparam_pt([l1_pt, l2_pt], epsilon)
577            l1_pt = min(radius_pt, arclens_pt[-1])
578            l2_pt = arclens_pt[-1]
579            params[-1] = pertinentnormsubpathitems[-1].arclentoparam_pt([l1_pt, l2_pt], epsilon)
580
581        newnormsubpath = normpath.normsubpath(epsilon=normsubpath.epsilon)
582        for i in range(len(pertinentnormsubpathitems)):
583            this = i
584            next = (i+1) % len(pertinentnormsubpathitems)
585            thisparams = params[this]
586            nextparams = params[next]
587            thisnormsubpathitem = pertinentnormsubpathitems[this]
588            nextnormsubpathitem = pertinentnormsubpathitems[next]
589            thisarclen_pt = arclens_pt[this]
590            nextarclen_pt = arclens_pt[next]
591
592            # insert the middle segment
593            newnormsubpath.append(thisnormsubpathitem.segments(thisparams)[0])
594
595            # insert replacement curves for the corners
596            if next or normsubpath.closed:
597
598                t1 = thisnormsubpathitem.rotation([thisparams[1]])[0].apply_pt(1, 0)
599                t2 = nextnormsubpathitem.rotation([nextparams[0]])[0].apply_pt(1, 0)
600                # TODO: normpath.invalid
601
602                if (isinstance(thisnormsubpathitem, normpath.normline_pt) and
603                    isinstance(nextnormsubpathitem, normpath.normline_pt)):
604
605                    # case of two lines -> replace by two curves
606                    d1, g1, f1, e, f2, g2, d2 = curvescontrols_from_endlines_pt(
607                        thisnormsubpathitem.atend_pt(), t1, t2,
608                        thisarclen_pt*(1-thisparams[1]), nextarclen_pt*(nextparams[0]), softness=self.softness)
609
610                    p1 = thisnormsubpathitem.at_pt([thisparams[1]])[0]
611                    p2 = nextnormsubpathitem.at_pt([nextparams[0]])[0]
612
613                    newnormsubpath.append(normpath.normcurve_pt(*(d1 + g1 + f1 + e)))
614                    newnormsubpath.append(normpath.normcurve_pt(*(e + f2 + g2 + d2)))
615
616                else:
617
618                    # generic case -> replace by a single curve with prescribed tangents and curvatures
619                    p1 = thisnormsubpathitem.at_pt([thisparams[1]])[0]
620                    p2 = nextnormsubpathitem.at_pt([nextparams[0]])[0]
621                    c1 = thisnormsubpathitem.curvature_pt([thisparams[1]])[0]
622                    c2 = nextnormsubpathitem.curvature_pt([nextparams[0]])[0]
623                    # TODO: normpath.invalid
624
625                    # TODO: more intelligent fallbacks:
626                    #   circle -> line
627                    #   circle -> circle
628
629                    if not self.obeycurv:
630                        # do not obey the sign of the curvature but
631                        # make the sign such that the curve smoothly passes to the next point
632                        # this results in a discontinuous curvature
633                        # (but the absolute value is still continuous)
634                        s1 = +mathutils.sign(t1[0] * (p2[1]-p1[1]) - t1[1] * (p2[0]-p1[0]))
635                        s2 = -mathutils.sign(t2[0] * (p2[1]-p1[1]) - t2[1] * (p2[0]-p1[0]))
636                        c1 = s1 * abs(c1)
637                        c2 = s2 * abs(c2)
638
639                    # get the length of the control "arms"
640                    controldists = controldists_from_endgeometry_pt(p1, p2, t1, t2, c1, c2)
641
642                    if controldists and (controldists[0][0] >= 0 and controldists[0][1] >= 0):
643                        # use the first entry in the controldists
644                        # this should be the "smallest" pair
645                        a, d = controldists[0]
646                        # avoid curves with invalid parameterization
647                        a = max(a, epsilon)
648                        d = max(d, epsilon)
649
650                        # avoid overshooting at the corners:
651                        # this changes not only the sign of the curvature
652                        # but also the magnitude
653                        if not self.obeycurv:
654                            t, s = intersection(p1, p2, t1, t2)
655                            if (t is not None and s is not None and
656                                t > 0 and s < 0):
657                                a = min(a, abs(t))
658                                d = min(d, abs(s))
659
660                    else:
661                        # use a fallback
662                        t, s = intersection(p1, p2, t1, t2)
663                        if t is not None and s is not None:
664                            a = 0.65 * abs(t)
665                            d = 0.65 * abs(s)
666                        else:
667                            # if there is no useful result:
668                            # take an arbitrary smoothing curve that does not obey
669                            # the curvature constraints
670                            dist = math.hypot(p1[0] - p2[0], p1[1] - p2[1])
671                            a = dist / (3.0 * math.hypot(*t1))
672                            d = dist / (3.0 * math.hypot(*t2))
673
674                    # calculate the two missing control points
675                    q1 = p1[0] + a * t1[0], p1[1] + a * t1[1]
676                    q2 = p2[0] - d * t2[0], p2[1] - d * t2[1]
677
678                    newnormsubpath.append(normpath.normcurve_pt(*(p1 + q1 + q2 + p2)))
679
680        if normsubpath.closed:
681            newnormsubpath.close()
682        return newnormsubpath
683
684# >>>
685
686cornersmoothed.clear = attr.clearclass(cornersmoothed)
687smoothed = cornersmoothed
688smoothed.clear = attr.clearclass(smoothed)
689
690class mynormpathparam(normpath.normpathparam): # <<<
691    """In the parallel deformer we use properties such as the curvature, which
692    are not continuous on a path (at points between normsubpathitems). We
693    therefore require a better parameter class which really resolves the
694    nsp-item """
695
696    # TODO: find reasonable values for these eps:
697    rounding_eps = 1.0e-8
698
699    def __init__(self, np, normsubpathindex, normsubpathitemindex, normsubpathitemparam):
700        normpath.normpathparam.__init__(self, np, normsubpathindex, normsubpathitemindex + normsubpathitemparam)
701        self.normsubpathitemindex = normsubpathitemindex
702        self.normsubpathitemparam = normsubpathitemparam
703        self.beg_nspitem_known = False
704        self.end_nspitem_known = False
705
706        # guarantee that normpath.normpathparam always gets the correct item:
707        if int(self.normsubpathparam) != self.normsubpathitemindex:
708            if int(self.normsubpathparam) == self.normsubpathitemindex - 1:
709                self.normsubpathparam = self.normsubpathitemindex + self.rounding_eps
710                self.beg_nspitem_known = True
711            elif int(self.normsubpathparam) == self.normsubpathitemindex + 1:
712                self.normsubpathparam = (self.normsubpathitemindex + 1) - self.rounding_eps
713                self.end_nspitem_known = True
714            else:
715                assert False
716        assert int(self.normsubpathparam) == self.normsubpathitemindex
717        #assert 0 <= self.normsubpathparam - self.normsubpathitemindex
718        #assert 1 >= self.normsubpathparam - self.normsubpathitemindex
719
720    def __str__(self):
721        return "npp(%d, %d, %.3f)" % (self.normsubpathindex, self.normsubpathitemindex, self.normsubpathitemparam)
722
723    def __eq__(self, other):
724        if isinstance(other, mynormpathparam):
725            assert self.normpath is other.normpath, "normpathparams have to belong to the same normpath"
726            return (self.normsubpathindex, self.normsubpathitemindex, self.normsubpathitemparam) == (other.normsubpathindex, other.normsubpathitemindex, other.normsubpathitemparam)
727        else:
728            normpath.normpathparam.__eq__(self, other)
729
730    def __lt__(self, other):
731        if isinstance(other, mynormpathparam):
732            assert self.normpath is other.normpath, "normpathparams have to belong to the same normpath"
733            return (self.normsubpathindex, self.normsubpathitemindex, self.normsubpathitemparam) < (other.normsubpathindex, other.normsubpathitemindex, other.normsubpathitemparam)
734        else:
735            normpath.normpathparam.__eq__(self, other)
736
737    def __hash__(self):
738        return id(self)
739
740    def smaller_equiv(self, epsilon=None):
741        """Returns smaller equivalent parameter, if self is between two nsp-items"""
742        if not self.is_beg_of_nspitem(epsilon):
743            return self
744        nsp = self.normpath[self.normsubpathindex]
745        nspi_index = self.normsubpathitemindex - 1
746        nspi_param = 1
747        if nsp.closed:
748            nspi_index = nspi_index % len(nsp)
749        elif nspi_index < 0:
750            nspi_index = 0
751            nspi_param = 0
752        other = mynormpathparam(self.normpath, self.normsubpathindex, nspi_index, nspi_param)
753        if self.is_equiv(other, epsilon):
754            return other
755        return self
756
757    def larger_equiv(self, epsilon=None):
758        """Returns smaller equivalent parameter, if self is between two nsp-items"""
759        if not self.is_end_of_nspitem(epsilon):
760            return self
761        nsp = self.normpath[self.normsubpathindex]
762        nspi_index = self.normsubpathitemindex + 1
763        nspi_param = 0
764        if nsp.closed:
765            nspi_index = nspi_index % len(nsp)
766        elif nspi_index >= len(nsp):
767            nspi_index = len(nsp) - 1
768            nspi_param = 1
769        other = mynormpathparam(self.normpath, self.normsubpathindex, nspi_index, nspi_param)
770        if self.is_equiv(other, epsilon):
771            return other
772        return self
773
774    def is_equiv(self, other, epsilon=None):
775        """Test whether the two params yield essentially the same point"""
776        assert self.normpath is other.normpath, "normpathparams have to belong to the same normpath"
777        if self.normsubpathindex != other.normsubpathindex:
778            return False
779        nsp = self.normpath[self.normsubpathindex]
780        if epsilon is None:
781            epsilon = nsp.epsilon
782        A, B = self.normpath.at_pt([self, other])
783        return math.hypot(A[0]-B[0], A[1]-B[1]) < epsilon
784
785    def is_beg_of_nspitem(self, epsilon=None):
786        if self.beg_nspitem_known:
787            return True
788        return self.is_equiv(mynormpathparam(self.normpath, self.normsubpathindex, self.normsubpathitemindex, 0), epsilon)
789
790    def is_end_of_nspitem(self, epsilon=None):
791        if self.end_nspitem_known:
792            return True
793        return self.is_equiv(mynormpathparam(self.normpath, self.normsubpathindex, self.normsubpathitemindex, 1), epsilon)
794
795    def is_beg_of_nsp(self, epsilon=None):
796        if self.normsubpathitemindex > 0:
797            return False
798        return self.is_equiv(mynormpathparam(self.normpath, self.normsubpathindex, 0, 0), epsilon)
799
800    def is_end_of_nsp(self, epsilon=None):
801        n = len(self.normpath[self.normsubpathindex]) - 1
802        if self.normsubpathitemindex < n:
803            return False
804        return self.is_equiv(mynormpathparam(self.normpath, self.normsubpathindex, n, 1), epsilon)
805
806# >>>
807def _length_pt(path, param1, param2): # <<<
808    point1, point2 = path.at_pt([param1, param2])
809    return math.hypot(point1[0] - point2[0], point1[1] - point2[1])
810# >>>
811class parallel(baseclasses.deformer): # <<<
812
813    """creates a parallel normpath with constant distance to the original normpath
814
815    A positive 'distance' results in a curve left of the original one -- and a
816    negative 'distance' in a curve at the right. Left/right are understood in
817    terms of the parameterization of the original curve. The construction of
818    the paralel path is done in two steps: First, an "extended" parallel path
819    is built. For each path element a parallel curve/line is constructed, which
820    can be too long or too short, depending on the presence of corners. At
821    corners, either a circular arc is drawn around the corner, or, if possible,
822    the parallel curve is cut in order to also exhibit a corner. In a second
823    step all self-intersection points are determined and unnecessary parts of
824    the path are cut away.
825
826    distance:            the distance of the parallel normpath
827    relerr:              distance*relerr is the maximal allowed error in the parallel distance
828    sharpoutercorners:   make the outer corners not round but sharp.
829                         The inner corners (corners after inflection points) will stay round
830    dointersection:      boolean for doing the intersection step (default: 1).
831                         Set this value to 0 if you want the whole parallel path
832    checkdistanceparams: a list of parameter values in the interval (0,1) where the
833                         parallel distance is checked on each normpathitem
834    lookforcurvatures:   number of points per normpathitem where is looked for
835                         a critical value of the curvature
836    """
837
838    # TODO:
839    # * DECIDE MEANING of arcs around corners (see case L in test/functional/test_deformer.py)
840    # * eliminate double, triple, ... pairs
841    # * implement self-intersection of normcurve_pt
842    # * implement _between_paths also for normcurve_pt
843
844
845    def __init__(self, distance, relerr=0.05, sharpoutercorners=False, dointersection=True,
846                       checkdistanceparams=[0.5], lookforcurvatures=11, searchstep=0.01, debug=None):
847        self.distance = distance
848        self.relerr = relerr
849        self.sharpoutercorners = sharpoutercorners
850        self.checkdistanceparams = checkdistanceparams
851        self.lookforcurvatures = lookforcurvatures
852        self.dointersection = dointersection
853        self.searchstep = searchstep
854        self.debug = debug
855
856    def __call__(self, distance=None, relerr=None, sharpoutercorners=None, dointersection=None,
857                       checkdistanceparams=None, lookforcurvatures=None, searchstep=None, debug=None):
858        # returns a copy of the deformer with different parameters
859        if distance is None:
860            distance = self.distance
861        if relerr is None:
862            relerr = self.relerr
863        if sharpoutercorners is None:
864            sharpoutercorners = self.sharpoutercorners
865        if dointersection is None:
866            dointersection = self.dointersection
867        if checkdistanceparams is None:
868            checkdistanceparams = self.checkdistanceparams
869        if lookforcurvatures is None:
870            lookforcurvatures = self.lookforcurvatures
871        if searchstep is None:
872            searchstep = self.searchstep
873        if debug is None:
874            debug = self.debug
875
876        return parallel(distance=distance, relerr=relerr,
877                        sharpoutercorners=sharpoutercorners,
878                        dointersection=dointersection,
879                        checkdistanceparams=checkdistanceparams,
880                        lookforcurvatures=lookforcurvatures,
881                        searchstep=searchstep,
882                        debug=debug)
883
884    def deform(self, basepath):
885        basepath = basepath.normpath()
886        self.dist_pt = unit.topt(self.distance)
887        resultnormsubpaths = []
888        par_to_orig = {}
889        for nsp in basepath.normsubpaths:
890            parallel_normpath, tmp1, tmp2, par2orig = self.deformsubpath(nsp)
891            resultnormsubpaths += parallel_normpath.normsubpaths
892            for key in par2orig:
893                par_to_orig[key] = par2orig[key]
894        result = normpath.normpath(resultnormsubpaths)
895
896        if self.dointersection:
897            result = self.rebuild_intersected_normpath(result, basepath, par_to_orig)
898
899        return result
900
901    def deformsubpath(self, orig_nsp): # <<<
902
903        """Performs the first step of the deformation: Returns a list of
904        normsubpaths building the parallel to the given normsubpath.
905        Then calls the intersection routine to do the second step."""
906        # the default case is not to join the result.
907
908        dist = self.dist_pt
909        epsilon = orig_nsp.epsilon
910
911        if len(orig_nsp.normsubpathitems) == 0:
912            return normpath.normpath([]), None, None, {}
913
914        # avoid too small dists: we would run into instabilities
915        if abs(dist) < abs(epsilon):
916            par_to_orig = {}
917            for nspitem in orig_nsp:
918                par_to_orig[nspitem] = nspitem
919            return normpath.normpath([orig_nsp]), None, None, par_to_orig
920
921        result = None
922        par_to_orig = None
923        join_begin = None
924        prev_joinend = None
925
926        # iterate over the normsubpath in the following way:
927        # * for each item first append the additional arc
928        #   and then add the next parallel piece
929        # * for the first item only add the parallel piece
930        #   (because this is done for curr_orig_nspitem, we need to start with i=0)
931        for i in range(len(orig_nsp.normsubpathitems)):
932            prev_orig_nspitem = orig_nsp.normsubpathitems[i-1]
933            curr_orig_nspitem = orig_nsp.normsubpathitems[i]
934
935            # get the next parallel piece for the normpath
936            next_parallel_normpath, joinbeg, joinend, par2orig = self.deformsubpathitem(curr_orig_nspitem, epsilon)
937            if result is None:
938                if join_begin is None:
939                    join_begin = joinbeg
940                else:
941                    join_begin = (join_begin and joinbeg)
942
943            if not (next_parallel_normpath.normsubpaths and next_parallel_normpath[0].normsubpathitems):
944                if prev_joinend is None:
945                    prev_joinend = joinend
946                else:
947                    prev_joinend = (prev_joinend and joinend)
948                continue
949
950            # this starts the whole normpath
951            if result is None:
952                result = next_parallel_normpath
953                par_to_orig = {}
954                for key in par2orig:
955                    par_to_orig[key] = par2orig[key]
956                prev_joinend = joinend
957                continue # there is nothing to join
958
959            prev_tangent, next_tangent, is_straight, is_concave = self._get_angles(prev_orig_nspitem, curr_orig_nspitem, epsilon)
960            if not (joinbeg and prev_joinend): # split due to loo large curvature
961                result += next_parallel_normpath
962            elif is_straight:
963                # The path is quite straight between prev and next item:
964                # normpath.normpath.join adds a straight line if necessary
965                result.join(next_parallel_normpath)
966            else:
967                # The parallel path can be extended continuously.
968                # We must add a corner or an arc around the corner:
969                cornerpath = self._path_around_corner(curr_orig_nspitem.atbegin_pt(), result.atend_pt(), next_parallel_normpath.atbegin_pt(),
970                                                      prev_tangent, next_tangent, is_concave, epsilon)
971                result.join(cornerpath)
972                assert len(cornerpath) == 1
973                corner = curr_orig_nspitem.atbegin_pt()
974                for cp_nspitem in cornerpath[0]:
975                    par_to_orig[cp_nspitem] = corner
976                # append the next parallel piece to the path
977                result.join(next_parallel_normpath)
978            for key in par2orig:
979                par_to_orig[key] = par2orig[key]
980            prev_joinend = joinend
981
982        # end here if nothing has been found so far
983        if result is None:
984            return normpath.normpath(), False, False, {}
985
986        # the curve around the closing corner may still be missing
987        if orig_nsp.closed:
988            prev_tangent, next_tangent, is_straight, is_concave = self._get_angles(orig_nsp.normsubpathitems[-1], orig_nsp.normsubpathitems[0], epsilon)
989            if not (joinend and join_begin): # do not close because of infinite curvature
990                do_close = False
991            elif is_straight:
992                # The path is quite straight at end and beginning
993                do_close = True
994            else:
995                # The parallel path can be extended continuously.
996                do_close = True
997                # We must add a corner or an arc around the corner:
998                cornerpath = self._path_around_corner(orig_nsp.atend_pt(), result.atend_pt(), result.atbegin_pt(),
999                                                      prev_tangent, next_tangent, is_concave, epsilon)
1000                result.join(cornerpath)
1001                corner = orig_nsp.atend_pt()
1002                assert len(cornerpath) == 1
1003                for cp_nspitem in cornerpath[0]:
1004                    par_to_orig[cp_nspitem] = corner
1005
1006            if do_close:
1007                if len(result) == 1:
1008                    result[0].close()
1009                else:
1010                    # if the parallel normpath is split into several subpaths anyway,
1011                    # then use the natural beginning and ending
1012                    # closing is not possible anymore
1013                    for nspitem in result[0]:
1014                        result[-1].append(nspitem)
1015                    result.normsubpaths = result.normsubpaths[1:]
1016                join_begin, joinend = False, False
1017
1018        return result, join_begin, joinend, par_to_orig
1019        # >>>
1020    def deformsubpathitem(self, nspitem, epsilon): # <<<
1021
1022        """Returns a parallel normpath for a single normsubpathitem
1023
1024        Analyzes the curvature of a normsubpathitem and returns a normpath with
1025        the appropriate number of normsubpaths. This must be a normpath because
1026        a normcurve can be strongly curved, such that the parallel path must
1027        contain a hole"""
1028        # the default case is to join the result. Only if there was an infinite
1029        # curvature at beginning/end, we return info not to join it.
1030
1031        dist = self.dist_pt
1032
1033        # for a simple line we return immediately
1034        if isinstance(nspitem, normpath.normline_pt):
1035            normal = nspitem.rotation([0])[0].apply_pt(0, 1)
1036            start = nspitem.atbegin_pt()
1037            end = nspitem.atend_pt()
1038            result = path.line_pt(start[0] + dist * normal[0], start[1] + dist * normal[1],
1039                                  end[0] + dist * normal[0], end[1] + dist * normal[1]).normpath(epsilon=epsilon)
1040            assert len(result) == 1 and len(result[0]) == 1
1041            return result, True, True, {result[0][0]:nspitem}
1042
1043        # for a curve we have to check if the curvatures
1044        # cross the singular value 1/dist
1045        crossings = list(self._distcrossingparameters(nspitem, epsilon))
1046        crossings.sort()
1047
1048        # depending on the number of crossings we must consider
1049        # three different cases:
1050        if crossings:
1051            # The curvature crosses the borderline 1/dist
1052            # the parallel curve contains points with infinite curvature!
1053            parallcurvs = [inf_curvature]*len(crossings)
1054
1055            result = normpath.normpath()
1056            join_begin, join_end = False, False
1057            par_to_orig = {}
1058
1059            # we need the endpoints of the nspitem
1060            if _length_pt(nspitem, crossings[0], 0) > epsilon:
1061                crossings.insert(0, 0)
1062                parallcurvs.insert(0, None)
1063            if _length_pt(nspitem, crossings[-1], 1) > epsilon:
1064                crossings.append(1)
1065                parallcurvs.append(None)
1066
1067            for i in range(len(crossings) - 1):
1068                middleparam = 0.5*(crossings[i] + crossings[i+1])
1069                middlecurv = nspitem.curvature_pt([middleparam])[0]
1070                # the radius is good if
1071                #  - middlecurv and dist have opposite signs : distance vector points "out" of the original curve
1072                #  - middlecurv is "smaller" than 1/dist : original curve is less curved than +-1/dist
1073                if dist*middlecurv < 0 or abs(dist*middlecurv) < 1:
1074                    if i == 0:
1075                        join_begin = True
1076                    elif i == len(crossings) - 2:
1077                        join_end = True
1078                    parallel_nsp, par2orig = self.deformnicecurve(nspitem.segments(crossings[i:i+2])[0], epsilon, curvA=parallcurvs[i], curvD=parallcurvs[i+1])
1079                    # never append empty normsubpaths
1080                    if parallel_nsp.normsubpathitems:
1081                        # infinite curvatures interrupt the path and start a new nsp
1082                        result.append(parallel_nsp)
1083                        for key in par2orig:
1084                            par_to_orig[key] = par2orig[key]
1085            if not (result.normsubpaths and result[0].normsubpathitems):
1086                return normpath.normpath(), True, True, {}
1087            return result, join_begin, join_end, par_to_orig
1088
1089        # the curvature is either bigger or smaller than 1/dist
1090        middlecurv = nspitem.curvature_pt([0.5])[0]
1091        if dist*middlecurv < 0 or abs(dist*middlecurv) < 1:
1092            # The curve is everywhere less curved than 1/dist
1093            # We can proceed finding the parallel curve for the whole piece
1094            parallel_nsp, par2orig = self.deformnicecurve(nspitem, epsilon)
1095            # never append empty normsubpaths
1096            if parallel_nsp.normsubpathitems:
1097                par_to_orig = {}
1098                for key in par2orig:
1099                    par_to_orig[key] = par2orig[key]
1100                return normpath.normpath([parallel_nsp]), True, True, par_to_orig
1101
1102        # the curve is everywhere stronger curved than 1/dist
1103        # There is nothing to be returned.
1104        return normpath.normpath(), False, False, {}
1105        # >>>
1106    def deformnicecurve(self, normcurve, epsilon, startparam=0.0, endparam=1.0, curvA=None, curvD=None): # <<<
1107        """Returns a parallel normsubpath for the normcurve.
1108
1109        This routine assumes that the normcurve is everywhere
1110        'less' curved than 1/dist. Only at the ends, the curvature
1111        can be exactly +-1/dist, which is marked by curvA and/or curvD.
1112        """
1113        dist = self.dist_pt
1114
1115        # normalized tangent directions
1116        tangA, tangD = normcurve.rotation([startparam, endparam])
1117        tangA = tangA.apply_pt(1, 0)
1118        tangD = tangD.apply_pt(1, 0)
1119
1120        # the new starting points
1121        orig_A, orig_D = normcurve.at_pt([startparam, endparam])
1122        A = orig_A[0] - dist * tangA[1], orig_A[1] + dist * tangA[0]
1123        D = orig_D[0] - dist * tangD[1], orig_D[1] + dist * tangD[0]
1124
1125        # we need to end this _before_ we will run into epsilon-problems
1126        # when creating curves we do not want to calculate the length of
1127        # or even split it for recursive calls
1128        if (math.hypot(A[0] - D[0], A[1] - D[1]) < epsilon and
1129            abs(dist)*(tangA[0]*tangD[1] - tangA[1]*tangD[0]) < epsilon):
1130            nl = normpath.normline_pt(A[0], A[1], D[0], D[1])
1131            return normpath.normsubpath([nl]), {nl:normcurve}
1132
1133        result = normpath.normsubpath(epsilon=epsilon)
1134        # is there enough space on the normals before they intersect?
1135        a, d = intersection(orig_A, orig_D, (-tangA[1], tangA[0]), (-tangD[1], tangD[0]))
1136        # a,d are the lengths to the intersection points:
1137        # for a (and equally for b) we can proceed in one of the following cases:
1138        #   a is None (means parallel normals)
1139        #   a and dist have opposite signs (and the same for b)
1140        #   a has the same sign but is bigger
1141        if ( (a is None or a*dist < 0 or abs(a) > abs(dist) + epsilon) or
1142             (d is None or d*dist < 0 or abs(d) > abs(dist) + epsilon) ):
1143            # the original path is long enough to draw a parallel piece
1144            # this is the generic case. Get the parallel curves
1145            orig_curvA, orig_curvD = normcurve.curvature_pt([startparam, endparam])
1146            if curvA is None:
1147                curvA = orig_curvA / (1.0 - dist*orig_curvA)
1148            if curvD is None:
1149                curvD = orig_curvD / (1.0 - dist*orig_curvD)
1150
1151            # first try to approximate the normcurve with a single item
1152            controldistpairs = controldists_from_endgeometry_pt(A, D, tangA, tangD, curvA, curvD)
1153
1154            if controldistpairs:
1155                # TODO: is it good enough to get the first entry here?
1156                #       from testing: this fails if there are loops in the original curve
1157                a, d = controldistpairs[0]
1158                if a >= 0 and d >= 0:
1159                    # we avoid to create curves with invalid parameterization
1160                    if a < epsilon and d < epsilon:
1161                        result = normpath.normsubpath([normpath.normline_pt(A[0], A[1], D[0], D[1])], epsilon=epsilon)
1162                    else:
1163                        a = max(a, epsilon)
1164                        d = max(d, epsilon)
1165                        result = normpath.normsubpath([normpath.normcurve_pt(
1166                            A[0], A[1],
1167                            A[0] + a * tangA[0], A[1] + a * tangA[1],
1168                            D[0] - d * tangD[0], D[1] - d * tangD[1],
1169                            D[0], D[1])], epsilon=epsilon)
1170
1171            # then try with two items, recursive call
1172            if ((not result.normsubpathitems) or
1173                (self.checkdistanceparams and result.normsubpathitems
1174                 and not self._distchecked(normcurve, result, epsilon, startparam, endparam))):
1175                # TODO: does this ever converge?
1176                # TODO: what if this hits epsilon?
1177                middleparam = 0.5*(startparam + endparam)
1178                firstnsp, first_par2orig = self.deformnicecurve(normcurve, epsilon, startparam, middleparam, curvA, None)
1179                secondnsp, second_par2orig = self.deformnicecurve(normcurve, epsilon, middleparam, endparam, None, curvD)
1180                if not (firstnsp.normsubpathitems and secondnsp.normsubpathitems):
1181                    result = normpath.normsubpath(
1182                        [normpath.normline_pt(A[0], A[1], D[0], D[1])], epsilon=epsilon)
1183                else:
1184                    result = firstnsp.joined(secondnsp)
1185
1186        par_to_orig = {}
1187        for key in result:
1188            par_to_orig[key] = normcurve
1189        return result, par_to_orig
1190        # >>>
1191
1192    def _path_around_corner(self, corner_pt, beg_pt, end_pt, beg_tangent, end_tangent, is_concave, epsilon): # <<<
1193        """Helper routine for parallel.deformsubpath: Draws an arc around a convex corner"""
1194        if self.sharpoutercorners and not is_concave:
1195            # straight lines:
1196            t1, t2 = intersection(beg_pt, end_pt, beg_tangent, end_tangent)
1197            B = beg_pt[0] + t1 * beg_tangent[0], beg_pt[1] + t1 * beg_tangent[1]
1198            return normpath.normpath([normpath.normsubpath([
1199                normpath.normline_pt(beg_pt[0], beg_pt[1], B[0], B[1]),
1200                normpath.normline_pt(B[0], B[1], end_pt[0], end_pt[1])
1201                ])])
1202
1203        # We append an arc around the corner
1204        # these asserts fail in test case "E"
1205        #assert abs(math.hypot(beg_pt[1] - corner_pt[1], beg_pt[0] - corner_pt[0]) - abs(self.dist_pt)) < epsilon
1206        #assert abs(math.hypot(end_pt[1] - corner_pt[1], end_pt[0] - corner_pt[0]) - abs(self.dist_pt)) < epsilon
1207        angle1 = math.atan2(beg_pt[1] - corner_pt[1], beg_pt[0] - corner_pt[0])
1208        angle2 = math.atan2(end_pt[1] - corner_pt[1], end_pt[0] - corner_pt[0])
1209
1210        # depending on the direction we have to use arc or arcn
1211        sinangle = beg_tangent[0]*end_tangent[1] - beg_tangent[1]*end_tangent[0] # >0 for left-turning, <0 for right-turning
1212        if self.dist_pt > 0:
1213            arcclass = path.arcn_pt
1214        else:
1215            arcclass = path.arc_pt
1216        return path.path(arcclass(
1217          corner_pt[0], corner_pt[1], abs(self.dist_pt),
1218          math.degrees(angle1), math.degrees(angle2))).normpath(epsilon=epsilon)
1219    # >>>
1220    def _distchecked(self, orig_normcurve, parallel_normsubpath, epsilon, tstart, tend): # <<<
1221        """Helper routine for parallel.deformnicecurve: Checks the distances between orig_normcurve and parallel_normsubpath.
1222
1223        The checking is done at parameters self.checkdistanceparams of orig_normcurve."""
1224
1225        dist = self.dist_pt
1226        # do not look closer than epsilon:
1227        dist_err = mathutils.sign(dist) * max(abs(self.relerr*dist), epsilon)
1228
1229        checkdistanceparams = [tstart + (tend-tstart)*t for t in self.checkdistanceparams]
1230
1231        for param, P, rotation in zip(checkdistanceparams,
1232                                      orig_normcurve.at_pt(checkdistanceparams),
1233                                      orig_normcurve.rotation(checkdistanceparams)):
1234            normal = rotation.apply_pt(0, 1)
1235
1236            # create a short cutline for intersection only:
1237            cutline = normpath.normsubpath([normpath.normline_pt(
1238              P[0] + (dist - 2*dist_err) * normal[0], P[1] + (dist - 2*dist_err) * normal[1],
1239              P[0] + (dist + 2*dist_err) * normal[0], P[1] + (dist + 2*dist_err) * normal[1])], epsilon=epsilon)
1240
1241            cutparams = parallel_normsubpath.intersect(cutline)
1242            distances = [math.hypot(P[0] - cutpoint[0], P[1] - cutpoint[1])
1243                         for cutpoint in cutline.at_pt(cutparams[1])]
1244
1245            if (not distances) or (abs(min(distances) - abs(dist)) > abs(dist_err)):
1246                return False
1247
1248        return True
1249    # >>>
1250    def _distcrossingparameters(self, normcurve, epsilon, tstart=0, tend=1): # <<<
1251        """Helper routine for parallel.deformsubpathitem: Returns a list of parameters where the curvature of normcurve is 1/distance"""
1252
1253        assert tstart < tend
1254        dist = self.dist_pt
1255
1256        # we _need_ to do this with the curvature, not with the radius
1257        # because the curvature is continuous at the straight line and the radius is not:
1258        # when passing from one slightly curved curve to the other with opposite curvature sign,
1259        # via the straight line, then the curvature changes its sign at curv=0, while the
1260        # radius changes its sign at +/-infinity
1261        # this causes instabilities for nearly straight curves
1262
1263        # include tstart and tend
1264        params = [tstart + i * (tend - tstart) / (self.lookforcurvatures - 1.0)
1265                  for i in range(self.lookforcurvatures)]
1266        curvs = normcurve.curvature_pt(params)
1267
1268        parampairs = list(zip(params[:-1], params[1:]))
1269        curvpairs = list(zip(curvs[:-1], curvs[1:]))
1270
1271        crossingparams = set()
1272        for parampair, curvpair in zip(parampairs, curvpairs):
1273            begparam, endparam = parampair
1274            begcurv, endcurv = curvpair
1275            begchange = begcurv*dist - 1
1276            endchange = endcurv*dist - 1
1277            if begchange*endchange < 0:
1278                # the curvature crosses the value 1/dist
1279                # get the parmeter value by linear interpolation:
1280                middleparam = (
1281                  (begparam * abs(begchange) + endparam * abs(endchange)) /
1282                  (abs(begchange) + abs(endchange)))
1283                try:
1284                    middleradius = 1/normcurve.curvature_pt([middleparam])[0]
1285                except ArithmeticError:
1286                    raise InvalidParamException(middleparam)
1287
1288                if abs(middleradius - dist) < epsilon or endparam-begparam < 1.0e-14:
1289                    # get the parmeter value by linear interpolation:
1290                    crossingparams.add(middleparam)
1291                else:
1292                    # call recursively:
1293                    for x in self._distcrossingparameters(normcurve, epsilon, tstart=begparam, tend=endparam):
1294                        crossingparams.add(x)
1295            else:
1296                if begchange == 0:
1297                    crossingparams.add(begparam)
1298                if endchange == 0:
1299                    crossingparams.add(endparam)
1300
1301        return crossingparams
1302        # >>>
1303    def _get_angles(self, prev_nspitem, next_nspitem, epsilon): # <<<
1304        prev_rotation = prev_nspitem.rotation([1])[0]
1305        next_rotation = next_nspitem.rotation([0])[0]
1306        prev_tangent = prev_rotation.apply_pt(1, 0)
1307        prev_orthogo = prev_rotation.apply_pt(0, self.dist_pt) # points towards parallel path (prev_nspitem is on original path)
1308        next_tangent = next_rotation.apply_pt(1, 0)
1309        #sinangle = prev_tangent[0]*next_tangent[1] - prev_tangent[1]*next_tangent[0] # >0 for left-turning, <0 for right-turning
1310        cosangle = prev_tangent[0]*next_tangent[0] + prev_tangent[1]*next_tangent[1]
1311        proj = prev_orthogo[0]*next_tangent[0] + prev_orthogo[1]*next_tangent[1]
1312        is_straight = (cosangle > 0 and abs(proj) < epsilon)
1313        is_concave = (proj > 0)
1314        return prev_tangent, next_tangent, is_straight, is_concave
1315    # >>>
1316
1317    def rebuild_intersected_normpath(self, par_np, orig_np, par2orig, epsilon=None): # <<<
1318
1319        dist = self.dist_pt
1320        if epsilon is None:
1321            epsilon = orig_np.normsubpaths[0].epsilon
1322        eps_comparepairs = 10*epsilon
1323
1324        # calculate the self-intersections of the par_np
1325        forwardpairs, backwardpairs = self.normpath_selfintersections(par_np, epsilon, eps_comparepairs)
1326        # calculate the intersections of the par_np with the original path
1327        origintparams, orig_origintparams = self.normpath_origintersections(orig_np, par_np, epsilon)
1328        if not forwardpairs:
1329            if origintparams:
1330                return normpath.normpath()
1331            else:
1332                return par_np
1333
1334        # parameters at begin and end of subnormpaths:
1335        # omit those which start/end on the original path
1336        beginparams = []
1337        endparams = []
1338        testparams = origintparams + list(forwardpairs.keys()) + list(forwardpairs.values())
1339        for i, nsp in enumerate(par_np):
1340            beginparam = mynormpathparam(par_np, i, 0, 0)
1341            is_new = True
1342            for param in testparams:
1343                if beginparam.is_equiv(param):
1344                    is_new = False
1345                    break
1346            if is_new:
1347                beginparams.append(beginparam)
1348
1349            endparam = mynormpathparam(par_np, i, len(nsp)-1, 1)
1350            is_new = True
1351            for param in testparams:
1352                if endparam.is_equiv(param):
1353                    is_new = False
1354                    break
1355            if is_new:
1356                endparams.append(endparam)
1357        beginparams.sort()
1358        endparams.sort()
1359
1360        # we need a way to get the "next" param on the normpath
1361        # XXX why + beginparams + endparams ?
1362        allparams = list(forwardpairs.keys()) + list(backwardpairs.keys()) + origintparams + beginparams + endparams
1363        allparams.sort()
1364        done = {}
1365        for param in allparams:
1366            done[param] = False
1367        nextp = {}
1368        for i, param in enumerate(allparams[:-1]):
1369            nextp[param] = allparams[i+1]
1370        for endparam in endparams:
1371            if par_np[endparam.normsubpathindex].closed:
1372                begparam = [p for p in allparams if p.normsubpathindex == endparam.normsubpathindex][0]
1373                assert begparam.normsubpathitemindex == 0
1374                assert begparam.normsubpathitemparam == 0
1375                nextp[endparam] = begparam
1376            else:
1377                nextp[endparam] = None
1378
1379        # exclude all intersections that are between the original and the parallel path:
1380        # See for example test/functional/test_deformer (parallel Z): There can
1381        # remain a little piece of the path (triangle) that lies between a lot
1382        # of intersection points. Simple intersection rules such as thoe in
1383        # trial_parampairs cannot exclude this piece.
1384        for param in forwardpairs:
1385            if done[param] or done[forwardpairs[param]]:
1386                done[param] = done[forwardpairs[param]] = True
1387            elif self._between_paths(par_np.at_pt(param), par2orig, 4*epsilon):
1388                done[param] = done[forwardpairs[param]] = True
1389        for param in beginparams + endparams:
1390            if self._between_paths(par_np.at_pt(param), par2orig, 4*epsilon):
1391                done[param] = True
1392
1393        # visualize the intersection points: # <<<
1394        if self.debug is not None:
1395            for param1, param2 in forwardpairs.items():
1396                point1, point2 = par_np.at([param1, param2])
1397                if not done[param1]:
1398                    self.debug.fill(path.circle(point1[0], point1[1], 0.05), [color.rgb.red])
1399                if not done[param2]:
1400                    self.debug.fill(path.circle(point2[0], point2[1], 0.03), [color.rgb.black])
1401            for param in origintparams:
1402                #assert done[param]
1403                point = par_np.at([param])[0]
1404                self.debug.fill(path.circle(point[0], point[1], 0.05), [color.rgb.green])
1405            for i, nsp in enumerate(par_np):
1406              for j, nspi in enumerate(nsp):
1407                x, y = nspi.at_pt([0.5])[0]
1408                self.debug.text_pt(x, y, "{}/{}".format(i,j))#, [text.halign.center, text.vshift.mathaxis])
1409            print("aborted path intersection due to debug")
1410            return par_np
1411        # >>>
1412
1413        def ptype(param): # <<<
1414            if param in forwardpairs : return "fw with partner %s" % (forwardpairs[param])
1415            if param in backwardpairs : return "bw with partner %s" % (backwardpairs[param])
1416            if param in origintparams: return "orig"
1417            if param in beginparams: return "begin"
1418            if param in endparams: return "end"
1419        # >>>
1420        def trial_parampairs(startp): # <<<
1421            """Starting at startp, try to find a valid series of intersection parameters"""
1422            tried = {} # a local copy of done
1423            for param in allparams:
1424                tried[param] = done[param]
1425
1426            previousp = startp
1427            currentp = nextp[previousp]
1428            result = []
1429
1430            while True:
1431                # successful and unsuccessful termination conditions:
1432                if tried[currentp]:
1433                    # we reached a branch that has already been treated
1434                    # ==> this is not a valid parallel path
1435                    return []
1436                if currentp in origintparams:
1437                    # we cross the original path
1438                    # ==> this is not a valid parallel path
1439                    return []
1440                if currentp in backwardpairs:
1441                    # we reached a branch that should be followed from another part
1442                    # ==> this is not a valid parallel path
1443                    return []
1444                if currentp is startp:
1445                    # we have reached again the starting point on a closed subpath.
1446                    assert startp in beginparams
1447                    assert previousp in endparams
1448                    return result
1449                if currentp in forwardpairs:
1450                    result.append((previousp, currentp))
1451                    if forwardpairs[currentp] is startp:
1452                        # we have found the same point as the startp (its pair partner)
1453                        return result
1454                    previousp = forwardpairs[currentp]
1455                if currentp in endparams:
1456                    result.append((previousp, currentp))
1457                    if nextp[currentp] is None: # open subpath
1458                        # we have found the end of a non-closed subpath
1459                        return result
1460                    previousp = currentp # closed subpath
1461                # follow the crossings on valid startpairs
1462                tried[currentp] = True
1463                tried[previousp] = True
1464                currentp = nextp[previousp]
1465            assert False # never reach this point
1466        # >>>
1467
1468        # first the paths that start at the beginning of a subnormpath:
1469        result = normpath.normpath()
1470        # paths can start on subnormpaths or crossings where we can "get away":
1471        bwkeys = list(backwardpairs.keys())
1472        bwkeys.sort()
1473        for startp in beginparams + bwkeys:
1474            if done[startp]:
1475                continue
1476
1477            # try to find a valid series of intersection points:
1478            parampairs = trial_parampairs(startp)
1479            if not parampairs:
1480                continue
1481
1482            # collect all the pieces between parampairs:
1483            add_nsp = normpath.normsubpath(epsilon=epsilon)
1484            for begin, end in parampairs:
1485                # check that trial_parampairs works correctly
1486                assert begin is not end
1487                for item in par_np[begin.normsubpathindex].segments(
1488                    [begin.normsubpathparam, end.normsubpathparam])[0].normsubpathitems:
1489                    # TODO: this should be obsolete with an improved intersection algorithm
1490                    #       guaranteeing epsilon
1491                    if add_nsp.normsubpathitems:
1492                        item = item.modifiedbegin_pt(*(add_nsp.atend_pt()))
1493                    add_nsp.append(item)
1494
1495                done[begin] = True
1496                done[end] = True
1497
1498            # close the path if necessary
1499            if add_nsp:
1500                if ((parampairs[-1][-1] in forwardpairs and forwardpairs[parampairs[-1][-1]] is parampairs[0][0]) or
1501                    (parampairs[-1][-1] in endparams and parampairs[0][0] in beginparams and parampairs[0][0] is nextp[parampairs[-1][-1]])):
1502                    add_nsp.normsubpathitems[-1] = add_nsp.normsubpathitems[-1].modifiedend_pt(*add_nsp.atbegin_pt())
1503                    add_nsp.close()
1504
1505            result.extend([add_nsp])
1506
1507        return result
1508    # >>>
1509    def normpath_selfintersections(self, np, epsilon, eps_comparepairs): # <<<
1510
1511        """Returns all self-intersection points of normpath np.
1512
1513        This does not include the intersections of a single normcurve with itself,
1514        but all intersections of one normpathitem with a different one in the path.
1515        The intersection pairs are such that the parallel path can be continued
1516        from the first to the second parameter, but not vice-versa."""
1517
1518        dist = self.dist_pt
1519
1520        n = len(np)
1521        forwardpairs = {}
1522        for nsp_i in range(n):
1523            for nsp_j in range(nsp_i, n):
1524                for nspitem_i in range(len(np[nsp_i])):
1525                    if nsp_j == nsp_i:
1526                        nspitem_j_range = list(range(nspitem_i+1, len(np[nsp_j])))
1527                    else:
1528                        nspitem_j_range = list(range(len(np[nsp_j])))
1529                    for nspitem_j in nspitem_j_range:
1530                        intsparams = np[nsp_i][nspitem_i].intersect(np[nsp_j][nspitem_j], epsilon)
1531                        if intsparams:
1532                            for intsparam_i, intsparam_j in intsparams:
1533                                npp_i = mynormpathparam(np, nsp_i, nspitem_i, intsparam_i)
1534                                npp_j = mynormpathparam(np, nsp_j, nspitem_j, intsparam_j)
1535
1536                                # skip successive nsp-items
1537                                if nsp_i == nsp_j:
1538                                    if nspitem_j == nspitem_i+1 and (npp_i.is_end_of_nspitem(epsilon) or npp_j.is_beg_of_nspitem(epsilon)):
1539                                        continue
1540                                    if np[nsp_i].closed and ((npp_i.is_beg_of_nsp(epsilon) and npp_j.is_end_of_nsp(epsilon)) or
1541                                                             (npp_j.is_beg_of_nsp(epsilon) and npp_i.is_end_of_nsp(epsilon))):
1542                                        continue
1543
1544                                # correct the order of the pair, such that we can use it to continue on the path
1545                                if not self._can_continue(npp_i, npp_j, epsilon):
1546                                    assert self._can_continue(npp_j, npp_i, epsilon)
1547                                    npp_i, npp_j = npp_j, npp_i
1548
1549                                # if the intersection is between two nsp-items, take the smallest -> largest
1550                                npp_i = npp_i.smaller_equiv(5*epsilon)
1551                                npp_j = npp_j.larger_equiv(5*epsilon)
1552
1553                                # because of the above change of npp_ij, and because there may be intersections between nsp-items,
1554                                # it may happen that we try to insert two times the same pair
1555                                if self._skip_intersection_doublet(npp_i, npp_j, forwardpairs, eps_comparepairs):
1556                                    continue
1557                                forwardpairs[npp_i] = npp_j
1558
1559        # this is partially done in _skip_intersection_doublet
1560        #forwardpairs = self._elim_intersection_doublets(forwardpairs, eps_comparepairs)
1561        # create the reverse mapping
1562        backwardpairs = {}
1563        for p, q in forwardpairs.items():
1564            backwardpairs[q] = p
1565        return forwardpairs, backwardpairs
1566
1567    # >>>
1568    def normpath_origintersections(self, orig_np, par_np, epsilon): # <<<
1569        """return all intersection points of the original path and the parallel path"""
1570
1571        # this code became necessary with introduction of mynormpathparam
1572        params = []
1573        oparams = []
1574        for nsp_i in range(len(orig_np)):
1575            for nsp_j in range(len(par_np)):
1576                for nspitem_i in range(len(orig_np[nsp_i])):
1577                    for nspitem_j in range(len(par_np[nsp_j])):
1578                        intsparams = orig_np[nsp_i][nspitem_i].intersect(par_np[nsp_j][nspitem_j], epsilon)
1579                        if intsparams:
1580                            for intsparam_i, intsparam_j in intsparams:
1581                                npp_i = mynormpathparam(orig_np, nsp_i, nspitem_i, intsparam_i)
1582                                npp_j = mynormpathparam(par_np, nsp_j, nspitem_j, intsparam_j)
1583
1584                                oparams.append(npp_i)
1585                                params.append(npp_j)
1586        return params, oparams
1587    # >>>
1588    def _can_continue(self, param1, param2, epsilon=None): # <<<
1589        """Test whether the parallel path can be continued at the param-pair (param1, param2)"""
1590        par_np = param1.normpath
1591        if epsilon is None:
1592            epsilon = par_np[0].epsilon
1593
1594        rot1, rot2 = par_np.rotation([param1, param2])
1595        orth1 = rot1.apply_pt(0, self.dist_pt) # directs away from original path (as seen from parallel path)
1596        tang2 = rot2.apply_pt(1, 0)
1597
1598        # the self-intersection is valid if the tangents
1599        # point into the correct direction or, for parallel tangents,
1600        # if the curvature is such that the on-going path does not
1601        # enter the region defined by dist
1602        proj = orth1[0]*tang2[0] + orth1[1]*tang2[1]
1603        if abs(proj) > epsilon: # the curves are not parallel
1604            # tang2 must go away from the original path
1605            return (proj > 0)
1606
1607        # tang1 and tang2 are parallel.
1608        curv1, curv2 = par_np.curvature_pt([param1, param2])
1609
1610        # We need to treat also cases where the params are nspitem-endpoints.
1611        # There, we know that the tangents are continuous, but the curvature is
1612        # not necessarily continuous. We have to test whether the curve *after*
1613        # param2 has curvature such that it enters the forbidden side of the
1614        # curve after param1
1615        if param1.is_end_of_nspitem(epsilon):
1616            curv1 = par_np.curvature_pt([param1.larger_equiv(epsilon)])[0]
1617        if param2.is_end_of_nspitem(epsilon):
1618            curv2 = par_np.curvature_pt([param2.larger_equiv(epsilon)])[0]
1619
1620        tang1 = rot1.apply_pt(1, 0)
1621        running_back = (tang1[0]*tang2[0] + tang1[1]*tang2[1] < 0)
1622        if running_back:
1623            # the second curve is running "back" -- the curvature sign appears to be switched
1624            curv2 = -curv2
1625            # endpoints of normsubpaths must be treated differently:
1626
1627        if (not running_back) and param1.is_end_of_nsp(epsilon):
1628            return True
1629
1630        if curv1 == curv2:
1631            raise IntersectionError("Cannot determine whether curves intersect (parallel and equally curved)")
1632
1633        if self.dist_pt > 0:
1634            return (curv2 > curv1)
1635        else:
1636            return (curv2 < curv1)
1637    # >>>
1638    def _skip_intersection_doublet(self, npp_i, npp_j, parampairs, epsilon): # <<<
1639        # An intersection point that lies exactly between two nsp-items can occur twice or more
1640        # times if we calculate all mutual intersections. We should take only
1641        # one such parameter pair, namely the one with smallest first and
1642        # largest last param.
1643        result = False
1644        delete_keys = []
1645        delete_values = []
1646        # TODO: improve complexity?
1647        for pi, pj in parampairs.items():
1648            if npp_i.is_equiv(pi, epsilon) and npp_j.is_equiv(pj, epsilon):
1649                #print("double pair: ", npp_i, npp_j, pi, pj)
1650                #print("... replacing ", pi, parampairs[pi], "by", min(npp_i, pi), max(npp_j, pj))
1651                delete_keys.append(pi)
1652                delete_values.append(pj)
1653                result = True # we have already added this one
1654        newkey = min([npp_i] + delete_keys)
1655        newval = max([npp_j] + delete_values)
1656        for pi in delete_keys:
1657            del parampairs[pi]
1658        parampairs[newkey] = newval
1659        return result
1660    # >>>
1661    def _elim_intersection_doublets(self, parampairs, epsilon): # <<<
1662        # It may always happen that three intersections coincide. (It will
1663        # occur often with degenerate distances for technical designs such as
1664        # those used in microfluidics). We then have two equivalent pairs in our
1665        # forward list, and we must throw away one of them.
1666        # One of them is indeed forbidden by the _can_continue of the other.
1667
1668        # TODO implement this
1669
1670        keys = list(parampairs.keys())
1671        n = len(keys)
1672        for i in range(n):
1673            start = "equivalent pairs\n"
1674            for j in range(i+1, n):
1675                key1, key2 = keys[i], keys[j]
1676                npp1 = parampairs[key1]
1677                npp2 = parampairs[key2]
1678                #assert key1.is_equiv(npp1, epsilon)
1679                #if not key2.is_equiv(npp2, epsilon):
1680                #    np = key2.normpath
1681                #    print(np.at_pt(key2), np.at_pt(npp2), _length_pt(np, key2, npp2)/epsilon)
1682                #assert key2.is_equiv(npp2, epsilon)
1683                if ((key1.is_equiv(key2, epsilon) and npp1.is_equiv(npp2, epsilon)) or
1684                    (key1.is_equiv(npp2, epsilon) and npp1.is_equiv(key2, epsilon))):
1685                    print(start,"pair: ", key1, npp1, " and ", key2, npp2)
1686                    start = ""
1687            if not start:
1688                print()
1689        return parampairs
1690    # >>>
1691    def _between_paths(self, pos, par2orig, epsilon): # <<<
1692        """Tests whether the given point (pos) is found in the forbidden zone between an original and a parallel nsp-item (these are in par2orig)
1693
1694        The test uses epsilon close to the original/parallel path, and sharp comparison at their ends."""
1695        dist = self.dist_pt
1696        for par_nspitem in par2orig:
1697            origobj = par2orig[par_nspitem]
1698            if isinstance(origobj, normpath.normline_pt):
1699                rot = origobj.rotation([0])[0]
1700                t, s = intersection(pos, origobj.atbegin_pt(), rot.apply_pt(0, mathutils.sign(dist)), rot.apply_pt(origobj.arclen_pt(epsilon), 0))
1701                if 0 <= s <= 1 and -abs(dist)+epsilon < t < -epsilon:
1702                    return True
1703            elif isinstance(origobj, normpath.normcurve_pt):
1704                # TODO: implement this
1705                # TODO: pre-sort par2orig as a list to fasten up this code
1706                pass
1707            else:
1708                cx, cy = origobj
1709                if math.hypot(pos[0]-cx, pos[1]-cy) < abs(dist) - epsilon:
1710                    if self.dist_pt > 0: # running around (cx,cy) in the negative sense (see _path_around_corner)
1711                        x0, y0 = par_nspitem.atend_pt()
1712                        x1, y1 = par_nspitem.atbegin_pt()
1713                    else: # running around (cx,cy) in the positive sense
1714                        x0, y0 = par_nspitem.atbegin_pt()
1715                        x1, y1 = par_nspitem.atend_pt()
1716                    t0, s0 = intersection(pos, (cx, cy), (-y0+cy, x0-cx), (x0-cx, y0-cy))
1717                    t1, s1 = intersection(pos, (cx, cy), ( y1-cy,-x1+cx), (x1-cx, y1-cy))
1718                    if t0 <= 0 and s0 >= 0 and t1 <= 0 and s1 >= 0:
1719                        return True
1720        return False
1721    # >>>
1722
1723# >>>
1724
1725parallel.clear = attr.clearclass(parallel)
1726
1727class linesmoothed(baseclasses.deformer): # <<<
1728
1729    def __init__(self, tension=1, atleast=False, lcurl=1, rcurl=1):
1730        """Tension and atleast control the tension of the replacement curves.
1731        l/rcurl control the curlynesses at (possible) endpoints. If a curl is
1732        set to None, the angle is taken from the original path."""
1733        if atleast:
1734            self.tension = -abs(tension)
1735        else:
1736            self.tension = abs(tension)
1737        self.lcurl = lcurl
1738        self.rcurl = rcurl
1739
1740    def __call__(self, tension=_marker, atleast=_marker, lcurl=_marker, rcurl=_marker):
1741        if tension is _marker:
1742            tension = self.tension
1743        if atleast is _marker:
1744            atleast = (self.tension < 0)
1745        if lcurl is _marker:
1746            lcurl = self.lcurl
1747        if rcurl is _marker:
1748            rcurl = self.rcurl
1749        return linesmoothed(tension, atleast, lcurl, rcurl)
1750
1751    def deform(self, basepath):
1752        newnp = normpath.normpath()
1753        for nsp in basepath.normpath().normsubpaths:
1754            newnp += self.deformsubpath(nsp)
1755        return newnp
1756
1757    def deformsubpath(self, nsp):
1758        from .metapost import path as mppath
1759        """Returns a path/normpath from the points in the given normsubpath"""
1760        # TODO: epsilon ?
1761        knots = []
1762
1763        # first point
1764        x_pt, y_pt = nsp.atbegin_pt()
1765        if nsp.closed:
1766            knots.append(mppath.smoothknot_pt(x_pt, y_pt))
1767        elif self.lcurl is None:
1768            rot = nsp.rotation([0])[0]
1769            dx, dy = rot.apply_pt(1, 0)
1770            angle = math.atan2(dy, dx)
1771            knots.append(mppath.beginknot_pt(x_pt, y_pt, angle=angle))
1772        else:
1773            knots.append(mppath.beginknot_pt(x_pt, y_pt, curl=self.lcurl))
1774
1775        # intermediate points:
1776        for npelem in nsp[:-1]:
1777            knots.append(mppath.tensioncurve(self.tension))
1778            knots.append(mppath.smoothknot_pt(*npelem.atend_pt()))
1779
1780        # last point
1781        knots.append(mppath.tensioncurve(self.tension))
1782        x_pt, y_pt = nsp.atend_pt()
1783        if nsp.closed:
1784            pass
1785        elif self.rcurl is None:
1786            rot = nsp.rotation([len(nsp)])[0]
1787            dx, dy = rot.apply_pt(1, 0)
1788            angle = math.atan2(dy, dx)
1789            knots.append(mppath.endknot_pt(x_pt, y_pt, angle=angle))
1790        else:
1791            knots.append(mppath.endknot_pt(x_pt, y_pt, curl=self.rcurl))
1792
1793        return mppath.path(knots)
1794# >>>
1795
1796linesmoothed.clear = attr.clearclass(linesmoothed)
1797
1798
1799# vim:foldmethod=marker:foldmarker=<<<,>>>
1800