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