2015-11-12 16:21:35 -08:00
|
|
|
# Copyright 2015 Google Inc. All Rights Reserved.
|
|
|
|
#
|
|
|
|
# Licensed under the Apache License, Version 2.0 (the "License");
|
|
|
|
# you may not use this file except in compliance with the License.
|
|
|
|
# You may obtain a copy of the License at
|
|
|
|
#
|
|
|
|
# http://www.apache.org/licenses/LICENSE-2.0
|
|
|
|
#
|
|
|
|
# Unless required by applicable law or agreed to in writing, software
|
|
|
|
# distributed under the License is distributed on an "AS IS" BASIS,
|
|
|
|
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
|
|
|
|
# See the License for the specific language governing permissions and
|
|
|
|
# limitations under the License.
|
|
|
|
|
|
|
|
|
2015-12-04 13:07:32 -08:00
|
|
|
from __future__ import print_function, division, absolute_import
|
|
|
|
|
2015-11-12 16:21:35 -08:00
|
|
|
from math import hypot
|
|
|
|
|
2015-12-08 12:41:33 -08:00
|
|
|
__all__ = ['curve_to_quadratic', 'curves_to_quadratic']
|
|
|
|
|
2015-12-08 12:44:29 -08:00
|
|
|
MAX_N = 100
|
|
|
|
|
2015-11-12 16:21:35 -08:00
|
|
|
|
2016-07-22 23:40:44 -07:00
|
|
|
def calcCubicPoints(a, b, c, d):
|
|
|
|
_1 = d
|
|
|
|
_2 = (c / 3.0) + d
|
|
|
|
_3 = (b + c) / 3.0 + _2
|
|
|
|
_4 = a + d + c + b
|
|
|
|
return _1, _2, _3, _4
|
|
|
|
|
|
|
|
def _splitCubicAtT(a, b, c, d, *ts):
|
|
|
|
ts = list(ts)
|
|
|
|
ts.insert(0, 0.0)
|
|
|
|
ts.append(1.0)
|
|
|
|
segments = []
|
|
|
|
for i in range(len(ts) - 1):
|
|
|
|
t1 = ts[i]
|
|
|
|
t2 = ts[i+1]
|
|
|
|
delta = (t2 - t1)
|
|
|
|
|
|
|
|
delta_2 = delta*delta
|
|
|
|
delta_3 = delta*delta_2
|
|
|
|
t1_2 = t1*t1
|
|
|
|
t1_3 = t1*t1_2
|
|
|
|
|
|
|
|
# calc new a, b, c and d
|
|
|
|
a1 = a * delta_3
|
|
|
|
b1 = (3*a*t1 + b) * delta_2
|
|
|
|
c1 = (2*b*t1 + c + 3*a*t1_2) * delta
|
|
|
|
d1 = a*t1_3 + b*t1_2 + c*t1 + d
|
|
|
|
pt1, pt2, pt3, pt4 = calcCubicPoints(a1, b1, c1, d1)
|
|
|
|
segments.append((pt1, pt2, pt3, pt4))
|
|
|
|
return segments
|
|
|
|
|
|
|
|
def calcCubicParameters(pt1, pt2, pt3, pt4):
|
|
|
|
c = (pt2 -pt1) * 3.0
|
|
|
|
b = (pt3 - pt2) * 3.0 - c
|
|
|
|
d = pt1
|
|
|
|
a = pt4 - d - c - b
|
|
|
|
return a, b, c, d
|
|
|
|
|
|
|
|
def splitCubicAtT(pt1, pt2, pt3, pt4, *ts):
|
|
|
|
a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)
|
|
|
|
return _splitCubicAtT(a, b, c, d, *ts)
|
|
|
|
|
2016-07-24 13:40:44 -04:00
|
|
|
def splitCubicIntoTwo(pt1, pt2, pt3, pt4):
|
|
|
|
mid = (pt1+3*(pt2+pt3)+pt4)*.125
|
|
|
|
deriv3 = (pt4+pt3-pt2-pt1)*.125
|
|
|
|
return ((pt1, (pt1+pt2)*.5, mid-deriv3, mid),
|
|
|
|
(mid, mid+deriv3, (pt3+pt4)*.5, pt4))
|
|
|
|
|
2016-07-24 14:10:41 -04:00
|
|
|
def splitCubicIntoThree(pt1, pt2, pt3, pt4, _27=1/27.):
|
|
|
|
mid1 = (pt1*8+pt2*12+pt3*6+pt4)*_27
|
|
|
|
deriv1 = (pt4+pt3*3-pt1*4)*_27
|
|
|
|
mid2 = (pt1+pt2*6+pt3*12+pt4*8)*_27
|
|
|
|
deriv2 = (pt4*4-pt2*3-pt1)*_27
|
|
|
|
return ((pt1, (pt1*2+pt2)/3., mid1-deriv1, mid1),
|
|
|
|
(mid1, mid1+deriv1, mid2-deriv2, mid2),
|
|
|
|
(mid2, mid2+deriv2, (pt3+pt4*2)/3., pt4))
|
|
|
|
|
2016-07-22 23:40:44 -07:00
|
|
|
|
2015-12-07 14:24:23 +00:00
|
|
|
class Cu2QuError(Exception):
|
|
|
|
pass
|
|
|
|
|
|
|
|
|
|
|
|
class ApproxNotFoundError(Cu2QuError):
|
|
|
|
def __init__(self, curve, error=None):
|
|
|
|
if error is None:
|
|
|
|
message = "no approximation found: %s" % curve
|
|
|
|
else:
|
|
|
|
message = ("approximation error exceeds max tolerance: %s, "
|
|
|
|
"error=%g" % (curve, error))
|
|
|
|
super(Cu2QuError, self).__init__(message)
|
|
|
|
self.curve = curve
|
|
|
|
self.error = error
|
|
|
|
|
|
|
|
|
2015-12-01 13:20:09 -08:00
|
|
|
def dot(v1, v2):
|
|
|
|
"""Return the dot product of two vectors."""
|
2016-07-22 23:40:44 -07:00
|
|
|
return v1.real * v2.real + v1.imag * v2.imag
|
2016-07-22 02:21:41 -07:00
|
|
|
|
2015-11-12 16:21:35 -08:00
|
|
|
|
2016-07-22 02:21:41 -07:00
|
|
|
def cubic_from_quadratic(p):
|
2016-07-22 23:40:44 -07:00
|
|
|
return (p[0], p[0]+(p[1]-p[0])*(2./3), p[2]+(p[1]-p[2])*(2./3), p[2])
|
2016-07-22 02:21:41 -07:00
|
|
|
|
|
|
|
|
2016-07-22 01:13:12 -07:00
|
|
|
def cubic_approx_control(p, t):
|
|
|
|
"""Approximate a cubic bezier curve with a quadratic one.
|
|
|
|
Returns the candidate control point."""
|
2015-11-12 16:21:35 -08:00
|
|
|
|
2016-07-22 23:40:44 -07:00
|
|
|
p1 = p[0]+(p[1]-p[0])*1.5
|
|
|
|
p2 = p[3]+(p[2]-p[3])*1.5
|
|
|
|
return p1+(p2-p1)*t
|
2015-11-12 16:21:35 -08:00
|
|
|
|
|
|
|
|
|
|
|
def calc_intersect(p):
|
|
|
|
"""Calculate the intersection of ab and cd, given [a, b, c, d]."""
|
|
|
|
|
|
|
|
a, b, c, d = p
|
2016-07-22 23:40:44 -07:00
|
|
|
ab = b - a
|
|
|
|
cd = d - c
|
|
|
|
p = ab * 1j
|
2015-11-12 16:21:35 -08:00
|
|
|
try:
|
2016-07-22 23:40:44 -07:00
|
|
|
h = dot(p, a - c) / dot(p, cd)
|
2015-11-12 16:21:35 -08:00
|
|
|
except ZeroDivisionError:
|
|
|
|
raise ValueError('Parallel vectors given to calc_intersect.')
|
2016-07-22 23:40:44 -07:00
|
|
|
return c + cd * h
|
2015-11-12 16:21:35 -08:00
|
|
|
|
|
|
|
|
2016-07-23 16:59:54 -07:00
|
|
|
def cubic_farthest_fit(p,tolerance):
|
|
|
|
"""Returns True if the cubic Bezier p entirely lies within a distance
|
|
|
|
tolerance of origin, False otherwise."""
|
2016-07-22 02:21:41 -07:00
|
|
|
|
2016-07-23 16:59:54 -07:00
|
|
|
if abs(p[0]) > tolerance or abs(p[3]) > tolerance:
|
|
|
|
return False
|
2016-07-22 02:21:41 -07:00
|
|
|
|
2016-07-23 16:59:54 -07:00
|
|
|
if abs(p[1]) <= tolerance and abs(p[2]) <= tolerance:
|
|
|
|
return True
|
2016-07-22 02:21:41 -07:00
|
|
|
|
|
|
|
# Split.
|
2016-07-24 13:40:44 -04:00
|
|
|
a,b = splitCubicIntoTwo(p[0], p[1], p[2], p[3])
|
2016-07-23 16:59:54 -07:00
|
|
|
return cubic_farthest_fit(a,tolerance) and cubic_farthest_fit(b,tolerance)
|
2016-07-22 02:21:41 -07:00
|
|
|
|
|
|
|
|
2016-07-23 16:59:54 -07:00
|
|
|
def cubic_cubic_fit(a,b,tolerance):
|
|
|
|
return cubic_farthest_fit((b[0] - a[0],
|
|
|
|
b[1] - a[1],
|
|
|
|
b[2] - a[2],
|
|
|
|
b[3] - a[3]), tolerance)
|
2016-07-22 02:21:41 -07:00
|
|
|
|
|
|
|
|
2016-07-23 16:59:54 -07:00
|
|
|
def cubic_quadratic_fit(a,b,tolerance):
|
|
|
|
return cubic_cubic_fit(a, cubic_from_quadratic(b), tolerance)
|
2016-07-22 02:21:41 -07:00
|
|
|
|
|
|
|
|
|
|
|
def cubic_approx_spline(p, n, tolerance):
|
2015-11-12 16:21:35 -08:00
|
|
|
"""Approximate a cubic bezier curve with a spline of n quadratics.
|
|
|
|
|
|
|
|
Returns None if n is 1 and the cubic's control vectors are parallel, since
|
|
|
|
no quadratic exists with this cubic's tangents.
|
|
|
|
"""
|
|
|
|
|
|
|
|
if n == 1:
|
|
|
|
try:
|
|
|
|
p1 = calc_intersect(p)
|
|
|
|
except ValueError:
|
|
|
|
return None
|
2016-07-22 02:21:41 -07:00
|
|
|
quad = (p[0], p1, p[3])
|
2016-07-23 16:59:54 -07:00
|
|
|
if not cubic_quadratic_fit(p, quad, tolerance):
|
2016-07-22 02:21:41 -07:00
|
|
|
return None
|
|
|
|
return quad
|
2015-11-12 16:21:35 -08:00
|
|
|
|
|
|
|
spline = [p[0]]
|
2016-07-24 13:42:21 -04:00
|
|
|
if n == 2:
|
|
|
|
segments = splitCubicIntoTwo(p[0], p[1], p[2], p[3])
|
2016-07-24 14:10:41 -04:00
|
|
|
elif n == 3:
|
|
|
|
segments = splitCubicIntoThree(p[0], p[1], p[2], p[3])
|
2016-07-24 13:42:21 -04:00
|
|
|
else:
|
|
|
|
ts = [i / n for i in range(1, n)]
|
|
|
|
segments = splitCubicAtT(p[0], p[1], p[2], p[3], *ts)
|
2015-11-12 16:21:35 -08:00
|
|
|
for i in range(len(segments)):
|
2016-07-22 01:13:12 -07:00
|
|
|
spline.append(cubic_approx_control(segments[i], i / (n - 1)))
|
2015-11-12 16:21:35 -08:00
|
|
|
spline.append(p[3])
|
|
|
|
|
2016-07-22 02:21:41 -07:00
|
|
|
for i in range(1,n+1):
|
|
|
|
if i == 1:
|
2016-07-23 17:11:45 -07:00
|
|
|
segment = (spline[0],spline[1],(spline[1]+spline[2])*.5)
|
|
|
|
elif i == n:
|
2016-07-22 23:40:44 -07:00
|
|
|
segment = (spline[-3]+spline[-2])*.5,spline[-2],spline[-1]
|
2016-07-23 17:11:45 -07:00
|
|
|
else:
|
2016-07-22 23:40:44 -07:00
|
|
|
segment = (spline[i-1]+spline[i])*.5, spline[i], (spline[i]+spline[i+1])*.5
|
2015-11-12 16:21:35 -08:00
|
|
|
|
2016-07-23 16:59:54 -07:00
|
|
|
if not cubic_quadratic_fit(segments[i-1], segment, tolerance):
|
2016-07-23 17:11:45 -07:00
|
|
|
return None
|
2015-11-12 16:21:35 -08:00
|
|
|
|
2016-07-22 02:21:41 -07:00
|
|
|
return spline
|
2015-11-12 16:21:35 -08:00
|
|
|
|
|
|
|
|
2015-12-08 12:44:29 -08:00
|
|
|
def curve_to_quadratic(p, max_err):
|
2015-12-07 14:24:23 +00:00
|
|
|
"""Return a quadratic spline approximating this cubic bezier, and
|
|
|
|
the error of approximation.
|
|
|
|
Raise 'ApproxNotFoundError' if no suitable approximation can be found
|
|
|
|
with the given parameters.
|
|
|
|
"""
|
2015-11-12 16:21:35 -08:00
|
|
|
|
2016-07-22 23:40:44 -07:00
|
|
|
p = [complex(*P) for P in p]
|
2015-12-07 14:24:23 +00:00
|
|
|
spline, error = None, None
|
2015-12-08 12:44:29 -08:00
|
|
|
for n in range(1, MAX_N + 1):
|
2016-07-22 02:21:41 -07:00
|
|
|
spline = cubic_approx_spline(p, n, max_err)
|
|
|
|
if spline is not None:
|
2015-11-12 16:21:35 -08:00
|
|
|
break
|
2015-12-07 14:24:23 +00:00
|
|
|
else:
|
|
|
|
# no break: approximation not found or error exceeds tolerance
|
|
|
|
raise ApproxNotFoundError(p, error)
|
|
|
|
return spline, error
|
2015-11-12 16:21:35 -08:00
|
|
|
|
|
|
|
|
2015-12-08 12:44:29 -08:00
|
|
|
def curves_to_quadratic(curves, max_errors):
|
2015-12-07 14:24:23 +00:00
|
|
|
"""Return quadratic splines approximating these cubic beziers, and
|
|
|
|
the respective errors of approximation.
|
|
|
|
Raise 'ApproxNotFoundError' if no suitable approximation can be found
|
|
|
|
for all curves with the given parameters.
|
|
|
|
"""
|
2015-11-12 16:21:35 -08:00
|
|
|
|
2016-07-22 23:40:44 -07:00
|
|
|
curves = [[complex(*P) for P in p] for p in curves]
|
2015-12-10 12:46:16 -08:00
|
|
|
num_curves = len(curves)
|
|
|
|
assert len(max_errors) == num_curves
|
|
|
|
|
|
|
|
splines = [None] * num_curves
|
|
|
|
errors = [None] * num_curves
|
2015-12-08 12:44:29 -08:00
|
|
|
for n in range(1, MAX_N + 1):
|
2016-07-22 02:21:41 -07:00
|
|
|
splines = [cubic_approx_spline(c, n, e) for c,e in zip(curves,max_errors)]
|
|
|
|
if all(splines):
|
2015-11-12 16:21:35 -08:00
|
|
|
break
|
2015-12-07 14:24:23 +00:00
|
|
|
else:
|
|
|
|
# no break: raise if any spline is None or error exceeds tolerance
|
|
|
|
for c, s, error, max_err in zip(curves, splines, errors, max_errors):
|
|
|
|
if s is None or error > max_err:
|
|
|
|
raise ApproxNotFoundError(c, error)
|
|
|
|
return splines, errors
|