cu2qu.py 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569
  1. # cython: language_level=3
  2. # distutils: define_macros=CYTHON_TRACE_NOGIL=1
  3. # Copyright 2015 Google Inc. All Rights Reserved.
  4. #
  5. # Licensed under the Apache License, Version 2.0 (the "License");
  6. # you may not use this file except in compliance with the License.
  7. # You may obtain a copy of the License at
  8. #
  9. # http://www.apache.org/licenses/LICENSE-2.0
  10. #
  11. # Unless required by applicable law or agreed to in writing, software
  12. # distributed under the License is distributed on an "AS IS" BASIS,
  13. # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14. # See the License for the specific language governing permissions and
  15. # limitations under the License.
  16. try:
  17. import cython
  18. except (AttributeError, ImportError):
  19. # if cython not installed, use mock module with no-op decorators and types
  20. from fontTools.misc import cython
  21. COMPILED = cython.compiled
  22. import math
  23. from .errors import Error as Cu2QuError, ApproxNotFoundError
  24. __all__ = ["curve_to_quadratic", "curves_to_quadratic"]
  25. MAX_N = 100
  26. NAN = float("NaN")
  27. @cython.cfunc
  28. @cython.inline
  29. @cython.returns(cython.double)
  30. @cython.locals(v1=cython.complex, v2=cython.complex, result=cython.double)
  31. def dot(v1, v2):
  32. """Return the dot product of two vectors.
  33. Args:
  34. v1 (complex): First vector.
  35. v2 (complex): Second vector.
  36. Returns:
  37. double: Dot product.
  38. """
  39. result = (v1 * v2.conjugate()).real
  40. # When vectors are perpendicular (i.e. dot product is 0), the above expression may
  41. # yield slightly different results when running in pure Python vs C/Cython,
  42. # both of which are correct within IEEE-754 floating-point precision.
  43. # It's probably due to the different order of operations and roundings in each
  44. # implementation. Because we are using the result in a denominator and catching
  45. # ZeroDivisionError (see `calc_intersect`), it's best to normalize the result here.
  46. if abs(result) < 1e-15:
  47. result = 0.0
  48. return result
  49. @cython.cfunc
  50. @cython.locals(z=cython.complex, den=cython.double)
  51. @cython.locals(zr=cython.double, zi=cython.double)
  52. def _complex_div_by_real(z, den):
  53. """Divide complex by real using Python's method (two separate divisions).
  54. This ensures bit-exact compatibility with Python's complex division,
  55. avoiding C's multiply-by-reciprocal optimization that can cause 1 ULP differences
  56. on some platforms/compilers (e.g. clang on macOS arm64).
  57. https://github.com/fonttools/fonttools/issues/3928
  58. """
  59. zr = z.real
  60. zi = z.imag
  61. return complex(zr / den, zi / den)
  62. @cython.cfunc
  63. @cython.inline
  64. @cython.locals(a=cython.complex, b=cython.complex, c=cython.complex, d=cython.complex)
  65. @cython.locals(
  66. _1=cython.complex, _2=cython.complex, _3=cython.complex, _4=cython.complex
  67. )
  68. def calc_cubic_points(a, b, c, d):
  69. _1 = d
  70. _2 = _complex_div_by_real(c, 3.0) + d
  71. _3 = _complex_div_by_real(b + c, 3.0) + _2
  72. _4 = a + d + c + b
  73. return _1, _2, _3, _4
  74. @cython.cfunc
  75. @cython.inline
  76. @cython.locals(
  77. p0=cython.complex, p1=cython.complex, p2=cython.complex, p3=cython.complex
  78. )
  79. @cython.locals(a=cython.complex, b=cython.complex, c=cython.complex, d=cython.complex)
  80. def calc_cubic_parameters(p0, p1, p2, p3):
  81. c = (p1 - p0) * 3.0
  82. b = (p2 - p1) * 3.0 - c
  83. d = p0
  84. a = p3 - d - c - b
  85. return a, b, c, d
  86. @cython.cfunc
  87. @cython.inline
  88. @cython.locals(
  89. p0=cython.complex, p1=cython.complex, p2=cython.complex, p3=cython.complex
  90. )
  91. def split_cubic_into_n_iter(p0, p1, p2, p3, n):
  92. """Split a cubic Bezier into n equal parts.
  93. Splits the curve into `n` equal parts by curve time.
  94. (t=0..1/n, t=1/n..2/n, ...)
  95. Args:
  96. p0 (complex): Start point of curve.
  97. p1 (complex): First handle of curve.
  98. p2 (complex): Second handle of curve.
  99. p3 (complex): End point of curve.
  100. Returns:
  101. An iterator yielding the control points (four complex values) of the
  102. subcurves.
  103. """
  104. # Hand-coded special-cases
  105. if n == 2:
  106. return iter(split_cubic_into_two(p0, p1, p2, p3))
  107. if n == 3:
  108. return iter(split_cubic_into_three(p0, p1, p2, p3))
  109. if n == 4:
  110. a, b = split_cubic_into_two(p0, p1, p2, p3)
  111. return iter(
  112. split_cubic_into_two(a[0], a[1], a[2], a[3])
  113. + split_cubic_into_two(b[0], b[1], b[2], b[3])
  114. )
  115. if n == 6:
  116. a, b = split_cubic_into_two(p0, p1, p2, p3)
  117. return iter(
  118. split_cubic_into_three(a[0], a[1], a[2], a[3])
  119. + split_cubic_into_three(b[0], b[1], b[2], b[3])
  120. )
  121. return _split_cubic_into_n_gen(p0, p1, p2, p3, n)
  122. @cython.locals(
  123. p0=cython.complex,
  124. p1=cython.complex,
  125. p2=cython.complex,
  126. p3=cython.complex,
  127. n=cython.int,
  128. )
  129. @cython.locals(a=cython.complex, b=cython.complex, c=cython.complex, d=cython.complex)
  130. @cython.locals(
  131. dt=cython.double, delta_2=cython.double, delta_3=cython.double, i=cython.int
  132. )
  133. @cython.locals(
  134. a1=cython.complex, b1=cython.complex, c1=cython.complex, d1=cython.complex
  135. )
  136. def _split_cubic_into_n_gen(p0, p1, p2, p3, n):
  137. a, b, c, d = calc_cubic_parameters(p0, p1, p2, p3)
  138. dt = 1 / n
  139. delta_2 = dt * dt
  140. delta_3 = dt * delta_2
  141. for i in range(n):
  142. t1 = i * dt
  143. t1_2 = t1 * t1
  144. # calc new a, b, c and d
  145. a1 = a * delta_3
  146. b1 = (3 * a * t1 + b) * delta_2
  147. c1 = (2 * b * t1 + c + 3 * a * t1_2) * dt
  148. d1 = a * t1 * t1_2 + b * t1_2 + c * t1 + d
  149. yield calc_cubic_points(a1, b1, c1, d1)
  150. @cython.cfunc
  151. @cython.inline
  152. @cython.locals(
  153. p0=cython.complex, p1=cython.complex, p2=cython.complex, p3=cython.complex
  154. )
  155. @cython.locals(mid=cython.complex, deriv3=cython.complex)
  156. def split_cubic_into_two(p0, p1, p2, p3):
  157. """Split a cubic Bezier into two equal parts.
  158. Splits the curve into two equal parts at t = 0.5
  159. Args:
  160. p0 (complex): Start point of curve.
  161. p1 (complex): First handle of curve.
  162. p2 (complex): Second handle of curve.
  163. p3 (complex): End point of curve.
  164. Returns:
  165. tuple: Two cubic Beziers (each expressed as a tuple of four complex
  166. values).
  167. """
  168. mid = (p0 + 3 * (p1 + p2) + p3) * 0.125
  169. deriv3 = (p3 + p2 - p1 - p0) * 0.125
  170. return (
  171. (p0, (p0 + p1) * 0.5, mid - deriv3, mid),
  172. (mid, mid + deriv3, (p2 + p3) * 0.5, p3),
  173. )
  174. @cython.cfunc
  175. @cython.inline
  176. @cython.locals(
  177. p0=cython.complex,
  178. p1=cython.complex,
  179. p2=cython.complex,
  180. p3=cython.complex,
  181. )
  182. @cython.locals(
  183. mid1=cython.complex,
  184. deriv1=cython.complex,
  185. mid2=cython.complex,
  186. deriv2=cython.complex,
  187. )
  188. def split_cubic_into_three(p0, p1, p2, p3):
  189. """Split a cubic Bezier into three equal parts.
  190. Splits the curve into three equal parts at t = 1/3 and t = 2/3
  191. Args:
  192. p0 (complex): Start point of curve.
  193. p1 (complex): First handle of curve.
  194. p2 (complex): Second handle of curve.
  195. p3 (complex): End point of curve.
  196. Returns:
  197. tuple: Three cubic Beziers (each expressed as a tuple of four complex
  198. values).
  199. """
  200. mid1 = (8 * p0 + 12 * p1 + 6 * p2 + p3) * (1 / 27)
  201. deriv1 = (p3 + 3 * p2 - 4 * p0) * (1 / 27)
  202. mid2 = (p0 + 6 * p1 + 12 * p2 + 8 * p3) * (1 / 27)
  203. deriv2 = (4 * p3 - 3 * p1 - p0) * (1 / 27)
  204. return (
  205. (p0, (2 * p0 + p1) / 3.0, mid1 - deriv1, mid1),
  206. (mid1, mid1 + deriv1, mid2 - deriv2, mid2),
  207. (mid2, mid2 + deriv2, (p2 + 2 * p3) / 3.0, p3),
  208. )
  209. @cython.cfunc
  210. @cython.inline
  211. @cython.returns(cython.complex)
  212. @cython.locals(
  213. t=cython.double,
  214. p0=cython.complex,
  215. p1=cython.complex,
  216. p2=cython.complex,
  217. p3=cython.complex,
  218. )
  219. @cython.locals(_p1=cython.complex, _p2=cython.complex)
  220. def cubic_approx_control(t, p0, p1, p2, p3):
  221. """Approximate a cubic Bezier using a quadratic one.
  222. Args:
  223. t (double): Position of control point.
  224. p0 (complex): Start point of curve.
  225. p1 (complex): First handle of curve.
  226. p2 (complex): Second handle of curve.
  227. p3 (complex): End point of curve.
  228. Returns:
  229. complex: Location of candidate control point on quadratic curve.
  230. """
  231. _p1 = p0 + (p1 - p0) * 1.5
  232. _p2 = p3 + (p2 - p3) * 1.5
  233. return _p1 + (_p2 - _p1) * t
  234. @cython.cfunc
  235. @cython.inline
  236. @cython.returns(cython.complex)
  237. @cython.locals(a=cython.complex, b=cython.complex, c=cython.complex, d=cython.complex)
  238. @cython.locals(ab=cython.complex, cd=cython.complex, p=cython.complex, h=cython.double)
  239. def calc_intersect(a, b, c, d):
  240. """Calculate the intersection of two lines.
  241. Args:
  242. a (complex): Start point of first line.
  243. b (complex): End point of first line.
  244. c (complex): Start point of second line.
  245. d (complex): End point of second line.
  246. Returns:
  247. complex: Location of intersection if one present, ``complex(NaN,NaN)``
  248. if no intersection was found.
  249. """
  250. ab = b - a
  251. cd = d - c
  252. p = ab * 1j
  253. try:
  254. h = dot(p, a - c) / dot(p, cd)
  255. except ZeroDivisionError:
  256. # if 3 or 4 points are equal, we do have an intersection despite the zero-div:
  257. # return one of the off-curves so that the algorithm can attempt a one-curve
  258. # solution if it's within tolerance:
  259. # https://github.com/linebender/kurbo/pull/484
  260. if b == c and (a == b or c == d):
  261. return b
  262. return complex(NAN, NAN)
  263. return c + cd * h
  264. @cython.cfunc
  265. @cython.returns(cython.int)
  266. @cython.locals(
  267. tolerance=cython.double,
  268. p0=cython.complex,
  269. p1=cython.complex,
  270. p2=cython.complex,
  271. p3=cython.complex,
  272. )
  273. @cython.locals(mid=cython.complex, deriv3=cython.complex)
  274. def cubic_farthest_fit_inside(p0, p1, p2, p3, tolerance):
  275. """Check if a cubic Bezier lies within a given distance of the origin.
  276. "Origin" means *the* origin (0,0), not the start of the curve. Note that no
  277. checks are made on the start and end positions of the curve; this function
  278. only checks the inside of the curve.
  279. Args:
  280. p0 (complex): Start point of curve.
  281. p1 (complex): First handle of curve.
  282. p2 (complex): Second handle of curve.
  283. p3 (complex): End point of curve.
  284. tolerance (double): Distance from origin.
  285. Returns:
  286. bool: True if the cubic Bezier ``p`` entirely lies within a distance
  287. ``tolerance`` of the origin, False otherwise.
  288. """
  289. # First check p2 then p1, as p2 has higher error early on.
  290. if abs(p2) <= tolerance and abs(p1) <= tolerance:
  291. return True
  292. # Split.
  293. mid = (p0 + 3 * (p1 + p2) + p3) * 0.125
  294. if abs(mid) > tolerance:
  295. return False
  296. deriv3 = (p3 + p2 - p1 - p0) * 0.125
  297. return cubic_farthest_fit_inside(
  298. p0, (p0 + p1) * 0.5, mid - deriv3, mid, tolerance
  299. ) and cubic_farthest_fit_inside(mid, mid + deriv3, (p2 + p3) * 0.5, p3, tolerance)
  300. @cython.cfunc
  301. @cython.inline
  302. @cython.locals(tolerance=cython.double)
  303. @cython.locals(
  304. q1=cython.complex,
  305. c0=cython.complex,
  306. c1=cython.complex,
  307. c2=cython.complex,
  308. c3=cython.complex,
  309. )
  310. def cubic_approx_quadratic(cubic, tolerance):
  311. """Approximate a cubic Bezier with a single quadratic within a given tolerance.
  312. Args:
  313. cubic (sequence): Four complex numbers representing control points of
  314. the cubic Bezier curve.
  315. tolerance (double): Permitted deviation from the original curve.
  316. Returns:
  317. Three complex numbers representing control points of the quadratic
  318. curve if it fits within the given tolerance, or ``None`` if no suitable
  319. curve could be calculated.
  320. """
  321. q1 = calc_intersect(cubic[0], cubic[1], cubic[2], cubic[3])
  322. if math.isnan(q1.imag):
  323. return None
  324. c0 = cubic[0]
  325. c3 = cubic[3]
  326. c1 = c0 + (q1 - c0) * (2 / 3)
  327. c2 = c3 + (q1 - c3) * (2 / 3)
  328. if not cubic_farthest_fit_inside(0, c1 - cubic[1], c2 - cubic[2], 0, tolerance):
  329. return None
  330. return c0, q1, c3
  331. @cython.cfunc
  332. @cython.locals(n=cython.int, tolerance=cython.double)
  333. @cython.locals(i=cython.int)
  334. @cython.locals(all_quadratic=cython.int)
  335. @cython.locals(
  336. c0=cython.complex, c1=cython.complex, c2=cython.complex, c3=cython.complex
  337. )
  338. @cython.locals(
  339. q0=cython.complex,
  340. q1=cython.complex,
  341. next_q1=cython.complex,
  342. q2=cython.complex,
  343. d1=cython.complex,
  344. )
  345. def cubic_approx_spline(cubic, n, tolerance, all_quadratic):
  346. """Approximate a cubic Bezier curve with a spline of n quadratics.
  347. Args:
  348. cubic (sequence): Four complex numbers representing control points of
  349. the cubic Bezier curve.
  350. n (int): Number of quadratic Bezier curves in the spline.
  351. tolerance (double): Permitted deviation from the original curve.
  352. Returns:
  353. A list of ``n+2`` complex numbers, representing control points of the
  354. quadratic spline if it fits within the given tolerance, or ``None`` if
  355. no suitable spline could be calculated.
  356. """
  357. if n == 1:
  358. return cubic_approx_quadratic(cubic, tolerance)
  359. if n == 2 and all_quadratic == False:
  360. return cubic
  361. cubics = split_cubic_into_n_iter(cubic[0], cubic[1], cubic[2], cubic[3], n)
  362. # calculate the spline of quadratics and check errors at the same time.
  363. next_cubic = next(cubics)
  364. next_q1 = cubic_approx_control(
  365. 0, next_cubic[0], next_cubic[1], next_cubic[2], next_cubic[3]
  366. )
  367. q2 = cubic[0]
  368. d1 = 0j
  369. spline = [cubic[0], next_q1]
  370. for i in range(1, n + 1):
  371. # Current cubic to convert
  372. c0, c1, c2, c3 = next_cubic
  373. # Current quadratic approximation of current cubic
  374. q0 = q2
  375. q1 = next_q1
  376. if i < n:
  377. next_cubic = next(cubics)
  378. next_q1 = cubic_approx_control(
  379. i / (n - 1), next_cubic[0], next_cubic[1], next_cubic[2], next_cubic[3]
  380. )
  381. spline.append(next_q1)
  382. q2 = (q1 + next_q1) * 0.5
  383. else:
  384. q2 = c3
  385. # End-point deltas
  386. d0 = d1
  387. d1 = q2 - c3
  388. if abs(d1) > tolerance or not cubic_farthest_fit_inside(
  389. d0,
  390. q0 + (q1 - q0) * (2 / 3) - c1,
  391. q2 + (q1 - q2) * (2 / 3) - c2,
  392. d1,
  393. tolerance,
  394. ):
  395. return None
  396. spline.append(cubic[3])
  397. return spline
  398. @cython.locals(max_err=cython.double)
  399. @cython.locals(n=cython.int)
  400. @cython.locals(all_quadratic=cython.int)
  401. def curve_to_quadratic(curve, max_err, all_quadratic=True):
  402. """Approximate a cubic Bezier curve with a spline of n quadratics.
  403. Args:
  404. cubic (sequence): Four 2D tuples representing control points of
  405. the cubic Bezier curve.
  406. max_err (double): Permitted deviation from the original curve.
  407. all_quadratic (bool): If True (default) returned value is a
  408. quadratic spline. If False, it's either a single quadratic
  409. curve or a single cubic curve.
  410. Returns:
  411. If all_quadratic is True: A list of 2D tuples representing
  412. control points of the quadratic spline.
  413. If all_quadratic is False: Either a quadratic curve (if length
  414. of output is 3), or a cubic curve (if length of output is 4).
  415. Raises:
  416. fontTools.cu2qu.errors.ApproxNotFoundError: if no suitable
  417. approximation can be found with the given parameters.
  418. """
  419. curve = [complex(*p) for p in curve]
  420. for n in range(1, MAX_N + 1):
  421. spline = cubic_approx_spline(curve, n, max_err, all_quadratic)
  422. if spline is not None:
  423. # done. go home
  424. return [(s.real, s.imag) for s in spline]
  425. raise ApproxNotFoundError(curve)
  426. @cython.locals(l=cython.int, last_i=cython.int, i=cython.int)
  427. @cython.locals(all_quadratic=cython.int)
  428. def curves_to_quadratic(curves, max_errors, all_quadratic=True):
  429. """Return quadratic Bezier splines approximating the input cubic Beziers.
  430. Args:
  431. curves: A sequence of *n* curves, each curve being a sequence of four
  432. 2D tuples.
  433. max_errors: A sequence of *n* floats representing the maximum permissible
  434. deviation from each of the cubic Bezier curves.
  435. all_quadratic (bool): If True (default) returned values are a
  436. quadratic spline. If False, they are either a single quadratic
  437. curve or a single cubic curve.
  438. Example::
  439. >>> curves_to_quadratic( [
  440. ... [ (50,50), (100,100), (150,100), (200,50) ],
  441. ... [ (75,50), (120,100), (150,75), (200,60) ]
  442. ... ], [1,1] )
  443. [[(50.0, 50.0), (75.0, 75.0), (125.0, 91.66666666666666), (175.0, 75.0), (200.0, 50.0)], [(75.0, 50.0), (97.5, 75.0), (135.41666666666666, 82.08333333333333), (175.0, 67.5), (200.0, 60.0)]]
  444. The returned splines have "implied oncurve points" suitable for use in
  445. TrueType ``glif`` outlines - i.e. in the first spline returned above,
  446. the first quadratic segment runs from (50,50) to
  447. ( (75 + 125)/2 , (120 + 91.666..)/2 ) = (100, 83.333...).
  448. Returns:
  449. If all_quadratic is True, a list of splines, each spline being a list
  450. of 2D tuples. If ``curves`` is empty, returns an empty list.
  451. If all_quadratic is False, a list of curves, each curve being a quadratic
  452. (length 3), or cubic (length 4).
  453. Raises:
  454. ValueError: if ``max_errors`` does not match the number of curves.
  455. fontTools.cu2qu.errors.ApproxNotFoundError: if no suitable approximation
  456. can be found for all curves with the given parameters.
  457. """
  458. curves = [[complex(*p) for p in curve] for curve in curves]
  459. if len(max_errors) != len(curves):
  460. raise ValueError("max_errors must match the number of curves")
  461. if not curves:
  462. return []
  463. l = len(curves)
  464. splines = [None] * l
  465. last_i = i = 0
  466. n = 1
  467. while True:
  468. spline = cubic_approx_spline(curves[i], n, max_errors[i], all_quadratic)
  469. if spline is None:
  470. if n == MAX_N:
  471. break
  472. n += 1
  473. last_i = i
  474. continue
  475. splines[i] = spline
  476. i = (i + 1) % l
  477. if i == last_i:
  478. # done. go home
  479. return [[(s.real, s.imag) for s in spline] for spline in splines]
  480. raise ApproxNotFoundError(curves)