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

21EPSILON = 1e-9

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

27__all__ = [

28 "approximateCubicArcLength",

29 "approximateCubicArcLengthC",

30 "approximateQuadraticArcLength",

31 "approximateQuadraticArcLengthC",

32 "calcCubicArcLength",

33 "calcCubicArcLengthC",

34 "calcQuadraticArcLength",

35 "calcQuadraticArcLengthC",

36 "calcCubicBounds",

37 "calcQuadraticBounds",

38 "splitLine",

39 "splitQuadratic",

40 "splitCubic",

41 "splitQuadraticAtT",

42 "splitCubicAtT",

43 "splitCubicAtTC",

44 "splitCubicIntoTwoAtTC",

45 "solveQuadratic",

46 "solveCubic",

47 "quadraticPointAtT",

48 "cubicPointAtT",

49 "cubicPointAtTC",

50 "linePointAtT",

51 "segmentPointAtT",

52 "lineLineIntersections",

53 "curveLineIntersections",

54 "curveCurveIntersections",

55 "segmentSegmentIntersections",

56]

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

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

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

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

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

66 Args:

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

68 tolerance: Controls the precision of the calcuation.

70 Returns:

71 Arc length value.

72 """

73 return calcCubicArcLengthC(

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

75 )

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

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

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

81 return (

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

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

84 )

87@cython.returns(cython.double)

88@cython.locals(

89 p0=cython.complex,

90 p1=cython.complex,

91 p2=cython.complex,

92 p3=cython.complex,

93)

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

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

96 arch = abs(p0 - p3)

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

98 if arch * mult + EPSILON >= box:

99 return (arch + box) * 0.5

100 else:

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

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

103 mult, *two

104 )

105

106

107@cython.returns(cython.double)

108@cython.locals(

109 pt1=cython.complex,

110 pt2=cython.complex,

111 pt3=cython.complex,

112 pt4=cython.complex,

113)

114@cython.locals(

115 tolerance=cython.double,

116 mult=cython.double,

117)

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

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

120

121 Args:

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

123 tolerance: Controls the precision of the calcuation.

124

125 Returns:

126 Arc length value.

127 """

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

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

130

131

132epsilonDigits = 6

133epsilon = 1e-10

134

135

136@cython.cfunc

137@cython.inline

138@cython.returns(cython.double)

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

140def _dot(v1, v2):

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

142

143

144@cython.cfunc

145@cython.inline

146@cython.returns(cython.double)

147@cython.locals(x=cython.double)

148def _intSecAtan(x):

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

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

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

152

153

154def calcQuadraticArcLength(pt1, pt2, pt3):

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

156

157 Args:

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

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

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

161

162 Returns:

163 Arc length value.

164

165 Example::

166

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

168 0.0

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

170 80.0

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

172 80.0

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

174 107.70329614269008

175 >>> calcQuadraticArcLength((0, 0), (0, 100), (100, 0))

176 154.02976155645263

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

178 120.21581243984076

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

180 102.53273816445825

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

182 66.66666666666667

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

184 40.0

185 """

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

187

188

189@cython.returns(cython.double)

190@cython.locals(

191 pt1=cython.complex,

192 pt2=cython.complex,

193 pt3=cython.complex,

194 d0=cython.complex,

195 d1=cython.complex,

196 d=cython.complex,

197 n=cython.complex,

198)

199@cython.locals(

200 scale=cython.double,

201 origDist=cython.double,

202 a=cython.double,

203 b=cython.double,

204 x0=cython.double,

205 x1=cython.double,

206 Len=cython.double,

207)

208def calcQuadraticArcLengthC(pt1, pt2, pt3):

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

210

211 Args:

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

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

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

215

216 Returns:

217 Arc length value.

218 """

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

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

221 d0 = pt2 - pt1

222 d1 = pt3 - pt2

223 d = d1 - d0

224 n = d * 1j

225 scale = abs(n)

226 if scale == 0.0:

227 return abs(pt3 - pt1)

228 origDist = _dot(n, d0)

229 if abs(origDist) < epsilon:

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

231 return abs(pt3 - pt1)

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

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

234 x0 = _dot(d, d0) / origDist

235 x1 = _dot(d, d1) / origDist

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

237 return Len

238

239

240def approximateQuadraticArcLength(pt1, pt2, pt3):

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

242

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

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

245

246 Args:

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

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

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

250

251 Returns:

252 Approximate arc length value.

253 """

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

255

256

257@cython.returns(cython.double)

258@cython.locals(

259 pt1=cython.complex,

260 pt2=cython.complex,

261 pt3=cython.complex,

262)

263@cython.locals(

264 v0=cython.double,

265 v1=cython.double,

266 v2=cython.double,

267)

268def approximateQuadraticArcLengthC(pt1, pt2, pt3):

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

270

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

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

273

274 Args:

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

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

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

278

279 Returns:

280 Approximate arc length value.

281 """

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

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

284 # approximation of it.

285 #

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

287

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

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

290 v0 = abs(

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

292 )

293 v1 = abs(pt3 - pt1) * 0.4444444444444444

294 v2 = abs(

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

296 )

297

298 return v0 + v1 + v2

299

300

301def calcQuadraticBounds(pt1, pt2, pt3):

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

303

304 Args:

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

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

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

308

309 Returns:

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

311

312 Example::

313

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

315 (0, 0, 100, 50.0)

316 >>> calcQuadraticBounds((0, 0), (100, 0), (100, 100))

317 (0.0, 0.0, 100, 100)

318 """

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

320 ax2 = ax * 2.0

321 ay2 = ay * 2.0

322 roots = []

323 if ax2 != 0:

324 roots.append(-bx / ax2)

325 if ay2 != 0:

326 roots.append(-by / ay2)

327 points = [

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

329 for t in roots

330 if 0 <= t < 1

331 ] + [pt1, pt3]

332 return calcBounds(points)

333

334

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

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

337

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

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

340

341 Args:

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

343

344 Returns:

345 Arc length value.

346

347 Example::

348

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

350 190.04332968932817

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

352 154.8852074945903

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

354 149.99999999999991

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

356 136.9267662156362

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

358 154.80848416537057

359 """

360 return approximateCubicArcLengthC(

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

362 )

363

364

365@cython.returns(cython.double)

366@cython.locals(

367 pt1=cython.complex,

368 pt2=cython.complex,

369 pt3=cython.complex,

370 pt4=cython.complex,

371)

372@cython.locals(

373 v0=cython.double,

374 v1=cython.double,

375 v2=cython.double,

376 v3=cython.double,

377 v4=cython.double,

378)

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

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

381

382 Args:

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

384

385 Returns:

386 Arc length value.

387 """

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

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

390 # approximation of it.

391 #

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

393

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

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

396 v0 = abs(pt2 - pt1) * 0.15

397 v1 = abs(

398 -0.558983582205757 * pt1

399 + 0.325650248872424 * pt2

400 + 0.208983582205757 * pt3

401 + 0.024349751127576 * pt4

402 )

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

404 v3 = abs(

405 -0.024349751127576 * pt1

406 - 0.208983582205757 * pt2

407 - 0.325650248872424 * pt3

408 + 0.558983582205757 * pt4

409 )

410 v4 = abs(pt4 - pt3) * 0.15

411

412 return v0 + v1 + v2 + v3 + v4

413

414

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

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

417

418 Args:

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

420

421 Returns:

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

423

424 Example::

425

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

427 (0, 0, 100, 75.0)

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

429 (0.0, 0.0, 100, 100)

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

431 35.566243 0.000000 64.433757 75.000000

432 """

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

434 # calc first derivative

435 ax3 = ax * 3.0

436 ay3 = ay * 3.0

437 bx2 = bx * 2.0

438 by2 = by * 2.0

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

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

441 roots = xRoots + yRoots

442

443 points = [

444 (

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

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

447 )

448 for t in roots

449 ] + [pt1, pt4]

450 return calcBounds(points)

451

452

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

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

455

456 Args:

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

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

459 where: Position at which to split the line.

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

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

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

463

464 Returns:

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

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

467 line.

468

469 Example::

470

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

472 ((0, 0), (50, 50))

473 ((50, 50), (100, 100))

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

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

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

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

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

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

480 ((0, 0), (0, 0))

481 ((0, 0), (100, 100))

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

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

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

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

486 ((0, 100), (0, 50))

487 ((0, 50), (0, 0))

488 """

489 pt1x, pt1y = pt1

490 pt2x, pt2y = pt2

491

492 ax = pt2x - pt1x

493 ay = pt2y - pt1y

494

495 bx = pt1x

496 by = pt1y

497

498 a = (ax, ay)[isHorizontal]

499

500 if a == 0:

501 return [(pt1, pt2)]

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

503 if 0 <= t < 1:

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

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

506 else:

507 return [(pt1, pt2)]

508

509

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

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

512

513 Args:

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

515 where: Position at which to split the curve.

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

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

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

519

520 Returns:

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

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

523 curve.

524

525 Example::

526

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

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

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

530 ((0, 0), (25, 50), (50, 50))

531 ((50, 50), (75, 50), (100, 0))

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

533 ((0, 0), (12.5, 25), (25, 37.5))

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

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

536 ((0, 0), (7.32233, 14.6447), (14.6447, 25))

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

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

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

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

541 ((0, 0), (25, 50), (50, 50))

542 ((50, 50), (50, 50), (50, 50))

543 ((50, 50), (75, 50), (100, 0))

544 """

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

546 solutions = solveQuadratic(

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

548 )

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

550 if not solutions:

551 return [(pt1, pt2, pt3)]

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

553

554

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

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

557

558 Args:

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

560 where: Position at which to split the curve.

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

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

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

564

565 Returns:

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

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

568 curve.

569

570 Example::

571

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

573 ((0, 0), (25, 100), (75, 100), (100, 0))

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

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

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

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

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

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

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

581 """

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

583 solutions = solveCubic(

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

585 )

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

587 if not solutions:

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

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

590

591

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

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

594

595 Args:

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

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

598

599 Returns:

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

601

602 Examples::

603

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

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

606 ((50, 50), (75, 50), (100, 0))

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

608 ((0, 0), (25, 50), (50, 50))

609 ((50, 50), (62.5, 50), (75, 37.5))

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

611 """

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

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

614

615

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

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

618

619 Args:

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

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

622

623 Returns:

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

625

626 Examples::

627

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

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

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

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

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

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

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

635 """

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

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

638

639

640@cython.locals(

641 pt1=cython.complex,

642 pt2=cython.complex,

643 pt3=cython.complex,

644 pt4=cython.complex,

645 a=cython.complex,

646 b=cython.complex,

647 c=cython.complex,

648 d=cython.complex,

649)

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

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

652

653 Args:

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

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

656

657 Yields:

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

659 """

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

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

662

663

664@cython.returns(cython.complex)

665@cython.locals(

666 t=cython.double,

667 pt1=cython.complex,

668 pt2=cython.complex,

669 pt3=cython.complex,

670 pt4=cython.complex,

671 pointAtT=cython.complex,

672 off1=cython.complex,

673 off2=cython.complex,

674)

675@cython.locals(

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

677)

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

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

680

681 Args:

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

683 t: Position at which to split the curve.

684

685 Returns:

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

687 """

688 t2 = t * t

689 _1_t = 1 - t

690 _1_t_2 = _1_t * _1_t

691 _2_t_1_t = 2 * t * _1_t

692 pointAtT = (

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

694 )

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

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

697

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

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

700

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

702

703

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

705 ts = list(ts)

706 segments = []

707 ts.insert(0, 0.0)

708 ts.append(1.0)

709 ax, ay = a

710 bx, by = b

711 cx, cy = c

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

713 t1 = ts[i]

714 t2 = ts[i + 1]

715 delta = t2 - t1

716 # calc new a, b and c

717 delta_2 = delta * delta

718 a1x = ax * delta_2

719 a1y = ay * delta_2

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

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

722 t1_2 = t1 * t1

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

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

725

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

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

728 return segments

729

730

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

732 ts = list(ts)

733 ts.insert(0, 0.0)

734 ts.append(1.0)

735 segments = []

736 ax, ay = a

737 bx, by = b

738 cx, cy = c

739 dx, dy = d

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

741 t1 = ts[i]

742 t2 = ts[i + 1]

743 delta = t2 - t1

744

745 delta_2 = delta * delta

746 delta_3 = delta * delta_2

747 t1_2 = t1 * t1

748 t1_3 = t1 * t1_2

749

750 # calc new a, b, c and d

751 a1x = ax * delta_3

752 a1y = ay * delta_3

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

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

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

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

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

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

759 pt1, pt2, pt3, pt4 = calcCubicPoints(

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

761 )

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

763 return segments

764

765

766@cython.locals(

767 a=cython.complex,

768 b=cython.complex,

769 c=cython.complex,

770 d=cython.complex,

771 t1=cython.double,

772 t2=cython.double,

773 delta=cython.double,

774 delta_2=cython.double,

775 delta_3=cython.double,

776 a1=cython.complex,

777 b1=cython.complex,

778 c1=cython.complex,

779 d1=cython.complex,

780)

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

782 ts = list(ts)

783 ts.insert(0, 0.0)

784 ts.append(1.0)

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

786 t1 = ts[i]

787 t2 = ts[i + 1]

788 delta = t2 - t1

789

790 delta_2 = delta * delta

791 delta_3 = delta * delta_2

792 t1_2 = t1 * t1

793 t1_3 = t1 * t1_2

794

795 # calc new a, b, c and d

796 a1 = a * delta_3

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

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

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

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

801 yield (pt1, pt2, pt3, pt4)

802

803

804#

805# Equation solvers.

806#

807

808from math import sqrt, acos, cos, pi

809

810

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

812 """Solve a quadratic equation.

813

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

815

816 Args:

817 a: coefficient of *x²*

818 b: coefficient of *x*

819 c: constant term

820

821 Returns:

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

823 be sorted nor to contain unique values!

824 """

825 if abs(a) < epsilon:

826 if abs(b) < epsilon:

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

828 roots = []

829 else:

830 # We have a linear equation with 1 root.

831 roots = [-c / b]

832 else:

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

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

835 if DD >= 0.0:

836 rDD = sqrt(DD)

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

838 else:

839 # complex roots, ignore

840 roots = []

841 return roots

842

843

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

845 """Solve a cubic equation.

846

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

848

849 Args:

850 a: coefficient of *x³*

851 b: coefficient of *x²*

852 c: coefficient of *x*

853 d: constant term

854

855 Returns:

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

857 be sorted nor to contain unique values!

858

859 Examples::

860

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

862 [-3.0, -0.0, 2.0]

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

864 [-2.9, 1.0, 1.0]

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

866 [-2.911392, 1.0, 1.0]

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

868 [1.5, 1.5, 1.5]

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

870 [0.5, 0.5, 0.5]

871 >>> solveCubic(

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

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

874 True

875 """

876 #

877 # adapted from:

878 # CUBIC.C - Solve a cubic polynomial

879 # public domain by Ross Cottrell

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

881 #

882 if abs(a) < epsilon:

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

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

885 return solveQuadratic(b, c, d)

886 a = float(a)

887 a1 = b / a

888 a2 = c / a

889 a3 = d / a

890

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

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

893

894 R2 = R * R

895 Q3 = Q * Q * Q

896 R2 = 0 if R2 < epsilon else R2

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

898

899 R2_Q3 = R2 - Q3

900

901 if R2 == 0.0 and Q3 == 0.0:

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

903 return [x, x, x]

904 elif R2_Q3 <= epsilon * 0.5:

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

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

907 rQ2 = -2.0 * sqrt(Q)

908 a1_3 = a1 / 3.0

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

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

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

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

913 # Merge roots that are close-enough

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

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

916 elif x1 - x0 < epsilon:

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

918 x2 = round(x2, epsilonDigits)

919 elif x2 - x1 < epsilon:

920 x0 = round(x0, epsilonDigits)

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

922 else:

923 x0 = round(x0, epsilonDigits)

924 x1 = round(x1, epsilonDigits)

925 x2 = round(x2, epsilonDigits)

926 return [x0, x1, x2]

927 else:

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

929 x = x + Q / x

930 if R >= 0.0:

931 x = -x

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

933 return [x]

934

935

936#

937# Conversion routines for points to parameters and vice versa

938#

939

940

941def calcQuadraticParameters(pt1, pt2, pt3):

942 x2, y2 = pt2

943 x3, y3 = pt3

944 cx, cy = pt1

945 bx = (x2 - cx) * 2.0

946 by = (y2 - cy) * 2.0

947 ax = x3 - cx - bx

948 ay = y3 - cy - by

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

950

951

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

953 x2, y2 = pt2

954 x3, y3 = pt3

955 x4, y4 = pt4

956 dx, dy = pt1

957 cx = (x2 - dx) * 3.0

958 cy = (y2 - dy) * 3.0

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

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

961 ax = x4 - dx - cx - bx

962 ay = y4 - dy - cy - by

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

964

965

966@cython.cfunc

967@cython.inline

968@cython.locals(

969 pt1=cython.complex,

970 pt2=cython.complex,

971 pt3=cython.complex,

972 pt4=cython.complex,

973 a=cython.complex,

974 b=cython.complex,

975 c=cython.complex,

976)

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

978 c = (pt2 - pt1) * 3.0

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

980 a = pt4 - pt1 - c - b

981 return (a, b, c, pt1)

982

983

984def calcQuadraticPoints(a, b, c):

985 ax, ay = a

986 bx, by = b

987 cx, cy = c

988 x1 = cx

989 y1 = cy

990 x2 = (bx * 0.5) + cx

991 y2 = (by * 0.5) + cy

992 x3 = ax + bx + cx

993 y3 = ay + by + cy

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

995

996

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

998 ax, ay = a

999 bx, by = b

1000 cx, cy = c

1001 dx, dy = d

1002 x1 = dx

1003 y1 = dy

1004 x2 = (cx / 3.0) + dx

1005 y2 = (cy / 3.0) + dy

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

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

1008 x4 = ax + dx + cx + bx

1009 y4 = ay + dy + cy + by

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

1011

1012

1013@cython.cfunc

1014@cython.inline

1015@cython.locals(

1016 a=cython.complex,

1017 b=cython.complex,

1018 c=cython.complex,

1019 d=cython.complex,

1020 p2=cython.complex,

1021 p3=cython.complex,

1022 p4=cython.complex,

1023)

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

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

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

1027 p4 = a + b + c + d

1028 return (d, p2, p3, p4)

1032# Point at time

1036def linePointAtT(pt1, pt2, t):

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

1038

1039 Args:

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

1041 t: The time along the line.

1042

1043 Returns:

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

1045 """

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

1047

1048

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

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

1051

1052 Args:

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

1054 t: The time along the curve.

1055

1056 Returns:

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

1058 """

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

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

1061 return (x, y)

1062

1063

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

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

1066

1067 Args:

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

1069 t: The time along the curve.

1070

1071 Returns:

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

1073 """

1074 t2 = t * t

1075 _1_t = 1 - t

1076 _1_t_2 = _1_t * _1_t

1077 x = (

1078 _1_t_2 * _1_t * pt1[0]

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

1080 + t2 * t * pt4[0]

1081 )

1082 y = (

1083 _1_t_2 * _1_t * pt1[1]

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

1085 + t2 * t * pt4[1]

1086 )

1087 return (x, y)

1088

1089

1090@cython.returns(cython.complex)

1091@cython.locals(

1092 t=cython.double,

1093 pt1=cython.complex,

1094 pt2=cython.complex,

1095 pt3=cython.complex,

1096 pt4=cython.complex,

1097)

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

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

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

1101

1102 Args:

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

1104 t: The time along the curve.

1105

1106 Returns:

1107 A complex number with the coordinates of the point.

1108 """

1109 t2 = t * t

1110 _1_t = 1 - t

1111 _1_t_2 = _1_t * _1_t

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

1113

1114

1115def segmentPointAtT(seg, t):

1116 if len(seg) == 2:

1117 return linePointAtT(*seg, t)

1118 elif len(seg) == 3:

1119 return quadraticPointAtT(*seg, t)

1120 elif len(seg) == 4:

1121 return cubicPointAtT(*seg, t)

1122 raise ValueError("Unknown curve degree")

1123

1124

1125#

1126# Intersection finders

1127#

1128

1129

1130def _line_t_of_pt(s, e, pt):

1131 sx, sy = s

1132 ex, ey = e

1133 px, py = pt

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

1135 # Line is a point!

1136 return -1

1137 # Use the largest

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

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

1140 else:

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

1142

1143

1144def _both_points_are_on_same_side_of_origin(a, b, origin):

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

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

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

1148

1149

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

1151 """Finds intersections between two line segments.

1152

1153 Args:

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

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

1156

1157 Returns:

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

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

1160 segment and time on second segment respectively.

1161

1162 Examples::

1163

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

1165 >>> len(a)

1166 1

1167 >>> intersection = a[0]

1168 >>> intersection.pt

1169 (374.44882952482897, 313.73458370177315)

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

1171 (0.45069111555824465, 0.5408153767394238)

1172 """

1173 s1x, s1y = s1

1174 e1x, e1y = e1

1175 s2x, s2y = s2

1176 e2x, e2y = e2

1177 if (

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

1179 ): # Parallel vertical

1180 return []

1181 if (

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

1183 ): # Parallel horizontal

1184 return []

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

1186 return []

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

1188 return []

1189 if math.isclose(e1x, s1x):

1190 x = s1x

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

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

1193 pt = (x, y)

1194 return [

1195 Intersection(

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

1197 )

1198 ]

1199 if math.isclose(s2x, e2x):

1200 x = s2x

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

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

1203 pt = (x, y)

1204 return [

1205 Intersection(

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

1207 )

1208 ]

1209

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

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

1212 if math.isclose(slope12, slope34):

1213 return []

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

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

1216 pt = (x, y)

1217 if _both_points_are_on_same_side_of_origin(

1218 pt, e1, s1

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

1220 return [

1221 Intersection(

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

1223 )

1224 ]

1225 return []

1226

1227

1228def _alignment_transformation(segment):

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

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

1231 # intersections with the segment.

1232 start = segment[0]

1233 end = segment[-1]

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

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

1236

1237

1238def _curve_line_intersections_t(curve, line):

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

1240 if len(curve) == 3:

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

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

1243 elif len(curve) == 4:

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

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

1246 else:

1247 raise ValueError("Unknown curve degree")

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

1249

1250

1251def curveLineIntersections(curve, line):

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

1253

1254 Args:

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

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

1257

1258 Returns:

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

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

1261 segment and time on second segment respectively.

1262

1263 Examples::

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

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

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

1267 >>> len(intersections)

1268 3

1269 >>> intersections[0].pt

1270 (84.9000930760723, 189.87306176459828)

1271 """

1272 if len(curve) == 3:

1273 pointFinder = quadraticPointAtT

1274 elif len(curve) == 4:

1275 pointFinder = cubicPointAtT

1276 else:

1277 raise ValueError("Unknown curve degree")

1278 intersections = []

1279 for t in _curve_line_intersections_t(curve, line):

1280 pt = pointFinder(*curve, t)

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

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

1283 line_t = _line_t_of_pt(*line, pt)

1284 pt = linePointAtT(*line, line_t)

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

1286 return intersections

1287

1288

1289def _curve_bounds(c):

1290 if len(c) == 3:

1291 return calcQuadraticBounds(*c)

1292 elif len(c) == 4:

1293 return calcCubicBounds(*c)

1294 raise ValueError("Unknown curve degree")

1295

1296

1297def _split_segment_at_t(c, t):

1298 if len(c) == 2:

1299 s, e = c

1300 midpoint = linePointAtT(s, e, t)

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

1302 if len(c) == 3:

1303 return splitQuadraticAtT(*c, t)

1304 elif len(c) == 4:

1305 return splitCubicAtT(*c, t)

1306 raise ValueError("Unknown curve degree")

1307

1308

1309def _curve_curve_intersections_t(

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

1311):

1312 bounds1 = _curve_bounds(curve1)

1313 bounds2 = _curve_bounds(curve2)

1314

1315 if not range1:

1316 range1 = (0.0, 1.0)

1317 if not range2:

1318 range2 = (0.0, 1.0)

1319

1320 # If bounds don't intersect, go home

1321 intersects, _ = sectRect(bounds1, bounds2)

1322 if not intersects:

1323 return []

1324

1325 def midpoint(r):

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

1327

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

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

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

1331

1332 c11, c12 = _split_segment_at_t(curve1, 0.5)

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

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

1335

1336 c21, c22 = _split_segment_at_t(curve2, 0.5)

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

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

1339

1340 found = []

1341 found.extend(

1342 _curve_curve_intersections_t(

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

1344 )

1345 )

1346 found.extend(

1347 _curve_curve_intersections_t(

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

1349 )

1350 )

1351 found.extend(

1352 _curve_curve_intersections_t(

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

1354 )

1355 )

1356 found.extend(

1357 _curve_curve_intersections_t(

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

1359 )

1360 )

1361

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

1363 seen = set()

1364 unique_values = []

1365

1366 for ts in found:

1367 key = unique_key(ts)

1368 if key in seen:

1369 continue

1370 seen.add(key)

1371 unique_values.append(ts)

1372

1373 return unique_values

1374

1375

1376def _is_linelike(segment):

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

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

1379

1380

1381def curveCurveIntersections(curve1, curve2):

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

1383

1384 Args:

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

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

1387

1388 Returns:

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

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

1391 segment and time on second segment respectively.

1392

1393 Examples::

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

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

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

1397 >>> len(intersections)

1398 3

1399 >>> intersections[0].pt

1400 (81.7831487395506, 109.88904552375288)

1401 """

1402 if _is_linelike(curve1):

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

1404 if _is_linelike(curve2):

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

1406 return lineLineIntersections(*line1, *line2)

1407 else:

1408 return curveLineIntersections(curve2, line1)

1409 elif _is_linelike(curve2):

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

1411 return curveLineIntersections(curve1, line2)

1412

1413 intersection_ts = _curve_curve_intersections_t(curve1, curve2)

1414 return [

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

1416 for ts in intersection_ts

1417 ]

1418

1419

1420def segmentSegmentIntersections(seg1, seg2):

1421 """Finds intersections between two segments.

1422

1423 Args:

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

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

1426

1427 Returns:

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

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

1430 segment and time on second segment respectively.

1431

1432 Examples::

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

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

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

1436 >>> len(intersections)

1437 3

1438 >>> intersections[0].pt

1439 (81.7831487395506, 109.88904552375288)

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

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

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

1443 >>> len(intersections)

1444 3

1445 >>> intersections[0].pt

1446 (84.9000930760723, 189.87306176459828)

1447

1448 """

1449 # Arrange by degree

1450 swapped = False

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

1452 seg2, seg1 = seg1, seg2

1453 swapped = True

1454 if len(seg1) > 2:

1455 if len(seg2) > 2:

1456 intersections = curveCurveIntersections(seg1, seg2)

1457 else:

1458 intersections = curveLineIntersections(seg1, seg2)

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

1460 intersections = lineLineIntersections(*seg1, *seg2)

1461 else:

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

1463 if not swapped:

1464 return intersections

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

1466

1467

1468def _segmentrepr(obj):

1469 """

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

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

1472 """

1473 try:

1474 it = iter(obj)

1475 except TypeError:

1476 return "%g" % obj

1477 else:

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

1479

1480

1481def printSegments(segments):

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

1483 segments on a single line as a tuple.

1484 """

1485 for segment in segments:

1486 print(_segmentrepr(segment))

1487

1488

1489if __name__ == "__main__":

1490 import sys

1491 import doctest

1492

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