set the raggedskip to 1fil for the text.flushalignments
[PyX/mjg.git] / pyx / deformer.py
blob9a6236d296a3ff667a4d1bb3d94160ebce2879f6
1 #!/usr/bin/env python
2 # -*- coding: ISO-8859-1 -*-
5 # Copyright (C) 2003-2005 Michael Schindler <m-schindler@users.sourceforge.net>
6 # Copyright (C) 2003-2005 André Wobst <wobsta@users.sourceforge.net>
8 # This file is part of PyX (http://pyx.sourceforge.net/).
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.
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.
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
25 import math, warnings
26 import attr, mathutils, path, normpath, unit, color
27 from path import degrees
29 # specific exception for an invalid parameterization point
30 # used in parallel
31 class InvalidParamException(Exception):
33 def __init__(self, param):
34 self.normsubpathitemparam = param
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
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
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)
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
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])
71 return (d1, g1, f1, e, f2, g2, d2)
72 # >>>
74 def controldists_from_endgeometry_pt(A, B, tangA, tangB, curvA, curvB, epsilon): # <<<
76 """For a curve with given tangents and curvatures at the endpoints this gives the distances between the controlpoints
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.
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 """
89 # these two thresholds are dimensionless, not lengths:
90 fallback_threshold = 1.0e-3
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])
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.
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)]
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)]
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)]
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)]
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)]
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)]
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]
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))
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)
202 solutions.sort(mycmp)
204 # XXX should the solutions's list also be unique'fied?
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."
211 return solutions
212 # >>>
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 # >>>
221 def intersection(A, D, tangA, tangD): # <<<
223 """returns the intersection parameters of two evens
225 they are defined by:
226 x(t) = A + t * tangA
227 x(s) = D + s * tangD
229 det = -tangA[0] * tangD[1] + tangA[1] * tangD[0]
230 try:
231 1.0 / det
232 except ArithmeticError:
233 return None, None
235 DA = D[0] - A[0], D[1] - A[1]
237 t = (-tangD[1]*DA[0] + tangD[0]*DA[1]) / det
238 s = (-tangA[1]*DA[0] + tangA[0]*DA[1]) / det
240 return t, s
241 # >>>
243 class deformer(attr.attr):
245 def deform (self, basepath):
246 return basepath
248 class cycloid(deformer): # <<<
249 """Wraps a cycloid around a path.
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
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
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
290 return cycloid(radius=radius, halfloops=halfloops, skipfirst=skipfirst, skiplast=skiplast,
291 curvesperhloop=curvesperhloop, sign=sign, turnangle=turnangle)
293 def deform(self, basepath):
294 resultnormsubpaths = [self.deformsubpath(nsp) for nsp in basepath.normpath().normsubpaths]
295 return normpath.normpath(resultnormsubpaths)
297 def deformsubpath(self, normsubpath):
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)
305 cosTurn = math.cos(turnangle)
306 sinTurn = math.sin(turnangle)
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
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))
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]
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]
340 # The point on the cycloid, in the basepath's local coordinate system
341 baseZ, baseY = 0, radius*math.sin(phi)
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
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
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
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))
370 if len(points) <= 1:
371 warnings.warn("normsubpath is too short for deformation with cycloid -- skipping...")
372 return normsubpath
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)
387 # That's it
388 return normpath.normsubpath(normsubpathitems, epsilon=normsubpath.epsilon)
389 # >>>
391 cycloid.clear = attr.clearclass(cycloid)
393 class smoothed(deformer): # <<<
395 """Bends corners in a normpath.
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
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
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)
429 def deform(self, basepath):
430 return normpath.normpath([self.deformsubpath(normsubpath)
431 for normsubpath in basepath.normpath().normsubpaths])
433 def deformsubpath(self, normsubpath):
434 radius_pt = unit.topt(self.radius)
435 epsilon = normsubpath.epsilon
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
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))
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)
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]
474 # insert the middle segment
475 newnormsubpath.append(thisnormsubpathitem.segments(thisparams)[0])
477 # insert replacement curves for the corners
478 if next or normsubpath.closed:
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
484 if (isinstance(thisnormsubpathitem, normpath.normline_pt) and
485 isinstance(nextnormsubpathitem, normpath.normline_pt)):
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)
492 p1 = thisnormsubpathitem.at_pt([thisparams[1]])[0]
493 p2 = nextnormsubpathitem.at_pt([nextparams[0]])[0]
495 newnormsubpath.append(normpath.normcurve_pt(*(d1 + g1 + f1 + e)))
496 newnormsubpath.append(normpath.normcurve_pt(*(e + f2 + g2 + d2)))
498 else:
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
507 # TODO: more intelligent fallbacks:
508 # circle -> line
509 # circle -> circle
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)
521 # get the length of the control "arms"
522 a, d = controldists_from_endgeometry_pt(p1, p2, t1, t2, c1, c2, epsilon=epsilon)[0]
524 if a >= 0 and d >= 0:
525 # avoid curves with invalid parameterization
526 a = max(a, epsilon)
527 d = max(d, epsilon)
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))
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))
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]
556 newnormsubpath.append(normpath.normcurve_pt(*(p1 + q1 + q2 + p2)))
558 if normsubpath.closed:
559 newnormsubpath.close()
560 return newnormsubpath
562 # >>>
564 smoothed.clear = attr.clearclass(smoothed)
566 class parallel(deformer): # <<<
568 """creates a parallel normpath with constant distance to the original normpath
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.
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
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?
595 # TODO for the test cases:
596 # * improve the intersection of nearly degenerate paths
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
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
627 return parallel(distance=distance, relerr=relerr,
628 sharpoutercorners=sharpoutercorners,
629 dointersection=dointersection,
630 checkdistanceparams=checkdistanceparams,
631 lookforcurvatures=lookforcurvatures,
632 debug=debug)
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
643 def deformsubpath(self, orig_nsp): # <<<
645 """returns a list of normsubpaths building the parallel curve"""
647 dist = self.dist_pt
648 epsilon = orig_nsp.epsilon
650 # avoid too small dists: we would run into instabilities
651 if abs(dist) < abs(epsilon):
652 return orig_nsp
654 result = normpath.normpath()
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]
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
672 prev_tangent = prev_rotation.apply_pt(1, 0)
673 next_tangent = next_rotation.apply_pt(1, 0)
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])
696 next_parallel_normpath = self.deformsubpath(normpath.normsubpath(segments, epsilon=epsilon))
698 if not (next_parallel_normpath.normsubpaths and next_parallel_normpath[0].normsubpathitems):
699 continue
701 # this starts the whole normpath
702 if not result.normsubpaths:
703 result = next_parallel_normpath
704 continue
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])
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)
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)
747 # end here if nothing has been found so far
748 if not (result.normsubpaths and result[-1].normsubpathitems):
749 return result
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
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]
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])
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)
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)
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:]
806 if self.dointersection:
807 result = self.rebuild_intersected_normpath(result, normpath.normpath([orig_nsp]), epsilon)
809 return result
810 # >>>
811 def deformsubpathitem(self, nspitem, epsilon): # <<<
813 """Returns a parallel normpath for a single normsubpathitem
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"""
820 dist = self.dist_pt
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)
831 # for a curve we have to check if the curvatures
832 # cross the singular value 1/dist
833 crossings = self.distcrossingparameters(nspitem, epsilon)
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()
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)
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)
862 return result
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()
881 # >>>
882 def deformnicecurve(self, normcurve, epsilon, startparam=0.0, endparam=1.0): # <<<
884 """Returns a parallel normsubpath for the normcurve.
886 This routine assumes that the normcurve is everywhere
887 'less' curved than 1/dist and contains no point with an
888 invalid parameterization
890 dist = self.dist_pt
891 T_threshold = 1.0e-5
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)
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]
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)
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)
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)
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)
968 return result
969 # >>>
971 def distchecked(self, orig_normcurve, parallel_normsubpath, epsilon, tstart, tend): # <<<
973 """Checks the distances between orig_normcurve and parallel_normsubpath
975 The checking is done at parameters self.checkdistanceparams of orig_normcurve."""
977 dist = self.dist_pt
978 # do not look closer than epsilon:
979 dist_relerr = mathutils.sign(dist) * max(abs(self.relerr*dist), epsilon)
981 checkdistanceparams = [tstart + (tend-tstart)*t for t in self.checkdistanceparams]
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)
991 normal = rotation.apply_pt(0, 1)
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)
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])]
1004 if (not distances) or (abs(min(distances) - abs(dist)) > abs(dist_relerr)):
1005 return 0
1007 return 1
1008 # >>>
1009 def distcrossingparameters(self, normcurve, epsilon, tstart=0, tend=1): # <<<
1011 """Returns a list of parameters where the curvature is 1/distance"""
1013 dist = self.dist_pt
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
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)
1027 # break everything at invalid curvatures
1028 for param, curv in zip(params, curvs):
1029 if curv is normpath.invalid:
1030 raise InvalidParamException(param)
1032 parampairs = zip(params[:-1], params[1:])
1033 curvpairs = zip(curvs[:-1], curvs[1:])
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]
1047 if middleradius is normpath.invalid:
1048 raise InvalidParamException(middleparam)
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
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 # >>>
1087 def normpath_selfintersections(self, np, epsilon): # <<<
1089 """return all self-intersection points of normpath np.
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"""
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
1127 # >>>
1128 def can_continue(self, par_np, param1, param2): # <<<
1129 dist = self.dist_pt
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])
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): # <<<
1158 dist = self.dist_pt
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]
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 # >>>
1176 result = normpath.normpath()
1177 if not selfintparams:
1178 if origintparams:
1179 return result
1180 else:
1181 return par_np
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])))
1189 allparams = selfintparams + origintparams + beginparams + endparams
1190 allparams.sort()
1191 allparamindices = {}
1192 for i, param in enumerate(allparams):
1193 allparamindices[param] = i
1195 done = {}
1196 for param in allparams:
1197 done[param] = 0
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]
1211 lastp = startp
1212 currentp = allparams[allparamindices[startp] + 1]
1213 result = []
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 # >>>
1247 # first the paths that start at the beginning of a subnormpath:
1248 for startp in beginparams + selfintparams:
1249 if done[startp]:
1250 continue
1252 parampairs = trial_parampairs(startp)
1253 if not parampairs:
1254 continue
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)
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
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()
1284 result.extend([add_nsp])
1286 return result
1288 # >>>
1290 # >>>
1292 parallel.clear = attr.clearclass(parallel)
1294 # vim:foldmethod=marker:foldmarker=<<<,>>>