Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/fontTools/misc/bezierTools.py: 19%

1# -*- coding: utf-8 -*-

2"""fontTools.misc.bezierTools.py -- tools for working with Bezier path segments.

3"""

5from fontTools.misc.arrayTools import calcBounds, sectRect, rectArea

6from fontTools.misc.transform import Identity

7import math

8from collections import namedtuple

10try:

11 import cython

13 COMPILED = cython.compiled

14except (AttributeError, ImportError):

15 # if cython not installed, use mock module with no-op decorators and types

16 from fontTools.misc import cython

18 COMPILED = False

21Intersection = namedtuple("Intersection", ["pt", "t1", "t2"])

24__all__ = [

25 "approximateCubicArcLength",

26 "approximateCubicArcLengthC",

27 "approximateQuadraticArcLength",

28 "approximateQuadraticArcLengthC",

29 "calcCubicArcLength",

30 "calcCubicArcLengthC",

31 "calcQuadraticArcLength",

32 "calcQuadraticArcLengthC",

33 "calcCubicBounds",

34 "calcQuadraticBounds",

35 "splitLine",

36 "splitQuadratic",

37 "splitCubic",

38 "splitQuadraticAtT",

39 "splitCubicAtT",

40 "splitCubicAtTC",

41 "splitCubicIntoTwoAtTC",

42 "solveQuadratic",

43 "solveCubic",

44 "quadraticPointAtT",

45 "cubicPointAtT",

46 "cubicPointAtTC",

47 "linePointAtT",

48 "segmentPointAtT",

49 "lineLineIntersections",

50 "curveLineIntersections",

51 "curveCurveIntersections",

52 "segmentSegmentIntersections",

53]

56def calcCubicArcLength(pt1, pt2, pt3, pt4, tolerance=0.005):

57 """Calculates the arc length for a cubic Bezier segment.

59 Whereas :func:`approximateCubicArcLength` approximates the length, this

60 function calculates it by "measuring", recursively dividing the curve

61 until the divided segments are shorter than ``tolerance``.

63 Args:

64 pt1,pt2,pt3,pt4: Control points of the Bezier as 2D tuples.

65 tolerance: Controls the precision of the calcuation.

67 Returns:

68 Arc length value.

69 """

70 return calcCubicArcLengthC(

71 complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4), tolerance

72 )

75def _split_cubic_into_two(p0, p1, p2, p3):

76 mid = (p0 + 3 * (p1 + p2) + p3) * 0.125

77 deriv3 = (p3 + p2 - p1 - p0) * 0.125

78 return (

79 (p0, (p0 + p1) * 0.5, mid - deriv3, mid),

80 (mid, mid + deriv3, (p2 + p3) * 0.5, p3),

81 )

84@cython.returns(cython.double)

85@cython.locals(

86 p0=cython.complex,

87 p1=cython.complex,

88 p2=cython.complex,

89 p3=cython.complex,

90)

91@cython.locals(mult=cython.double, arch=cython.double, box=cython.double)

92def _calcCubicArcLengthCRecurse(mult, p0, p1, p2, p3):

93 arch = abs(p0 - p3)

94 box = abs(p0 - p1) + abs(p1 - p2) + abs(p2 - p3)

95 if arch * mult >= box:

96 return (arch + box) * 0.5

97 else:

98 one, two = _split_cubic_into_two(p0, p1, p2, p3)

99 return _calcCubicArcLengthCRecurse(mult, *one) + _calcCubicArcLengthCRecurse(

100 mult, *two

101 )

102

103

104@cython.returns(cython.double)

105@cython.locals(

106 pt1=cython.complex,

107 pt2=cython.complex,

108 pt3=cython.complex,

109 pt4=cython.complex,

110)

111@cython.locals(

112 tolerance=cython.double,

113 mult=cython.double,

114)

115def calcCubicArcLengthC(pt1, pt2, pt3, pt4, tolerance=0.005):

116 """Calculates the arc length for a cubic Bezier segment.

117

118 Args:

119 pt1,pt2,pt3,pt4: Control points of the Bezier as complex numbers.

120 tolerance: Controls the precision of the calcuation.

121

122 Returns:

123 Arc length value.

124 """

125 mult = 1.0 + 1.5 * tolerance # The 1.5 is a empirical hack; no math

126 return _calcCubicArcLengthCRecurse(mult, pt1, pt2, pt3, pt4)

127

128

129epsilonDigits = 6

130epsilon = 1e-10

131

132

133@cython.cfunc

134@cython.inline

135@cython.returns(cython.double)

136@cython.locals(v1=cython.complex, v2=cython.complex)

137def _dot(v1, v2):

138 return (v1 * v2.conjugate()).real

139

140

141@cython.cfunc

142@cython.inline

143@cython.returns(cython.double)

144@cython.locals(x=cython.complex)

145def _intSecAtan(x):

146 # In : sympy.integrate(sp.sec(sp.atan(x)))

147 # Out: x*sqrt(x**2 + 1)/2 + asinh(x)/2

148 return x * math.sqrt(x**2 + 1) / 2 + math.asinh(x) / 2

149

150

151def calcQuadraticArcLength(pt1, pt2, pt3):

152 """Calculates the arc length for a quadratic Bezier segment.

153

154 Args:

155 pt1: Start point of the Bezier as 2D tuple.

156 pt2: Handle point of the Bezier as 2D tuple.

157 pt3: End point of the Bezier as 2D tuple.

158

159 Returns:

160 Arc length value.

161

162 Example::

163

164 >>> calcQuadraticArcLength((0, 0), (0, 0), (0, 0)) # empty segment

165 0.0

166 >>> calcQuadraticArcLength((0, 0), (50, 0), (80, 0)) # collinear points

167 80.0

168 >>> calcQuadraticArcLength((0, 0), (0, 50), (0, 80)) # collinear points vertical

169 80.0

170 >>> calcQuadraticArcLength((0, 0), (50, 20), (100, 40)) # collinear points

171 107.70329614269008

172 >>> calcQuadraticArcLength((0, 0), (0, 100), (100, 0))

173 154.02976155645263

174 >>> calcQuadraticArcLength((0, 0), (0, 50), (100, 0))

175 120.21581243984076

176 >>> calcQuadraticArcLength((0, 0), (50, -10), (80, 50))

177 102.53273816445825

178 >>> calcQuadraticArcLength((0, 0), (40, 0), (-40, 0)) # collinear points, control point outside

179 66.66666666666667

180 >>> calcQuadraticArcLength((0, 0), (40, 0), (0, 0)) # collinear points, looping back

181 40.0

182 """

183 return calcQuadraticArcLengthC(complex(*pt1), complex(*pt2), complex(*pt3))

184

185

186@cython.returns(cython.double)

187@cython.locals(

188 pt1=cython.complex,

189 pt2=cython.complex,

190 pt3=cython.complex,

191 d0=cython.complex,

192 d1=cython.complex,

193 d=cython.complex,

194 n=cython.complex,

195)

196@cython.locals(

197 scale=cython.double,

198 origDist=cython.double,

199 a=cython.double,

200 b=cython.double,

201 x0=cython.double,

202 x1=cython.double,

203 Len=cython.double,

204)

205def calcQuadraticArcLengthC(pt1, pt2, pt3):

206 """Calculates the arc length for a quadratic Bezier segment.

207

208 Args:

209 pt1: Start point of the Bezier as a complex number.

210 pt2: Handle point of the Bezier as a complex number.

211 pt3: End point of the Bezier as a complex number.

212

213 Returns:

214 Arc length value.

215 """

216 # Analytical solution to the length of a quadratic bezier.

217 # Documentation: https://github.com/fonttools/fonttools/issues/3055

218 d0 = pt2 - pt1

219 d1 = pt3 - pt2

220 d = d1 - d0

221 n = d * 1j

222 scale = abs(n)

223 if scale == 0.0:

224 return abs(pt3 - pt1)

225 origDist = _dot(n, d0)

226 if abs(origDist) < epsilon:

227 if _dot(d0, d1) >= 0:

228 return abs(pt3 - pt1)

229 a, b = abs(d0), abs(d1)

230 return (a * a + b * b) / (a + b)

231 x0 = _dot(d, d0) / origDist

232 x1 = _dot(d, d1) / origDist

233 Len = abs(2 * (_intSecAtan(x1) - _intSecAtan(x0)) * origDist / (scale * (x1 - x0)))

234 return Len

235

236

237def approximateQuadraticArcLength(pt1, pt2, pt3):

238 """Calculates the arc length for a quadratic Bezier segment.

239

240 Uses Gauss-Legendre quadrature for a branch-free approximation.

241 See :func:`calcQuadraticArcLength` for a slower but more accurate result.

242

243 Args:

244 pt1: Start point of the Bezier as 2D tuple.

245 pt2: Handle point of the Bezier as 2D tuple.

246 pt3: End point of the Bezier as 2D tuple.

247

248 Returns:

249 Approximate arc length value.

250 """

251 return approximateQuadraticArcLengthC(complex(*pt1), complex(*pt2), complex(*pt3))

252

253

254@cython.returns(cython.double)

255@cython.locals(

256 pt1=cython.complex,

257 pt2=cython.complex,

258 pt3=cython.complex,

259)

260@cython.locals(

261 v0=cython.double,

262 v1=cython.double,

263 v2=cython.double,

264)

265def approximateQuadraticArcLengthC(pt1, pt2, pt3):

266 """Calculates the arc length for a quadratic Bezier segment.

267

268 Uses Gauss-Legendre quadrature for a branch-free approximation.

269 See :func:`calcQuadraticArcLength` for a slower but more accurate result.

270

271 Args:

272 pt1: Start point of the Bezier as a complex number.

273 pt2: Handle point of the Bezier as a complex number.

274 pt3: End point of the Bezier as a complex number.

275

276 Returns:

277 Approximate arc length value.

278 """

279 # This, essentially, approximates the length-of-derivative function

280 # to be integrated with the best-matching fifth-degree polynomial

281 # approximation of it.

282 #

283 # https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Legendre_quadrature

284

285 # abs(BezierCurveC[2].diff(t).subs({t:T})) for T in sorted(.5, .5±sqrt(3/5)/2),

286 # weighted 5/18, 8/18, 5/18 respectively.

287 v0 = abs(

288 -0.492943519233745 * pt1 + 0.430331482911935 * pt2 + 0.0626120363218102 * pt3

289 )

290 v1 = abs(pt3 - pt1) * 0.4444444444444444

291 v2 = abs(

292 -0.0626120363218102 * pt1 - 0.430331482911935 * pt2 + 0.492943519233745 * pt3

293 )

294

295 return v0 + v1 + v2

296

297

298def calcQuadraticBounds(pt1, pt2, pt3):

299 """Calculates the bounding rectangle for a quadratic Bezier segment.

300

301 Args:

302 pt1: Start point of the Bezier as a 2D tuple.

303 pt2: Handle point of the Bezier as a 2D tuple.

304 pt3: End point of the Bezier as a 2D tuple.

305

306 Returns:

307 A four-item tuple representing the bounding rectangle ``(xMin, yMin, xMax, yMax)``.

308

309 Example::

310

311 >>> calcQuadraticBounds((0, 0), (50, 100), (100, 0))

312 (0, 0, 100, 50.0)

313 >>> calcQuadraticBounds((0, 0), (100, 0), (100, 100))

314 (0.0, 0.0, 100, 100)

315 """

316 (ax, ay), (bx, by), (cx, cy) = calcQuadraticParameters(pt1, pt2, pt3)

317 ax2 = ax * 2.0

318 ay2 = ay * 2.0

319 roots = []

320 if ax2 != 0:

321 roots.append(-bx / ax2)

322 if ay2 != 0:

323 roots.append(-by / ay2)

324 points = [

325 (ax * t * t + bx * t + cx, ay * t * t + by * t + cy)

326 for t in roots

327 if 0 <= t < 1

328 ] + [pt1, pt3]

329 return calcBounds(points)

330

331

332def approximateCubicArcLength(pt1, pt2, pt3, pt4):

333 """Approximates the arc length for a cubic Bezier segment.

334

335 Uses Gauss-Lobatto quadrature with n=5 points to approximate arc length.

336 See :func:`calcCubicArcLength` for a slower but more accurate result.

337

338 Args:

339 pt1,pt2,pt3,pt4: Control points of the Bezier as 2D tuples.

340

341 Returns:

342 Arc length value.

343

344 Example::

345

346 >>> approximateCubicArcLength((0, 0), (25, 100), (75, 100), (100, 0))

347 190.04332968932817

348 >>> approximateCubicArcLength((0, 0), (50, 0), (100, 50), (100, 100))

349 154.8852074945903

350 >>> approximateCubicArcLength((0, 0), (50, 0), (100, 0), (150, 0)) # line; exact result should be 150.

351 149.99999999999991

352 >>> approximateCubicArcLength((0, 0), (50, 0), (100, 0), (-50, 0)) # cusp; exact result should be 150.

353 136.9267662156362

354 >>> approximateCubicArcLength((0, 0), (50, 0), (100, -50), (-50, 0)) # cusp

355 154.80848416537057

356 """

357 return approximateCubicArcLengthC(

358 complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4)

359 )

360

361

362@cython.returns(cython.double)

363@cython.locals(

364 pt1=cython.complex,

365 pt2=cython.complex,

366 pt3=cython.complex,

367 pt4=cython.complex,

368)

369@cython.locals(

370 v0=cython.double,

371 v1=cython.double,

372 v2=cython.double,

373 v3=cython.double,

374 v4=cython.double,

375)

376def approximateCubicArcLengthC(pt1, pt2, pt3, pt4):

377 """Approximates the arc length for a cubic Bezier segment.

378

379 Args:

380 pt1,pt2,pt3,pt4: Control points of the Bezier as complex numbers.

381

382 Returns:

383 Arc length value.

384 """

385 # This, essentially, approximates the length-of-derivative function

386 # to be integrated with the best-matching seventh-degree polynomial

387 # approximation of it.

388 #

389 # https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules

390

391 # abs(BezierCurveC[3].diff(t).subs({t:T})) for T in sorted(0, .5±(3/7)**.5/2, .5, 1),

392 # weighted 1/20, 49/180, 32/90, 49/180, 1/20 respectively.

393 v0 = abs(pt2 - pt1) * 0.15

394 v1 = abs(

395 -0.558983582205757 * pt1

396 + 0.325650248872424 * pt2

397 + 0.208983582205757 * pt3

398 + 0.024349751127576 * pt4

399 )

400 v2 = abs(pt4 - pt1 + pt3 - pt2) * 0.26666666666666666

401 v3 = abs(

402 -0.024349751127576 * pt1

403 - 0.208983582205757 * pt2

404 - 0.325650248872424 * pt3

405 + 0.558983582205757 * pt4

406 )

407 v4 = abs(pt4 - pt3) * 0.15

408

409 return v0 + v1 + v2 + v3 + v4

410

411

412def calcCubicBounds(pt1, pt2, pt3, pt4):

413 """Calculates the bounding rectangle for a quadratic Bezier segment.

414

415 Args:

416 pt1,pt2,pt3,pt4: Control points of the Bezier as 2D tuples.

417

418 Returns:

419 A four-item tuple representing the bounding rectangle ``(xMin, yMin, xMax, yMax)``.

420

421 Example::

422

423 >>> calcCubicBounds((0, 0), (25, 100), (75, 100), (100, 0))

424 (0, 0, 100, 75.0)

425 >>> calcCubicBounds((0, 0), (50, 0), (100, 50), (100, 100))

426 (0.0, 0.0, 100, 100)

427 >>> print("%f %f %f %f" % calcCubicBounds((50, 0), (0, 100), (100, 100), (50, 0)))

428 35.566243 0.000000 64.433757 75.000000

429 """

430 (ax, ay), (bx, by), (cx, cy), (dx, dy) = calcCubicParameters(pt1, pt2, pt3, pt4)

431 # calc first derivative

432 ax3 = ax * 3.0

433 ay3 = ay * 3.0

434 bx2 = bx * 2.0

435 by2 = by * 2.0

436 xRoots = [t for t in solveQuadratic(ax3, bx2, cx) if 0 <= t < 1]

437 yRoots = [t for t in solveQuadratic(ay3, by2, cy) if 0 <= t < 1]

438 roots = xRoots + yRoots

439

440 points = [

441 (

442 ax * t * t * t + bx * t * t + cx * t + dx,

443 ay * t * t * t + by * t * t + cy * t + dy,

444 )

445 for t in roots

446 ] + [pt1, pt4]

447 return calcBounds(points)

448

449

450def splitLine(pt1, pt2, where, isHorizontal):

451 """Split a line at a given coordinate.

452

453 Args:

454 pt1: Start point of line as 2D tuple.

455 pt2: End point of line as 2D tuple.

456 where: Position at which to split the line.

457 isHorizontal: Direction of the ray splitting the line. If true,

458 ``where`` is interpreted as a Y coordinate; if false, then

459 ``where`` is interpreted as an X coordinate.

460

461 Returns:

462 A list of two line segments (each line segment being two 2D tuples)

463 if the line was successfully split, or a list containing the original

464 line.

465

466 Example::

467

468 >>> printSegments(splitLine((0, 0), (100, 100), 50, True))

469 ((0, 0), (50, 50))

470 ((50, 50), (100, 100))

471 >>> printSegments(splitLine((0, 0), (100, 100), 100, True))

472 ((0, 0), (100, 100))

473 >>> printSegments(splitLine((0, 0), (100, 100), 0, True))

474 ((0, 0), (0, 0))

475 ((0, 0), (100, 100))

476 >>> printSegments(splitLine((0, 0), (100, 100), 0, False))

477 ((0, 0), (0, 0))

478 ((0, 0), (100, 100))

479 >>> printSegments(splitLine((100, 0), (0, 0), 50, False))

480 ((100, 0), (50, 0))

481 ((50, 0), (0, 0))

482 >>> printSegments(splitLine((0, 100), (0, 0), 50, True))

483 ((0, 100), (0, 50))

484 ((0, 50), (0, 0))

485 """

486 pt1x, pt1y = pt1

487 pt2x, pt2y = pt2

488

489 ax = pt2x - pt1x

490 ay = pt2y - pt1y

491

492 bx = pt1x

493 by = pt1y

494

495 a = (ax, ay)[isHorizontal]

496

497 if a == 0:

498 return [(pt1, pt2)]

499 t = (where - (bx, by)[isHorizontal]) / a

500 if 0 <= t < 1:

501 midPt = ax * t + bx, ay * t + by

502 return [(pt1, midPt), (midPt, pt2)]

503 else:

504 return [(pt1, pt2)]

505

506

507def splitQuadratic(pt1, pt2, pt3, where, isHorizontal):

508 """Split a quadratic Bezier curve at a given coordinate.

509

510 Args:

511 pt1,pt2,pt3: Control points of the Bezier as 2D tuples.

512 where: Position at which to split the curve.

513 isHorizontal: Direction of the ray splitting the curve. If true,

514 ``where`` is interpreted as a Y coordinate; if false, then

515 ``where`` is interpreted as an X coordinate.

516

517 Returns:

518 A list of two curve segments (each curve segment being three 2D tuples)

519 if the curve was successfully split, or a list containing the original

520 curve.

521

522 Example::

523

524 >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 150, False))

525 ((0, 0), (50, 100), (100, 0))

526 >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 50, False))

527 ((0, 0), (25, 50), (50, 50))

528 ((50, 50), (75, 50), (100, 0))

529 >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 25, False))

530 ((0, 0), (12.5, 25), (25, 37.5))

531 ((25, 37.5), (62.5, 75), (100, 0))

532 >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 25, True))

533 ((0, 0), (7.32233, 14.6447), (14.6447, 25))

534 ((14.6447, 25), (50, 75), (85.3553, 25))

535 ((85.3553, 25), (92.6777, 14.6447), (100, -7.10543e-15))

536 >>> # XXX I'm not at all sure if the following behavior is desirable:

537 >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 50, True))

538 ((0, 0), (25, 50), (50, 50))

539 ((50, 50), (50, 50), (50, 50))

540 ((50, 50), (75, 50), (100, 0))

541 """

542 a, b, c = calcQuadraticParameters(pt1, pt2, pt3)

543 solutions = solveQuadratic(

544 a[isHorizontal], b[isHorizontal], c[isHorizontal] - where

545 )

546 solutions = sorted(t for t in solutions if 0 <= t < 1)

547 if not solutions:

548 return [(pt1, pt2, pt3)]

549 return _splitQuadraticAtT(a, b, c, *solutions)

550

551

552def splitCubic(pt1, pt2, pt3, pt4, where, isHorizontal):

553 """Split a cubic Bezier curve at a given coordinate.

554

555 Args:

556 pt1,pt2,pt3,pt4: Control points of the Bezier as 2D tuples.

557 where: Position at which to split the curve.

558 isHorizontal: Direction of the ray splitting the curve. If true,

559 ``where`` is interpreted as a Y coordinate; if false, then

560 ``where`` is interpreted as an X coordinate.

561

562 Returns:

563 A list of two curve segments (each curve segment being four 2D tuples)

564 if the curve was successfully split, or a list containing the original

565 curve.

566

567 Example::

568

569 >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 150, False))

570 ((0, 0), (25, 100), (75, 100), (100, 0))

571 >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 50, False))

572 ((0, 0), (12.5, 50), (31.25, 75), (50, 75))

573 ((50, 75), (68.75, 75), (87.5, 50), (100, 0))

574 >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 25, True))

575 ((0, 0), (2.29379, 9.17517), (4.79804, 17.5085), (7.47414, 25))

576 ((7.47414, 25), (31.2886, 91.6667), (68.7114, 91.6667), (92.5259, 25))

577 ((92.5259, 25), (95.202, 17.5085), (97.7062, 9.17517), (100, 1.77636e-15))

578 """

579 a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)

580 solutions = solveCubic(

581 a[isHorizontal], b[isHorizontal], c[isHorizontal], d[isHorizontal] - where

582 )

583 solutions = sorted(t for t in solutions if 0 <= t < 1)

584 if not solutions:

585 return [(pt1, pt2, pt3, pt4)]

586 return _splitCubicAtT(a, b, c, d, *solutions)

587

588

589def splitQuadraticAtT(pt1, pt2, pt3, *ts):

590 """Split a quadratic Bezier curve at one or more values of t.

591

592 Args:

593 pt1,pt2,pt3: Control points of the Bezier as 2D tuples.

594 *ts: Positions at which to split the curve.

595

596 Returns:

597 A list of curve segments (each curve segment being three 2D tuples).

598

599 Examples::

600

601 >>> printSegments(splitQuadraticAtT((0, 0), (50, 100), (100, 0), 0.5))

602 ((0, 0), (25, 50), (50, 50))

603 ((50, 50), (75, 50), (100, 0))

604 >>> printSegments(splitQuadraticAtT((0, 0), (50, 100), (100, 0), 0.5, 0.75))

605 ((0, 0), (25, 50), (50, 50))

606 ((50, 50), (62.5, 50), (75, 37.5))

607 ((75, 37.5), (87.5, 25), (100, 0))

608 """

609 a, b, c = calcQuadraticParameters(pt1, pt2, pt3)

610 return _splitQuadraticAtT(a, b, c, *ts)

611

612

613def splitCubicAtT(pt1, pt2, pt3, pt4, *ts):

614 """Split a cubic Bezier curve at one or more values of t.

615

616 Args:

617 pt1,pt2,pt3,pt4: Control points of the Bezier as 2D tuples.

618 *ts: Positions at which to split the curve.

619

620 Returns:

621 A list of curve segments (each curve segment being four 2D tuples).

622

623 Examples::

624

625 >>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5))

626 ((0, 0), (12.5, 50), (31.25, 75), (50, 75))

627 ((50, 75), (68.75, 75), (87.5, 50), (100, 0))

628 >>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5, 0.75))

629 ((0, 0), (12.5, 50), (31.25, 75), (50, 75))

630 ((50, 75), (59.375, 75), (68.75, 68.75), (77.3438, 56.25))

631 ((77.3438, 56.25), (85.9375, 43.75), (93.75, 25), (100, 0))

632 """

633 a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)

634 return _splitCubicAtT(a, b, c, d, *ts)

635

636

637@cython.locals(

638 pt1=cython.complex,

639 pt2=cython.complex,

640 pt3=cython.complex,

641 pt4=cython.complex,

642 a=cython.complex,

643 b=cython.complex,

644 c=cython.complex,

645 d=cython.complex,

646)

647def splitCubicAtTC(pt1, pt2, pt3, pt4, *ts):

648 """Split a cubic Bezier curve at one or more values of t.

649

650 Args:

651 pt1,pt2,pt3,pt4: Control points of the Bezier as complex numbers..

652 *ts: Positions at which to split the curve.

653

654 Yields:

655 Curve segments (each curve segment being four complex numbers).

656 """

657 a, b, c, d = calcCubicParametersC(pt1, pt2, pt3, pt4)

658 yield from _splitCubicAtTC(a, b, c, d, *ts)

659

660

661@cython.returns(cython.complex)

662@cython.locals(

663 t=cython.double,

664 pt1=cython.complex,

665 pt2=cython.complex,

666 pt3=cython.complex,

667 pt4=cython.complex,

668 pointAtT=cython.complex,

669 off1=cython.complex,

670 off2=cython.complex,

671)

672@cython.locals(

673 t2=cython.double, _1_t=cython.double, _1_t_2=cython.double, _2_t_1_t=cython.double

674)

675def splitCubicIntoTwoAtTC(pt1, pt2, pt3, pt4, t):

676 """Split a cubic Bezier curve at t.

677

678 Args:

679 pt1,pt2,pt3,pt4: Control points of the Bezier as complex numbers.

680 t: Position at which to split the curve.

681

682 Returns:

683 A tuple of two curve segments (each curve segment being four complex numbers).

684 """

685 t2 = t * t

686 _1_t = 1 - t

687 _1_t_2 = _1_t * _1_t

688 _2_t_1_t = 2 * t * _1_t

689 pointAtT = (

690 _1_t_2 * _1_t * pt1 + 3 * (_1_t_2 * t * pt2 + _1_t * t2 * pt3) + t2 * t * pt4

691 )

692 off1 = _1_t_2 * pt1 + _2_t_1_t * pt2 + t2 * pt3

693 off2 = _1_t_2 * pt2 + _2_t_1_t * pt3 + t2 * pt4

694

695 pt2 = pt1 + (pt2 - pt1) * t

696 pt3 = pt4 + (pt3 - pt4) * _1_t

697

698 return ((pt1, pt2, off1, pointAtT), (pointAtT, off2, pt3, pt4))

699

700

701def _splitQuadraticAtT(a, b, c, *ts):

702 ts = list(ts)

703 segments = []

704 ts.insert(0, 0.0)

705 ts.append(1.0)

706 ax, ay = a

707 bx, by = b

708 cx, cy = c

709 for i in range(len(ts) - 1):

710 t1 = ts[i]

711 t2 = ts[i + 1]

712 delta = t2 - t1

713 # calc new a, b and c

714 delta_2 = delta * delta

715 a1x = ax * delta_2

716 a1y = ay * delta_2

717 b1x = (2 * ax * t1 + bx) * delta

718 b1y = (2 * ay * t1 + by) * delta

719 t1_2 = t1 * t1

720 c1x = ax * t1_2 + bx * t1 + cx

721 c1y = ay * t1_2 + by * t1 + cy

722

723 pt1, pt2, pt3 = calcQuadraticPoints((a1x, a1y), (b1x, b1y), (c1x, c1y))

724 segments.append((pt1, pt2, pt3))

725 return segments

726

727

728def _splitCubicAtT(a, b, c, d, *ts):

729 ts = list(ts)

730 ts.insert(0, 0.0)

731 ts.append(1.0)

732 segments = []

733 ax, ay = a

734 bx, by = b

735 cx, cy = c

736 dx, dy = d

737 for i in range(len(ts) - 1):

738 t1 = ts[i]

739 t2 = ts[i + 1]

740 delta = t2 - t1

741

742 delta_2 = delta * delta

743 delta_3 = delta * delta_2

744 t1_2 = t1 * t1

745 t1_3 = t1 * t1_2

746

747 # calc new a, b, c and d

748 a1x = ax * delta_3

749 a1y = ay * delta_3

750 b1x = (3 * ax * t1 + bx) * delta_2

751 b1y = (3 * ay * t1 + by) * delta_2

752 c1x = (2 * bx * t1 + cx + 3 * ax * t1_2) * delta

753 c1y = (2 * by * t1 + cy + 3 * ay * t1_2) * delta

754 d1x = ax * t1_3 + bx * t1_2 + cx * t1 + dx

755 d1y = ay * t1_3 + by * t1_2 + cy * t1 + dy

756 pt1, pt2, pt3, pt4 = calcCubicPoints(

757 (a1x, a1y), (b1x, b1y), (c1x, c1y), (d1x, d1y)

758 )

759 segments.append((pt1, pt2, pt3, pt4))

760 return segments

761

762

763@cython.locals(

764 a=cython.complex,

765 b=cython.complex,

766 c=cython.complex,

767 d=cython.complex,

768 t1=cython.double,

769 t2=cython.double,

770 delta=cython.double,

771 delta_2=cython.double,

772 delta_3=cython.double,

773 a1=cython.complex,

774 b1=cython.complex,

775 c1=cython.complex,

776 d1=cython.complex,

777)

778def _splitCubicAtTC(a, b, c, d, *ts):

779 ts = list(ts)

780 ts.insert(0, 0.0)

781 ts.append(1.0)

782 for i in range(len(ts) - 1):

783 t1 = ts[i]

784 t2 = ts[i + 1]

785 delta = t2 - t1

786

787 delta_2 = delta * delta

788 delta_3 = delta * delta_2

789 t1_2 = t1 * t1

790 t1_3 = t1 * t1_2

791

792 # calc new a, b, c and d

793 a1 = a * delta_3

794 b1 = (3 * a * t1 + b) * delta_2

795 c1 = (2 * b * t1 + c + 3 * a * t1_2) * delta

796 d1 = a * t1_3 + b * t1_2 + c * t1 + d

797 pt1, pt2, pt3, pt4 = calcCubicPointsC(a1, b1, c1, d1)

798 yield (pt1, pt2, pt3, pt4)

799

800

801#

802# Equation solvers.

803#

804

805from math import sqrt, acos, cos, pi

806

807

808def solveQuadratic(a, b, c, sqrt=sqrt):

809 """Solve a quadratic equation.

810

811 Solves *a*x*x + b*x + c = 0* where a, b and c are real.

812

813 Args:

814 a: coefficient of *x²*

815 b: coefficient of *x*

816 c: constant term

817

818 Returns:

819 A list of roots. Note that the returned list is neither guaranteed to

820 be sorted nor to contain unique values!

821 """

822 if abs(a) < epsilon:

823 if abs(b) < epsilon:

824 # We have a non-equation; therefore, we have no valid solution

825 roots = []

826 else:

827 # We have a linear equation with 1 root.

828 roots = [-c / b]

829 else:

830 # We have a true quadratic equation. Apply the quadratic formula to find two roots.

831 DD = b * b - 4.0 * a * c

832 if DD >= 0.0:

833 rDD = sqrt(DD)

834 roots = [(-b + rDD) / 2.0 / a, (-b - rDD) / 2.0 / a]

835 else:

836 # complex roots, ignore

837 roots = []

838 return roots

839

840

841def solveCubic(a, b, c, d):

842 """Solve a cubic equation.

843

844 Solves *a*x*x*x + b*x*x + c*x + d = 0* where a, b, c and d are real.

845

846 Args:

847 a: coefficient of *x³*

848 b: coefficient of *x²*

849 c: coefficient of *x*

850 d: constant term

851

852 Returns:

853 A list of roots. Note that the returned list is neither guaranteed to

854 be sorted nor to contain unique values!

855

856 Examples::

857

858 >>> solveCubic(1, 1, -6, 0)

859 [-3.0, -0.0, 2.0]

860 >>> solveCubic(-10.0, -9.0, 48.0, -29.0)

861 [-2.9, 1.0, 1.0]

862 >>> solveCubic(-9.875, -9.0, 47.625, -28.75)

863 [-2.911392, 1.0, 1.0]

864 >>> solveCubic(1.0, -4.5, 6.75, -3.375)

865 [1.5, 1.5, 1.5]

866 >>> solveCubic(-12.0, 18.0, -9.0, 1.50023651123)

867 [0.5, 0.5, 0.5]

868 >>> solveCubic(

869 ... 9.0, 0.0, 0.0, -7.62939453125e-05

870 ... ) == [-0.0, -0.0, -0.0]

871 True

872 """

873 #

874 # adapted from:

875 # CUBIC.C - Solve a cubic polynomial

876 # public domain by Ross Cottrell

877 # found at: http://www.strangecreations.com/library/snippets/Cubic.C

878 #

879 if abs(a) < epsilon:

880 # don't just test for zero; for very small values of 'a' solveCubic()

881 # returns unreliable results, so we fall back to quad.

882 return solveQuadratic(b, c, d)

883 a = float(a)

884 a1 = b / a

885 a2 = c / a

886 a3 = d / a

887

888 Q = (a1 * a1 - 3.0 * a2) / 9.0

889 R = (2.0 * a1 * a1 * a1 - 9.0 * a1 * a2 + 27.0 * a3) / 54.0

890

891 R2 = R * R

892 Q3 = Q * Q * Q

893 R2 = 0 if R2 < epsilon else R2

894 Q3 = 0 if abs(Q3) < epsilon else Q3

895

896 R2_Q3 = R2 - Q3

897

898 if R2 == 0.0 and Q3 == 0.0:

899 x = round(-a1 / 3.0, epsilonDigits)

900 return [x, x, x]

901 elif R2_Q3 <= epsilon * 0.5:

902 # The epsilon * .5 above ensures that Q3 is not zero.

903 theta = acos(max(min(R / sqrt(Q3), 1.0), -1.0))

904 rQ2 = -2.0 * sqrt(Q)

905 a1_3 = a1 / 3.0

906 x0 = rQ2 * cos(theta / 3.0) - a1_3

907 x1 = rQ2 * cos((theta + 2.0 * pi) / 3.0) - a1_3

908 x2 = rQ2 * cos((theta + 4.0 * pi) / 3.0) - a1_3

909 x0, x1, x2 = sorted([x0, x1, x2])

910 # Merge roots that are close-enough

911 if x1 - x0 < epsilon and x2 - x1 < epsilon:

912 x0 = x1 = x2 = round((x0 + x1 + x2) / 3.0, epsilonDigits)

913 elif x1 - x0 < epsilon:

914 x0 = x1 = round((x0 + x1) / 2.0, epsilonDigits)

915 x2 = round(x2, epsilonDigits)

916 elif x2 - x1 < epsilon:

917 x0 = round(x0, epsilonDigits)

918 x1 = x2 = round((x1 + x2) / 2.0, epsilonDigits)

919 else:

920 x0 = round(x0, epsilonDigits)

921 x1 = round(x1, epsilonDigits)

922 x2 = round(x2, epsilonDigits)

923 return [x0, x1, x2]

924 else:

925 x = pow(sqrt(R2_Q3) + abs(R), 1 / 3.0)

926 x = x + Q / x

927 if R >= 0.0:

928 x = -x

929 x = round(x - a1 / 3.0, epsilonDigits)

930 return [x]

931

932

933#

934# Conversion routines for points to parameters and vice versa

935#

936

937

938def calcQuadraticParameters(pt1, pt2, pt3):

939 x2, y2 = pt2

940 x3, y3 = pt3

941 cx, cy = pt1

942 bx = (x2 - cx) * 2.0

943 by = (y2 - cy) * 2.0

944 ax = x3 - cx - bx

945 ay = y3 - cy - by

946 return (ax, ay), (bx, by), (cx, cy)

947

948

949def calcCubicParameters(pt1, pt2, pt3, pt4):

950 x2, y2 = pt2

951 x3, y3 = pt3

952 x4, y4 = pt4

953 dx, dy = pt1

954 cx = (x2 - dx) * 3.0

955 cy = (y2 - dy) * 3.0

956 bx = (x3 - x2) * 3.0 - cx

957 by = (y3 - y2) * 3.0 - cy

958 ax = x4 - dx - cx - bx

959 ay = y4 - dy - cy - by

960 return (ax, ay), (bx, by), (cx, cy), (dx, dy)

961

962

963@cython.cfunc

964@cython.inline

965@cython.locals(

966 pt1=cython.complex,

967 pt2=cython.complex,

968 pt3=cython.complex,

969 pt4=cython.complex,

970 a=cython.complex,

971 b=cython.complex,

972 c=cython.complex,

973)

974def calcCubicParametersC(pt1, pt2, pt3, pt4):

975 c = (pt2 - pt1) * 3.0

976 b = (pt3 - pt2) * 3.0 - c

977 a = pt4 - pt1 - c - b

978 return (a, b, c, pt1)

979

980

981def calcQuadraticPoints(a, b, c):

982 ax, ay = a

983 bx, by = b

984 cx, cy = c

985 x1 = cx

986 y1 = cy

987 x2 = (bx * 0.5) + cx

988 y2 = (by * 0.5) + cy

989 x3 = ax + bx + cx

990 y3 = ay + by + cy

991 return (x1, y1), (x2, y2), (x3, y3)

992

993

994def calcCubicPoints(a, b, c, d):

995 ax, ay = a

996 bx, by = b

997 cx, cy = c

998 dx, dy = d

999 x1 = dx

1000 y1 = dy

1001 x2 = (cx / 3.0) + dx

1002 y2 = (cy / 3.0) + dy

1003 x3 = (bx + cx) / 3.0 + x2

1004 y3 = (by + cy) / 3.0 + y2

1005 x4 = ax + dx + cx + bx

1006 y4 = ay + dy + cy + by

1007 return (x1, y1), (x2, y2), (x3, y3), (x4, y4)

1008

1009

1010@cython.cfunc

1011@cython.inline

1012@cython.locals(

1013 a=cython.complex,

1014 b=cython.complex,

1015 c=cython.complex,

1016 d=cython.complex,

1017 p2=cython.complex,

1018 p3=cython.complex,

1019 p4=cython.complex,

1020)

1021def calcCubicPointsC(a, b, c, d):

1022 p2 = c * (1 / 3) + d

1023 p3 = (b + c) * (1 / 3) + p2

1024 p4 = a + b + c + d

1025 return (d, p2, p3, p4)

1029# Point at time

1033def linePointAtT(pt1, pt2, t):

1034 """Finds the point at time `t` on a line.

1035

1036 Args:

1037 pt1, pt2: Coordinates of the line as 2D tuples.

1038 t: The time along the line.

1039

1040 Returns:

1041 A 2D tuple with the coordinates of the point.

1042 """

1043 return ((pt1[0] * (1 - t) + pt2[0] * t), (pt1[1] * (1 - t) + pt2[1] * t))

1044

1045

1046def quadraticPointAtT(pt1, pt2, pt3, t):

1047 """Finds the point at time `t` on a quadratic curve.

1048

1049 Args:

1050 pt1, pt2, pt3: Coordinates of the curve as 2D tuples.

1051 t: The time along the curve.

1052

1053 Returns:

1054 A 2D tuple with the coordinates of the point.

1055 """

1056 x = (1 - t) * (1 - t) * pt1[0] + 2 * (1 - t) * t * pt2[0] + t * t * pt3[0]

1057 y = (1 - t) * (1 - t) * pt1[1] + 2 * (1 - t) * t * pt2[1] + t * t * pt3[1]

1058 return (x, y)

1059

1060

1061def cubicPointAtT(pt1, pt2, pt3, pt4, t):

1062 """Finds the point at time `t` on a cubic curve.

1063

1064 Args:

1065 pt1, pt2, pt3, pt4: Coordinates of the curve as 2D tuples.

1066 t: The time along the curve.

1067

1068 Returns:

1069 A 2D tuple with the coordinates of the point.

1070 """

1071 t2 = t * t

1072 _1_t = 1 - t

1073 _1_t_2 = _1_t * _1_t

1074 x = (

1075 _1_t_2 * _1_t * pt1[0]

1076 + 3 * (_1_t_2 * t * pt2[0] + _1_t * t2 * pt3[0])

1077 + t2 * t * pt4[0]

1078 )

1079 y = (

1080 _1_t_2 * _1_t * pt1[1]

1081 + 3 * (_1_t_2 * t * pt2[1] + _1_t * t2 * pt3[1])

1082 + t2 * t * pt4[1]

1083 )

1084 return (x, y)

1085

1086

1087@cython.returns(cython.complex)

1088@cython.locals(

1089 t=cython.double,

1090 pt1=cython.complex,

1091 pt2=cython.complex,

1092 pt3=cython.complex,

1093 pt4=cython.complex,

1094)

1095@cython.locals(t2=cython.double, _1_t=cython.double, _1_t_2=cython.double)

1096def cubicPointAtTC(pt1, pt2, pt3, pt4, t):

1097 """Finds the point at time `t` on a cubic curve.

1098

1099 Args:

1100 pt1, pt2, pt3, pt4: Coordinates of the curve as complex numbers.

1101 t: The time along the curve.

1102

1103 Returns:

1104 A complex number with the coordinates of the point.

1105 """

1106 t2 = t * t

1107 _1_t = 1 - t

1108 _1_t_2 = _1_t * _1_t

1109 return _1_t_2 * _1_t * pt1 + 3 * (_1_t_2 * t * pt2 + _1_t * t2 * pt3) + t2 * t * pt4

1110

1111

1112def segmentPointAtT(seg, t):

1113 if len(seg) == 2:

1114 return linePointAtT(*seg, t)

1115 elif len(seg) == 3:

1116 return quadraticPointAtT(*seg, t)

1117 elif len(seg) == 4:

1118 return cubicPointAtT(*seg, t)

1119 raise ValueError("Unknown curve degree")

1120

1121

1122#

1123# Intersection finders

1124#

1125

1126

1127def _line_t_of_pt(s, e, pt):

1128 sx, sy = s

1129 ex, ey = e

1130 px, py = pt

1131 if abs(sx - ex) < epsilon and abs(sy - ey) < epsilon:

1132 # Line is a point!

1133 return -1

1134 # Use the largest

1135 if abs(sx - ex) > abs(sy - ey):

1136 return (px - sx) / (ex - sx)

1137 else:

1138 return (py - sy) / (ey - sy)

1139

1140

1141def _both_points_are_on_same_side_of_origin(a, b, origin):

1142 xDiff = (a[0] - origin[0]) * (b[0] - origin[0])

1143 yDiff = (a[1] - origin[1]) * (b[1] - origin[1])

1144 return not (xDiff <= 0.0 and yDiff <= 0.0)

1145

1146

1147def lineLineIntersections(s1, e1, s2, e2):

1148 """Finds intersections between two line segments.

1149

1150 Args:

1151 s1, e1: Coordinates of the first line as 2D tuples.

1152 s2, e2: Coordinates of the second line as 2D tuples.

1153

1154 Returns:

1155 A list of ``Intersection`` objects, each object having ``pt``, ``t1``

1156 and ``t2`` attributes containing the intersection point, time on first

1157 segment and time on second segment respectively.

1158

1159 Examples::

1160

1161 >>> a = lineLineIntersections( (310,389), (453, 222), (289, 251), (447, 367))

1162 >>> len(a)

1163 1

1164 >>> intersection = a[0]

1165 >>> intersection.pt

1166 (374.44882952482897, 313.73458370177315)

1167 >>> (intersection.t1, intersection.t2)

1168 (0.45069111555824465, 0.5408153767394238)

1169 """

1170 s1x, s1y = s1

1171 e1x, e1y = e1

1172 s2x, s2y = s2

1173 e2x, e2y = e2

1174 if (

1175 math.isclose(s2x, e2x) and math.isclose(s1x, e1x) and not math.isclose(s1x, s2x)

1176 ): # Parallel vertical

1177 return []

1178 if (

1179 math.isclose(s2y, e2y) and math.isclose(s1y, e1y) and not math.isclose(s1y, s2y)

1180 ): # Parallel horizontal

1181 return []

1182 if math.isclose(s2x, e2x) and math.isclose(s2y, e2y): # Line segment is tiny

1183 return []

1184 if math.isclose(s1x, e1x) and math.isclose(s1y, e1y): # Line segment is tiny

1185 return []

1186 if math.isclose(e1x, s1x):

1187 x = s1x

1188 slope34 = (e2y - s2y) / (e2x - s2x)

1189 y = slope34 * (x - s2x) + s2y

1190 pt = (x, y)

1191 return [

1192 Intersection(

1193 pt=pt, t1=_line_t_of_pt(s1, e1, pt), t2=_line_t_of_pt(s2, e2, pt)

1194 )

1195 ]

1196 if math.isclose(s2x, e2x):

1197 x = s2x

1198 slope12 = (e1y - s1y) / (e1x - s1x)

1199 y = slope12 * (x - s1x) + s1y

1200 pt = (x, y)

1201 return [

1202 Intersection(

1203 pt=pt, t1=_line_t_of_pt(s1, e1, pt), t2=_line_t_of_pt(s2, e2, pt)

1204 )

1205 ]

1206

1207 slope12 = (e1y - s1y) / (e1x - s1x)

1208 slope34 = (e2y - s2y) / (e2x - s2x)

1209 if math.isclose(slope12, slope34):

1210 return []

1211 x = (slope12 * s1x - s1y - slope34 * s2x + s2y) / (slope12 - slope34)

1212 y = slope12 * (x - s1x) + s1y

1213 pt = (x, y)

1214 if _both_points_are_on_same_side_of_origin(

1215 pt, e1, s1

1216 ) and _both_points_are_on_same_side_of_origin(pt, s2, e2):

1217 return [

1218 Intersection(

1219 pt=pt, t1=_line_t_of_pt(s1, e1, pt), t2=_line_t_of_pt(s2, e2, pt)

1220 )

1221 ]

1222 return []

1223

1224

1225def _alignment_transformation(segment):

1226 # Returns a transformation which aligns a segment horizontally at the

1227 # origin. Apply this transformation to curves and root-find to find

1228 # intersections with the segment.

1229 start = segment[0]

1230 end = segment[-1]

1231 angle = math.atan2(end[1] - start[1], end[0] - start[0])

1232 return Identity.rotate(-angle).translate(-start[0], -start[1])

1233

1234

1235def _curve_line_intersections_t(curve, line):

1236 aligned_curve = _alignment_transformation(line).transformPoints(curve)

1237 if len(curve) == 3:

1238 a, b, c = calcQuadraticParameters(*aligned_curve)

1239 intersections = solveQuadratic(a[1], b[1], c[1])

1240 elif len(curve) == 4:

1241 a, b, c, d = calcCubicParameters(*aligned_curve)

1242 intersections = solveCubic(a[1], b[1], c[1], d[1])

1243 else:

1244 raise ValueError("Unknown curve degree")

1245 return sorted(i for i in intersections if 0.0 <= i <= 1)

1246

1247

1248def curveLineIntersections(curve, line):

1249 """Finds intersections between a curve and a line.

1250

1251 Args:

1252 curve: List of coordinates of the curve segment as 2D tuples.

1253 line: List of coordinates of the line segment as 2D tuples.

1254

1255 Returns:

1256 A list of ``Intersection`` objects, each object having ``pt``, ``t1``

1257 and ``t2`` attributes containing the intersection point, time on first

1258 segment and time on second segment respectively.

1259

1260 Examples::

1261 >>> curve = [ (100, 240), (30, 60), (210, 230), (160, 30) ]

1262 >>> line = [ (25, 260), (230, 20) ]

1263 >>> intersections = curveLineIntersections(curve, line)

1264 >>> len(intersections)

1265 3

1266 >>> intersections[0].pt

1267 (84.9000930760723, 189.87306176459828)

1268 """

1269 if len(curve) == 3:

1270 pointFinder = quadraticPointAtT

1271 elif len(curve) == 4:

1272 pointFinder = cubicPointAtT

1273 else:

1274 raise ValueError("Unknown curve degree")

1275 intersections = []

1276 for t in _curve_line_intersections_t(curve, line):

1277 pt = pointFinder(*curve, t)

1278 # Back-project the point onto the line, to avoid problems with

1279 # numerical accuracy in the case of vertical and horizontal lines

1280 line_t = _line_t_of_pt(*line, pt)

1281 pt = linePointAtT(*line, line_t)

1282 intersections.append(Intersection(pt=pt, t1=t, t2=line_t))

1283 return intersections

1284

1285

1286def _curve_bounds(c):

1287 if len(c) == 3:

1288 return calcQuadraticBounds(*c)

1289 elif len(c) == 4:

1290 return calcCubicBounds(*c)

1291 raise ValueError("Unknown curve degree")

1292

1293

1294def _split_segment_at_t(c, t):

1295 if len(c) == 2:

1296 s, e = c

1297 midpoint = linePointAtT(s, e, t)

1298 return [(s, midpoint), (midpoint, e)]

1299 if len(c) == 3:

1300 return splitQuadraticAtT(*c, t)

1301 elif len(c) == 4:

1302 return splitCubicAtT(*c, t)

1303 raise ValueError("Unknown curve degree")

1304

1305

1306def _curve_curve_intersections_t(

1307 curve1, curve2, precision=1e-3, range1=None, range2=None

1308):

1309 bounds1 = _curve_bounds(curve1)

1310 bounds2 = _curve_bounds(curve2)

1311

1312 if not range1:

1313 range1 = (0.0, 1.0)

1314 if not range2:

1315 range2 = (0.0, 1.0)

1316

1317 # If bounds don't intersect, go home

1318 intersects, _ = sectRect(bounds1, bounds2)

1319 if not intersects:

1320 return []

1321

1322 def midpoint(r):

1323 return 0.5 * (r[0] + r[1])

1324

1325 # If they do overlap but they're tiny, approximate

1326 if rectArea(bounds1) < precision and rectArea(bounds2) < precision:

1327 return [(midpoint(range1), midpoint(range2))]

1328

1329 c11, c12 = _split_segment_at_t(curve1, 0.5)

1330 c11_range = (range1[0], midpoint(range1))

1331 c12_range = (midpoint(range1), range1[1])

1332

1333 c21, c22 = _split_segment_at_t(curve2, 0.5)

1334 c21_range = (range2[0], midpoint(range2))

1335 c22_range = (midpoint(range2), range2[1])

1336

1337 found = []

1338 found.extend(

1339 _curve_curve_intersections_t(

1340 c11, c21, precision, range1=c11_range, range2=c21_range

1341 )

1342 )

1343 found.extend(

1344 _curve_curve_intersections_t(

1345 c12, c21, precision, range1=c12_range, range2=c21_range

1346 )

1347 )

1348 found.extend(

1349 _curve_curve_intersections_t(

1350 c11, c22, precision, range1=c11_range, range2=c22_range

1351 )

1352 )

1353 found.extend(

1354 _curve_curve_intersections_t(

1355 c12, c22, precision, range1=c12_range, range2=c22_range

1356 )

1357 )

1358

1359 unique_key = lambda ts: (int(ts[0] / precision), int(ts[1] / precision))

1360 seen = set()

1361 unique_values = []

1362

1363 for ts in found:

1364 key = unique_key(ts)

1365 if key in seen:

1366 continue

1367 seen.add(key)

1368 unique_values.append(ts)

1369

1370 return unique_values

1371

1372

1373def curveCurveIntersections(curve1, curve2):

1374 """Finds intersections between a curve and a curve.

1375

1376 Args:

1377 curve1: List of coordinates of the first curve segment as 2D tuples.

1378 curve2: List of coordinates of the second curve segment as 2D tuples.

1379

1380 Returns:

1381 A list of ``Intersection`` objects, each object having ``pt``, ``t1``

1382 and ``t2`` attributes containing the intersection point, time on first

1383 segment and time on second segment respectively.

1384

1385 Examples::

1386 >>> curve1 = [ (10,100), (90,30), (40,140), (220,220) ]

1387 >>> curve2 = [ (5,150), (180,20), (80,250), (210,190) ]

1388 >>> intersections = curveCurveIntersections(curve1, curve2)

1389 >>> len(intersections)

1390 3

1391 >>> intersections[0].pt

1392 (81.7831487395506, 109.88904552375288)

1393 """

1394 intersection_ts = _curve_curve_intersections_t(curve1, curve2)

1395 return [

1396 Intersection(pt=segmentPointAtT(curve1, ts[0]), t1=ts[0], t2=ts[1])

1397 for ts in intersection_ts

1398 ]

1399

1400

1401def segmentSegmentIntersections(seg1, seg2):

1402 """Finds intersections between two segments.

1403

1404 Args:

1405 seg1: List of coordinates of the first segment as 2D tuples.

1406 seg2: List of coordinates of the second segment as 2D tuples.

1407

1408 Returns:

1409 A list of ``Intersection`` objects, each object having ``pt``, ``t1``

1410 and ``t2`` attributes containing the intersection point, time on first

1411 segment and time on second segment respectively.

1412

1413 Examples::

1414 >>> curve1 = [ (10,100), (90,30), (40,140), (220,220) ]

1415 >>> curve2 = [ (5,150), (180,20), (80,250), (210,190) ]

1416 >>> intersections = segmentSegmentIntersections(curve1, curve2)

1417 >>> len(intersections)

1418 3

1419 >>> intersections[0].pt

1420 (81.7831487395506, 109.88904552375288)

1421 >>> curve3 = [ (100, 240), (30, 60), (210, 230), (160, 30) ]

1422 >>> line = [ (25, 260), (230, 20) ]

1423 >>> intersections = segmentSegmentIntersections(curve3, line)

1424 >>> len(intersections)

1425 3

1426 >>> intersections[0].pt

1427 (84.9000930760723, 189.87306176459828)

1428

1429 """

1430 # Arrange by degree

1431 swapped = False

1432 if len(seg2) > len(seg1):

1433 seg2, seg1 = seg1, seg2

1434 swapped = True

1435 if len(seg1) > 2:

1436 if len(seg2) > 2:

1437 intersections = curveCurveIntersections(seg1, seg2)

1438 else:

1439 intersections = curveLineIntersections(seg1, seg2)

1440 elif len(seg1) == 2 and len(seg2) == 2:

1441 intersections = lineLineIntersections(*seg1, *seg2)

1442 else:

1443 raise ValueError("Couldn't work out which intersection function to use")

1444 if not swapped:

1445 return intersections

1446 return [Intersection(pt=i.pt, t1=i.t2, t2=i.t1) for i in intersections]

1447

1448

1449def _segmentrepr(obj):

1450 """

1451 >>> _segmentrepr([1, [2, 3], [], [[2, [3, 4], [0.1, 2.2]]]])

1452 '(1, (2, 3), (), ((2, (3, 4), (0.1, 2.2))))'

1453 """

1454 try:

1455 it = iter(obj)

1456 except TypeError:

1457 return "%g" % obj

1458 else:

1459 return "(%s)" % ", ".join(_segmentrepr(x) for x in it)

1460

1461

1462def printSegments(segments):

1463 """Helper for the doctests, displaying each segment in a list of

1464 segments on a single line as a tuple.

1465 """

1466 for segment in segments:

1467 print(_segmentrepr(segment))

1468

1469

1470if __name__ == "__main__":

1471 import sys

1472 import doctest

1473

1474 sys.exit(doctest.testmod().failed)