Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.11/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

12except (AttributeError, ImportError):

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

14 from fontTools.misc import cython

15COMPILED = cython.compiled

18EPSILON = 1e-9

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 + EPSILON >= 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.double)

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 split = _splitCubicAtT(a, b, c, d, *ts)

635

636 # the split impl can introduce floating point errors; we know the first

637 # segment should always start at pt1 and the last segment should end at pt4,

638 # so we set those values directly before returning.

639 split[0] = (pt1, *split[0][1:])

640 split[-1] = (*split[-1][:-1], pt4)

641 return split

642

643

644@cython.locals(

645 pt1=cython.complex,

646 pt2=cython.complex,

647 pt3=cython.complex,

648 pt4=cython.complex,

649 a=cython.complex,

650 b=cython.complex,

651 c=cython.complex,

652 d=cython.complex,

653)

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

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

656

657 Args:

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

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

660

661 Yields:

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

663 """

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

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

666

667

668@cython.returns(cython.complex)

669@cython.locals(

670 t=cython.double,

671 pt1=cython.complex,

672 pt2=cython.complex,

673 pt3=cython.complex,

674 pt4=cython.complex,

675 pointAtT=cython.complex,

676 off1=cython.complex,

677 off2=cython.complex,

678)

679@cython.locals(

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

681)

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

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

684

685 Args:

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

687 t: Position at which to split the curve.

688

689 Returns:

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

691 """

692 t2 = t * t

693 _1_t = 1 - t

694 _1_t_2 = _1_t * _1_t

695 _2_t_1_t = 2 * t * _1_t

696 pointAtT = (

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

698 )

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

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

701

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

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

704

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

706

707

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

709 ts = list(ts)

710 segments = []

711 ts.insert(0, 0.0)

712 ts.append(1.0)

713 ax, ay = a

714 bx, by = b

715 cx, cy = c

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

717 t1 = ts[i]

718 t2 = ts[i + 1]

719 delta = t2 - t1

720 # calc new a, b and c

721 delta_2 = delta * delta

722 a1x = ax * delta_2

723 a1y = ay * delta_2

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

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

726 t1_2 = t1 * t1

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

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

729

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

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

732 return segments

733

734

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

736 ts = list(ts)

737 ts.insert(0, 0.0)

738 ts.append(1.0)

739 segments = []

740 ax, ay = a

741 bx, by = b

742 cx, cy = c

743 dx, dy = d

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

745 t1 = ts[i]

746 t2 = ts[i + 1]

747 delta = t2 - t1

748

749 delta_2 = delta * delta

750 delta_3 = delta * delta_2

751 t1_2 = t1 * t1

752 t1_3 = t1 * t1_2

753

754 # calc new a, b, c and d

755 a1x = ax * delta_3

756 a1y = ay * delta_3

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

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

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

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

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

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

763 pt1, pt2, pt3, pt4 = calcCubicPoints(

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

765 )

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

767 return segments

768

769

770@cython.locals(

771 a=cython.complex,

772 b=cython.complex,

773 c=cython.complex,

774 d=cython.complex,

775 t1=cython.double,

776 t2=cython.double,

777 delta=cython.double,

778 delta_2=cython.double,

779 delta_3=cython.double,

780 a1=cython.complex,

781 b1=cython.complex,

782 c1=cython.complex,

783 d1=cython.complex,

784)

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

786 ts = list(ts)

787 ts.insert(0, 0.0)

788 ts.append(1.0)

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

790 t1 = ts[i]

791 t2 = ts[i + 1]

792 delta = t2 - t1

793

794 delta_2 = delta * delta

795 delta_3 = delta * delta_2

796 t1_2 = t1 * t1

797 t1_3 = t1 * t1_2

798

799 # calc new a, b, c and d

800 a1 = a * delta_3

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

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

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

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

805 yield (pt1, pt2, pt3, pt4)

806

807

808#

809# Equation solvers.

810#

811

812from math import sqrt, acos, cos, pi

813

814

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

816 """Solve a quadratic equation.

817

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

819

820 Args:

821 a: coefficient of *x²*

822 b: coefficient of *x*

823 c: constant term

824

825 Returns:

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

827 be sorted nor to contain unique values!

828 """

829 if abs(a) < epsilon:

830 if abs(b) < epsilon:

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

832 roots = []

833 else:

834 # We have a linear equation with 1 root.

835 roots = [-c / b]

836 else:

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

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

839 if DD >= 0.0:

840 rDD = sqrt(DD)

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

842 else:

843 # complex roots, ignore

844 roots = []

845 return roots

846

847

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

849 """Solve a cubic equation.

850

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

852

853 Args:

854 a: coefficient of *x³*

855 b: coefficient of *x²*

856 c: coefficient of *x*

857 d: constant term

858

859 Returns:

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

861 be sorted nor to contain unique values!

862

863 Examples::

864

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

866 [-3.0, -0.0, 2.0]

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

868 [-2.9, 1.0, 1.0]

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

870 [-2.911392, 1.0, 1.0]

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

872 [1.5, 1.5, 1.5]

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

874 [0.5, 0.5, 0.5]

875 >>> solveCubic(

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

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

878 True

879 """

880 #

881 # adapted from:

882 # CUBIC.C - Solve a cubic polynomial

883 # public domain by Ross Cottrell

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

885 #

886 if abs(a) < epsilon:

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

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

889 return solveQuadratic(b, c, d)

890 a = float(a)

891 a1 = b / a

892 a2 = c / a

893 a3 = d / a

894

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

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

897

898 R2 = R * R

899 Q3 = Q * Q * Q

900 R2 = 0 if R2 < epsilon else R2

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

902

903 R2_Q3 = R2 - Q3

904

905 if R2 == 0.0 and Q3 == 0.0:

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

907 return [x, x, x]

908 elif R2_Q3 <= epsilon * 0.5:

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

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

911 rQ2 = -2.0 * sqrt(Q)

912 a1_3 = a1 / 3.0

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

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

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

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

917 # Merge roots that are close-enough

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

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

920 elif x1 - x0 < epsilon:

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

922 x2 = round(x2, epsilonDigits)

923 elif x2 - x1 < epsilon:

924 x0 = round(x0, epsilonDigits)

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

926 else:

927 x0 = round(x0, epsilonDigits)

928 x1 = round(x1, epsilonDigits)

929 x2 = round(x2, epsilonDigits)

930 return [x0, x1, x2]

931 else:

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

933 x = x + Q / x

934 if R >= 0.0:

935 x = -x

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

937 return [x]

938

939

940#

941# Conversion routines for points to parameters and vice versa

942#

943

944

945def calcQuadraticParameters(pt1, pt2, pt3):

946 x2, y2 = pt2

947 x3, y3 = pt3

948 cx, cy = pt1

949 bx = (x2 - cx) * 2.0

950 by = (y2 - cy) * 2.0

951 ax = x3 - cx - bx

952 ay = y3 - cy - by

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

954

955

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

957 x2, y2 = pt2

958 x3, y3 = pt3

959 x4, y4 = pt4

960 dx, dy = pt1

961 cx = (x2 - dx) * 3.0

962 cy = (y2 - dy) * 3.0

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

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

965 ax = x4 - dx - cx - bx

966 ay = y4 - dy - cy - by

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

968

969

970@cython.cfunc

971@cython.inline

972@cython.locals(

973 pt1=cython.complex,

974 pt2=cython.complex,

975 pt3=cython.complex,

976 pt4=cython.complex,

977 a=cython.complex,

978 b=cython.complex,

979 c=cython.complex,

980)

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

982 c = (pt2 - pt1) * 3.0

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

984 a = pt4 - pt1 - c - b

985 return (a, b, c, pt1)

986

987

988def calcQuadraticPoints(a, b, c):

989 ax, ay = a

990 bx, by = b

991 cx, cy = c

992 x1 = cx

993 y1 = cy

994 x2 = (bx * 0.5) + cx

995 y2 = (by * 0.5) + cy

996 x3 = ax + bx + cx

997 y3 = ay + by + cy

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

999

1000

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

1002 ax, ay = a

1003 bx, by = b

1004 cx, cy = c

1005 dx, dy = d

1006 x1 = dx

1007 y1 = dy

1008 x2 = (cx / 3.0) + dx

1009 y2 = (cy / 3.0) + dy

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

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

1012 x4 = ax + dx + cx + bx

1013 y4 = ay + dy + cy + by

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

1015

1016

1017@cython.cfunc

1018@cython.inline

1019@cython.locals(

1020 a=cython.complex,

1021 b=cython.complex,

1022 c=cython.complex,

1023 d=cython.complex,

1024 p2=cython.complex,

1025 p3=cython.complex,

1026 p4=cython.complex,

1027)

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

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

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

1031 p4 = a + b + c + d

1032 return (d, p2, p3, p4)

1036# Point at time

1040def linePointAtT(pt1, pt2, t):

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

1042

1043 Args:

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

1045 t: The time along the line.

1046

1047 Returns:

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

1049 """

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

1051

1052

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

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

1055

1056 Args:

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

1058 t: The time along the curve.

1059

1060 Returns:

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

1062 """

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

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

1065 return (x, y)

1066

1067

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

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

1070

1071 Args:

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

1073 t: The time along the curve.

1074

1075 Returns:

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

1077 """

1078 t2 = t * t

1079 _1_t = 1 - t

1080 _1_t_2 = _1_t * _1_t

1081 x = (

1082 _1_t_2 * _1_t * pt1[0]

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

1084 + t2 * t * pt4[0]

1085 )

1086 y = (

1087 _1_t_2 * _1_t * pt1[1]

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

1089 + t2 * t * pt4[1]

1090 )

1091 return (x, y)

1092

1093

1094@cython.returns(cython.complex)

1095@cython.locals(

1096 t=cython.double,

1097 pt1=cython.complex,

1098 pt2=cython.complex,

1099 pt3=cython.complex,

1100 pt4=cython.complex,

1101)

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

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

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

1105

1106 Args:

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

1108 t: The time along the curve.

1109

1110 Returns:

1111 A complex number with the coordinates of the point.

1112 """

1113 t2 = t * t

1114 _1_t = 1 - t

1115 _1_t_2 = _1_t * _1_t

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

1117

1118

1119def segmentPointAtT(seg, t):

1120 if len(seg) == 2:

1121 return linePointAtT(*seg, t)

1122 elif len(seg) == 3:

1123 return quadraticPointAtT(*seg, t)

1124 elif len(seg) == 4:

1125 return cubicPointAtT(*seg, t)

1126 raise ValueError("Unknown curve degree")

1127

1128

1129#

1130# Intersection finders

1131#

1132

1133

1134def _line_t_of_pt(s, e, pt):

1135 sx, sy = s

1136 ex, ey = e

1137 px, py = pt

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

1139 # Line is a point!

1140 return -1

1141 # Use the largest

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

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

1144 else:

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

1146

1147

1148def _both_points_are_on_same_side_of_origin(a, b, origin):

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

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

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

1152

1153

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

1155 """Finds intersections between two line segments.

1156

1157 Args:

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

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

1160

1161 Returns:

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

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

1164 segment and time on second segment respectively.

1165

1166 Examples::

1167

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

1169 >>> len(a)

1170 1

1171 >>> intersection = a[0]

1172 >>> intersection.pt

1173 (374.44882952482897, 313.73458370177315)

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

1175 (0.45069111555824465, 0.5408153767394238)

1176 """

1177 s1x, s1y = s1

1178 e1x, e1y = e1

1179 s2x, s2y = s2

1180 e2x, e2y = e2

1181 if (

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

1183 ): # Parallel vertical

1184 return []

1185 if (

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

1187 ): # Parallel horizontal

1188 return []

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

1190 return []

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

1192 return []

1193 if math.isclose(e1x, s1x):

1194 x = s1x

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

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

1197 pt = (x, y)

1198 return [

1199 Intersection(

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

1201 )

1202 ]

1203 if math.isclose(s2x, e2x):

1204 x = s2x

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

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

1207 pt = (x, y)

1208 return [

1209 Intersection(

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

1211 )

1212 ]

1213

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

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

1216 if math.isclose(slope12, slope34):

1217 return []

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

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

1220 pt = (x, y)

1221 if _both_points_are_on_same_side_of_origin(

1222 pt, e1, s1

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

1224 return [

1225 Intersection(

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

1227 )

1228 ]

1229 return []

1230

1231

1232def _alignment_transformation(segment):

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

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

1235 # intersections with the segment.

1236 start = segment[0]

1237 end = segment[-1]

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

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

1240

1241

1242def _curve_line_intersections_t(curve, line):

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

1244 if len(curve) == 3:

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

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

1247 elif len(curve) == 4:

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

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

1250 else:

1251 raise ValueError("Unknown curve degree")

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

1253

1254

1255def curveLineIntersections(curve, line):

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

1257

1258 Args:

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

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

1261

1262 Returns:

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

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

1265 segment and time on second segment respectively.

1266

1267 Examples::

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

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

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

1271 >>> len(intersections)

1272 3

1273 >>> intersections[0].pt

1274 (84.9000930760723, 189.87306176459828)

1275 """

1276 if len(curve) == 3:

1277 pointFinder = quadraticPointAtT

1278 elif len(curve) == 4:

1279 pointFinder = cubicPointAtT

1280 else:

1281 raise ValueError("Unknown curve degree")

1282 intersections = []

1283 for t in _curve_line_intersections_t(curve, line):

1284 pt = pointFinder(*curve, t)

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

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

1287 line_t = _line_t_of_pt(*line, pt)

1288 pt = linePointAtT(*line, line_t)

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

1290 return intersections

1291

1292

1293def _curve_bounds(c):

1294 if len(c) == 3:

1295 return calcQuadraticBounds(*c)

1296 elif len(c) == 4:

1297 return calcCubicBounds(*c)

1298 raise ValueError("Unknown curve degree")

1299

1300

1301def _split_segment_at_t(c, t):

1302 if len(c) == 2:

1303 s, e = c

1304 midpoint = linePointAtT(s, e, t)

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

1306 if len(c) == 3:

1307 return splitQuadraticAtT(*c, t)

1308 elif len(c) == 4:

1309 return splitCubicAtT(*c, t)

1310 raise ValueError("Unknown curve degree")

1311

1312

1313def _curve_curve_intersections_t(

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

1315):

1316 bounds1 = _curve_bounds(curve1)

1317 bounds2 = _curve_bounds(curve2)

1318

1319 if not range1:

1320 range1 = (0.0, 1.0)

1321 if not range2:

1322 range2 = (0.0, 1.0)

1323

1324 # If bounds don't intersect, go home

1325 intersects, _ = sectRect(bounds1, bounds2)

1326 if not intersects:

1327 return []

1328

1329 def midpoint(r):

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

1331

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

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

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

1335

1336 c11, c12 = _split_segment_at_t(curve1, 0.5)

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

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

1339

1340 c21, c22 = _split_segment_at_t(curve2, 0.5)

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

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

1343

1344 found = []

1345 found.extend(

1346 _curve_curve_intersections_t(

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

1348 )

1349 )

1350 found.extend(

1351 _curve_curve_intersections_t(

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

1353 )

1354 )

1355 found.extend(

1356 _curve_curve_intersections_t(

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

1358 )

1359 )

1360 found.extend(

1361 _curve_curve_intersections_t(

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

1363 )

1364 )

1365

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

1367 seen = set()

1368 unique_values = []

1369

1370 for ts in found:

1371 key = unique_key(ts)

1372 if key in seen:

1373 continue

1374 seen.add(key)

1375 unique_values.append(ts)

1376

1377 return unique_values

1378

1379

1380def _is_linelike(segment):

1381 maybeline = _alignment_transformation(segment).transformPoints(segment)

1382 return all(math.isclose(p[1], 0.0) for p in maybeline)

1383

1384

1385def curveCurveIntersections(curve1, curve2):

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

1387

1388 Args:

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

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

1391

1392 Returns:

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

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

1395 segment and time on second segment respectively.

1396

1397 Examples::

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

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

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

1401 >>> len(intersections)

1402 3

1403 >>> intersections[0].pt

1404 (81.7831487395506, 109.88904552375288)

1405 """

1406 if _is_linelike(curve1):

1407 line1 = curve1[0], curve1[-1]

1408 if _is_linelike(curve2):

1409 line2 = curve2[0], curve2[-1]

1410 return lineLineIntersections(*line1, *line2)

1411 else:

1412 return curveLineIntersections(curve2, line1)

1413 elif _is_linelike(curve2):

1414 line2 = curve2[0], curve2[-1]

1415 return curveLineIntersections(curve1, line2)

1416

1417 intersection_ts = _curve_curve_intersections_t(curve1, curve2)

1418 return [

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

1420 for ts in intersection_ts

1421 ]

1422

1423

1424def segmentSegmentIntersections(seg1, seg2):

1425 """Finds intersections between two segments.

1426

1427 Args:

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

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

1430

1431 Returns:

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

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

1434 segment and time on second segment respectively.

1435

1436 Examples::

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

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

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

1440 >>> len(intersections)

1441 3

1442 >>> intersections[0].pt

1443 (81.7831487395506, 109.88904552375288)

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

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

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

1447 >>> len(intersections)

1448 3

1449 >>> intersections[0].pt

1450 (84.9000930760723, 189.87306176459828)

1451

1452 """

1453 # Arrange by degree

1454 swapped = False

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

1456 seg2, seg1 = seg1, seg2

1457 swapped = True

1458 if len(seg1) > 2:

1459 if len(seg2) > 2:

1460 intersections = curveCurveIntersections(seg1, seg2)

1461 else:

1462 intersections = curveLineIntersections(seg1, seg2)

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

1464 intersections = lineLineIntersections(*seg1, *seg2)

1465 else:

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

1467 if not swapped:

1468 return intersections

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

1470

1471

1472def _segmentrepr(obj):

1473 """

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

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

1476 """

1477 try:

1478 it = iter(obj)

1479 except TypeError:

1480 return "%g" % obj

1481 else:

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

1483

1484

1485def printSegments(segments):

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

1487 segments on a single line as a tuple.

1488 """

1489 for segment in segments:

1490 print(_segmentrepr(segment))

1491

1492

1493if __name__ == "__main__":

1494 import sys

1495 import doctest

1496

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