Black whole module
This commit is contained in:
Simon Cozens
parent
1a3478da0e
commit
60c1ee0107
@ -54,23 +54,31 @@ def calcCubicArcLength(pt1, pt2, pt3, pt4, tolerance=0.005):
|
||||
Returns:
|
||||
Arc length value.
|
||||
"""
|
||||
return calcCubicArcLengthC(complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4), tolerance)
|
||||
return calcCubicArcLengthC(
|
||||
complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4), tolerance
|
||||
)
|
||||
|
||||
|
||||
def _split_cubic_into_two(p0, p1, p2, p3):
|
||||
mid = (p0 + 3 * (p1 + p2) + p3) * .125
|
||||
deriv3 = (p3 + p2 - p1 - p0) * .125
|
||||
return ((p0, (p0 + p1) * .5, mid - deriv3, mid),
|
||||
(mid, mid + deriv3, (p2 + p3) * .5, p3))
|
||||
mid = (p0 + 3 * (p1 + p2) + p3) * 0.125
|
||||
deriv3 = (p3 + p2 - p1 - p0) * 0.125
|
||||
return (
|
||||
(p0, (p0 + p1) * 0.5, mid - deriv3, mid),
|
||||
(mid, mid + deriv3, (p2 + p3) * 0.5, p3),
|
||||
)
|
||||
|
||||
|
||||
def _calcCubicArcLengthCRecurse(mult, p0, p1, p2, p3):
|
||||
arch = abs(p0-p3)
|
||||
box = abs(p0-p1) + abs(p1-p2) + abs(p2-p3)
|
||||
if arch * mult >= box:
|
||||
return (arch + box) * .5
|
||||
else:
|
||||
one,two = _split_cubic_into_two(p0,p1,p2,p3)
|
||||
return _calcCubicArcLengthCRecurse(mult, *one) + _calcCubicArcLengthCRecurse(mult, *two)
|
||||
arch = abs(p0 - p3)
|
||||
box = abs(p0 - p1) + abs(p1 - p2) + abs(p2 - p3)
|
||||
if arch * mult >= box:
|
||||
return (arch + box) * 0.5
|
||||
else:
|
||||
one, two = _split_cubic_into_two(p0, p1, p2, p3)
|
||||
return _calcCubicArcLengthCRecurse(mult, *one) + _calcCubicArcLengthCRecurse(
|
||||
mult, *two
|
||||
)
|
||||
|
||||
|
||||
def calcCubicArcLengthC(pt1, pt2, pt3, pt4, tolerance=0.005):
|
||||
"""Calculates the arc length for a cubic Bezier segment.
|
||||
@ -82,7 +90,7 @@ def calcCubicArcLengthC(pt1, pt2, pt3, pt4, tolerance=0.005):
|
||||
Returns:
|
||||
Arc length value.
|
||||
"""
|
||||
mult = 1. + 1.5 * tolerance # The 1.5 is a empirical hack; no math
|
||||
mult = 1.0 + 1.5 * tolerance # The 1.5 is a empirical hack; no math
|
||||
return _calcCubicArcLengthCRecurse(mult, pt1, pt2, pt3, pt4)
|
||||
|
||||
|
||||
@ -97,7 +105,7 @@ def _dot(v1, v2):
|
||||
def _intSecAtan(x):
|
||||
# In : sympy.integrate(sp.sec(sp.atan(x)))
|
||||
# Out: x*sqrt(x**2 + 1)/2 + asinh(x)/2
|
||||
return x * math.sqrt(x**2 + 1)/2 + math.asinh(x)/2
|
||||
return x * math.sqrt(x ** 2 + 1) / 2 + math.asinh(x) / 2
|
||||
|
||||
|
||||
def calcQuadraticArcLength(pt1, pt2, pt3):
|
||||
@ -153,16 +161,16 @@ def calcQuadraticArcLengthC(pt1, pt2, pt3):
|
||||
d = d1 - d0
|
||||
n = d * 1j
|
||||
scale = abs(n)
|
||||
if scale == 0.:
|
||||
return abs(pt3-pt1)
|
||||
origDist = _dot(n,d0)
|
||||
if scale == 0.0:
|
||||
return abs(pt3 - pt1)
|
||||
origDist = _dot(n, d0)
|
||||
if abs(origDist) < epsilon:
|
||||
if _dot(d0,d1) >= 0:
|
||||
return abs(pt3-pt1)
|
||||
if _dot(d0, d1) >= 0:
|
||||
return abs(pt3 - pt1)
|
||||
a, b = abs(d0), abs(d1)
|
||||
return (a*a + b*b) / (a+b)
|
||||
x0 = _dot(d,d0) / origDist
|
||||
x1 = _dot(d,d1) / origDist
|
||||
return (a * a + b * b) / (a + b)
|
||||
x0 = _dot(d, d0) / origDist
|
||||
x1 = _dot(d, d1) / origDist
|
||||
Len = abs(2 * (_intSecAtan(x1) - _intSecAtan(x0)) * origDist / (scale * (x1 - x0)))
|
||||
return Len
|
||||
|
||||
@ -202,13 +210,17 @@ def approximateQuadraticArcLengthC(pt1, pt2, pt3):
|
||||
# to be integrated with the best-matching fifth-degree polynomial
|
||||
# approximation of it.
|
||||
#
|
||||
#https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Legendre_quadrature
|
||||
# https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Legendre_quadrature
|
||||
|
||||
# abs(BezierCurveC[2].diff(t).subs({t:T})) for T in sorted(.5, .5±sqrt(3/5)/2),
|
||||
# weighted 5/18, 8/18, 5/18 respectively.
|
||||
v0 = abs(-0.492943519233745*pt1 + 0.430331482911935*pt2 + 0.0626120363218102*pt3)
|
||||
v1 = abs(pt3-pt1)*0.4444444444444444
|
||||
v2 = abs(-0.0626120363218102*pt1 - 0.430331482911935*pt2 + 0.492943519233745*pt3)
|
||||
v0 = abs(
|
||||
-0.492943519233745 * pt1 + 0.430331482911935 * pt2 + 0.0626120363218102 * pt3
|
||||
)
|
||||
v1 = abs(pt3 - pt1) * 0.4444444444444444
|
||||
v2 = abs(
|
||||
-0.0626120363218102 * pt1 - 0.430331482911935 * pt2 + 0.492943519233745 * pt3
|
||||
)
|
||||
|
||||
return v0 + v1 + v2
|
||||
|
||||
@ -232,14 +244,18 @@ def calcQuadraticBounds(pt1, pt2, pt3):
|
||||
(0.0, 0.0, 100, 100)
|
||||
"""
|
||||
(ax, ay), (bx, by), (cx, cy) = calcQuadraticParameters(pt1, pt2, pt3)
|
||||
ax2 = ax*2.0
|
||||
ay2 = ay*2.0
|
||||
ax2 = ax * 2.0
|
||||
ay2 = ay * 2.0
|
||||
roots = []
|
||||
if ax2 != 0:
|
||||
roots.append(-bx/ax2)
|
||||
roots.append(-bx / ax2)
|
||||
if ay2 != 0:
|
||||
roots.append(-by/ay2)
|
||||
points = [(ax*t*t + bx*t + cx, ay*t*t + by*t + cy) for t in roots if 0 <= t < 1] + [pt1, pt3]
|
||||
roots.append(-by / ay2)
|
||||
points = [
|
||||
(ax * t * t + bx * t + cx, ay * t * t + by * t + cy)
|
||||
for t in roots
|
||||
if 0 <= t < 1
|
||||
] + [pt1, pt3]
|
||||
return calcBounds(points)
|
||||
|
||||
|
||||
@ -268,7 +284,9 @@ def approximateCubicArcLength(pt1, pt2, pt3, pt4):
|
||||
>>> approximateCubicArcLength((0, 0), (50, 0), (100, -50), (-50, 0)) # cusp
|
||||
154.80848416537057
|
||||
"""
|
||||
return approximateCubicArcLengthC(complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4))
|
||||
return approximateCubicArcLengthC(
|
||||
complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4)
|
||||
)
|
||||
|
||||
|
||||
def approximateCubicArcLengthC(pt1, pt2, pt3, pt4):
|
||||
@ -288,11 +306,21 @@ def approximateCubicArcLengthC(pt1, pt2, pt3, pt4):
|
||||
|
||||
# abs(BezierCurveC[3].diff(t).subs({t:T})) for T in sorted(0, .5±(3/7)**.5/2, .5, 1),
|
||||
# weighted 1/20, 49/180, 32/90, 49/180, 1/20 respectively.
|
||||
v0 = abs(pt2-pt1)*.15
|
||||
v1 = abs(-0.558983582205757*pt1 + 0.325650248872424*pt2 + 0.208983582205757*pt3 + 0.024349751127576*pt4)
|
||||
v2 = abs(pt4-pt1+pt3-pt2)*0.26666666666666666
|
||||
v3 = abs(-0.024349751127576*pt1 - 0.208983582205757*pt2 - 0.325650248872424*pt3 + 0.558983582205757*pt4)
|
||||
v4 = abs(pt4-pt3)*.15
|
||||
v0 = abs(pt2 - pt1) * 0.15
|
||||
v1 = abs(
|
||||
-0.558983582205757 * pt1
|
||||
+ 0.325650248872424 * pt2
|
||||
+ 0.208983582205757 * pt3
|
||||
+ 0.024349751127576 * pt4
|
||||
)
|
||||
v2 = abs(pt4 - pt1 + pt3 - pt2) * 0.26666666666666666
|
||||
v3 = abs(
|
||||
-0.024349751127576 * pt1
|
||||
- 0.208983582205757 * pt2
|
||||
- 0.325650248872424 * pt3
|
||||
+ 0.558983582205757 * pt4
|
||||
)
|
||||
v4 = abs(pt4 - pt3) * 0.15
|
||||
|
||||
return v0 + v1 + v2 + v3 + v4
|
||||
|
||||
@ -325,7 +353,13 @@ def calcCubicBounds(pt1, pt2, pt3, pt4):
|
||||
yRoots = [t for t in solveQuadratic(ay3, by2, cy) if 0 <= t < 1]
|
||||
roots = xRoots + yRoots
|
||||
|
||||
points = [(ax*t*t*t + bx*t*t + cx * t + dx, ay*t*t*t + by*t*t + cy * t + dy) for t in roots] + [pt1, pt4]
|
||||
points = [
|
||||
(
|
||||
ax * t * t * t + bx * t * t + cx * t + dx,
|
||||
ay * t * t * t + by * t * t + cy * t + dy,
|
||||
)
|
||||
for t in roots
|
||||
] + [pt1, pt4]
|
||||
return calcBounds(points)
|
||||
|
||||
|
||||
@ -368,8 +402,8 @@ def splitLine(pt1, pt2, where, isHorizontal):
|
||||
pt1x, pt1y = pt1
|
||||
pt2x, pt2y = pt2
|
||||
|
||||
ax = (pt2x - pt1x)
|
||||
ay = (pt2y - pt1y)
|
||||
ax = pt2x - pt1x
|
||||
ay = pt2y - pt1y
|
||||
|
||||
bx = pt1x
|
||||
by = pt1y
|
||||
@ -422,8 +456,9 @@ def splitQuadratic(pt1, pt2, pt3, where, isHorizontal):
|
||||
((50, 50), (75, 50), (100, 0))
|
||||
"""
|
||||
a, b, c = calcQuadraticParameters(pt1, pt2, pt3)
|
||||
solutions = solveQuadratic(a[isHorizontal], b[isHorizontal],
|
||||
c[isHorizontal] - where)
|
||||
solutions = solveQuadratic(
|
||||
a[isHorizontal], b[isHorizontal], c[isHorizontal] - where
|
||||
)
|
||||
solutions = sorted([t for t in solutions if 0 <= t < 1])
|
||||
if not solutions:
|
||||
return [(pt1, pt2, pt3)]
|
||||
@ -458,8 +493,9 @@ def splitCubic(pt1, pt2, pt3, pt4, where, isHorizontal):
|
||||
((92.5259, 25), (95.202, 17.5085), (97.7062, 9.17517), (100, 1.77636e-15))
|
||||
"""
|
||||
a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)
|
||||
solutions = solveCubic(a[isHorizontal], b[isHorizontal], c[isHorizontal],
|
||||
d[isHorizontal] - where)
|
||||
solutions = solveCubic(
|
||||
a[isHorizontal], b[isHorizontal], c[isHorizontal], d[isHorizontal] - where
|
||||
)
|
||||
solutions = sorted([t for t in solutions if 0 <= t < 1])
|
||||
if not solutions:
|
||||
return [(pt1, pt2, pt3, pt4)]
|
||||
@ -524,17 +560,17 @@ def _splitQuadraticAtT(a, b, c, *ts):
|
||||
cx, cy = c
|
||||
for i in range(len(ts) - 1):
|
||||
t1 = ts[i]
|
||||
t2 = ts[i+1]
|
||||
delta = (t2 - t1)
|
||||
t2 = ts[i + 1]
|
||||
delta = t2 - t1
|
||||
# calc new a, b and c
|
||||
delta_2 = delta*delta
|
||||
delta_2 = delta * delta
|
||||
a1x = ax * delta_2
|
||||
a1y = ay * delta_2
|
||||
b1x = (2*ax*t1 + bx) * delta
|
||||
b1y = (2*ay*t1 + by) * delta
|
||||
t1_2 = t1*t1
|
||||
c1x = ax*t1_2 + bx*t1 + cx
|
||||
c1y = ay*t1_2 + by*t1 + cy
|
||||
b1x = (2 * ax * t1 + bx) * delta
|
||||
b1y = (2 * ay * t1 + by) * delta
|
||||
t1_2 = t1 * t1
|
||||
c1x = ax * t1_2 + bx * t1 + cx
|
||||
c1y = ay * t1_2 + by * t1 + cy
|
||||
|
||||
pt1, pt2, pt3 = calcQuadraticPoints((a1x, a1y), (b1x, b1y), (c1x, c1y))
|
||||
segments.append((pt1, pt2, pt3))
|
||||
@ -552,24 +588,26 @@ def _splitCubicAtT(a, b, c, d, *ts):
|
||||
dx, dy = d
|
||||
for i in range(len(ts) - 1):
|
||||
t1 = ts[i]
|
||||
t2 = ts[i+1]
|
||||
delta = (t2 - t1)
|
||||
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
|
||||
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
|
||||
a1x = ax * delta_3
|
||||
a1y = ay * delta_3
|
||||
b1x = (3*ax*t1 + bx) * delta_2
|
||||
b1y = (3*ay*t1 + by) * delta_2
|
||||
c1x = (2*bx*t1 + cx + 3*ax*t1_2) * delta
|
||||
c1y = (2*by*t1 + cy + 3*ay*t1_2) * delta
|
||||
d1x = ax*t1_3 + bx*t1_2 + cx*t1 + dx
|
||||
d1y = ay*t1_3 + by*t1_2 + cy*t1 + dy
|
||||
pt1, pt2, pt3, pt4 = calcCubicPoints((a1x, a1y), (b1x, b1y), (c1x, c1y), (d1x, d1y))
|
||||
b1x = (3 * ax * t1 + bx) * delta_2
|
||||
b1y = (3 * ay * t1 + by) * delta_2
|
||||
c1x = (2 * bx * t1 + cx + 3 * ax * t1_2) * delta
|
||||
c1y = (2 * by * t1 + cy + 3 * ay * t1_2) * delta
|
||||
d1x = ax * t1_3 + bx * t1_2 + cx * t1 + dx
|
||||
d1y = ay * t1_3 + by * t1_2 + cy * t1 + dy
|
||||
pt1, pt2, pt3, pt4 = calcCubicPoints(
|
||||
(a1x, a1y), (b1x, b1y), (c1x, c1y), (d1x, d1y)
|
||||
)
|
||||
segments.append((pt1, pt2, pt3, pt4))
|
||||
return segments
|
||||
|
||||
@ -581,8 +619,7 @@ def _splitCubicAtT(a, b, c, d, *ts):
|
||||
from math import sqrt, acos, cos, pi
|
||||
|
||||
|
||||
def solveQuadratic(a, b, c,
|
||||
sqrt=sqrt):
|
||||
def solveQuadratic(a, b, c, sqrt=sqrt):
|
||||
"""Solve a quadratic equation.
|
||||
|
||||
Solves *a*x*x + b*x + c = 0* where a, b and c are real.
|
||||
@ -602,13 +639,13 @@ def solveQuadratic(a, b, c,
|
||||
roots = []
|
||||
else:
|
||||
# We have a linear equation with 1 root.
|
||||
roots = [-c/b]
|
||||
roots = [-c / b]
|
||||
else:
|
||||
# We have a true quadratic equation. Apply the quadratic formula to find two roots.
|
||||
DD = b*b - 4.0*a*c
|
||||
DD = b * b - 4.0 * a * c
|
||||
if DD >= 0.0:
|
||||
rDD = sqrt(DD)
|
||||
roots = [(-b+rDD)/2.0/a, (-b-rDD)/2.0/a]
|
||||
roots = [(-b + rDD) / 2.0 / a, (-b - rDD) / 2.0 / a]
|
||||
else:
|
||||
# complex roots, ignore
|
||||
roots = []
|
||||
@ -658,52 +695,52 @@ def solveCubic(a, b, c, d):
|
||||
# returns unreliable results, so we fall back to quad.
|
||||
return solveQuadratic(b, c, d)
|
||||
a = float(a)
|
||||
a1 = b/a
|
||||
a2 = c/a
|
||||
a3 = d/a
|
||||
a1 = b / a
|
||||
a2 = c / a
|
||||
a3 = d / a
|
||||
|
||||
Q = (a1*a1 - 3.0*a2)/9.0
|
||||
R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0
|
||||
Q = (a1 * a1 - 3.0 * a2) / 9.0
|
||||
R = (2.0 * a1 * a1 * a1 - 9.0 * a1 * a2 + 27.0 * a3) / 54.0
|
||||
|
||||
R2 = R*R
|
||||
Q3 = Q*Q*Q
|
||||
R2 = R * R
|
||||
Q3 = Q * Q * Q
|
||||
R2 = 0 if R2 < epsilon else R2
|
||||
Q3 = 0 if abs(Q3) < epsilon else Q3
|
||||
|
||||
R2_Q3 = R2 - Q3
|
||||
|
||||
if R2 == 0. and Q3 == 0.:
|
||||
x = round(-a1/3.0, epsilonDigits)
|
||||
if R2 == 0.0 and Q3 == 0.0:
|
||||
x = round(-a1 / 3.0, epsilonDigits)
|
||||
return [x, x, x]
|
||||
elif R2_Q3 <= epsilon * .5:
|
||||
elif R2_Q3 <= epsilon * 0.5:
|
||||
# The epsilon * .5 above ensures that Q3 is not zero.
|
||||
theta = acos(max(min(R/sqrt(Q3), 1.0), -1.0))
|
||||
rQ2 = -2.0*sqrt(Q)
|
||||
a1_3 = a1/3.0
|
||||
x0 = rQ2*cos(theta/3.0) - a1_3
|
||||
x1 = rQ2*cos((theta+2.0*pi)/3.0) - a1_3
|
||||
x2 = rQ2*cos((theta+4.0*pi)/3.0) - a1_3
|
||||
theta = acos(max(min(R / sqrt(Q3), 1.0), -1.0))
|
||||
rQ2 = -2.0 * sqrt(Q)
|
||||
a1_3 = a1 / 3.0
|
||||
x0 = rQ2 * cos(theta / 3.0) - a1_3
|
||||
x1 = rQ2 * cos((theta + 2.0 * pi) / 3.0) - a1_3
|
||||
x2 = rQ2 * cos((theta + 4.0 * pi) / 3.0) - a1_3
|
||||
x0, x1, x2 = sorted([x0, x1, x2])
|
||||
# Merge roots that are close-enough
|
||||
if x1 - x0 < epsilon and x2 - x1 < epsilon:
|
||||
x0 = x1 = x2 = round((x0 + x1 + x2) / 3., epsilonDigits)
|
||||
x0 = x1 = x2 = round((x0 + x1 + x2) / 3.0, epsilonDigits)
|
||||
elif x1 - x0 < epsilon:
|
||||
x0 = x1 = round((x0 + x1) / 2., epsilonDigits)
|
||||
x0 = x1 = round((x0 + x1) / 2.0, epsilonDigits)
|
||||
x2 = round(x2, epsilonDigits)
|
||||
elif x2 - x1 < epsilon:
|
||||
x0 = round(x0, epsilonDigits)
|
||||
x1 = x2 = round((x1 + x2) / 2., epsilonDigits)
|
||||
x1 = x2 = round((x1 + x2) / 2.0, epsilonDigits)
|
||||
else:
|
||||
x0 = round(x0, epsilonDigits)
|
||||
x1 = round(x1, epsilonDigits)
|
||||
x2 = round(x2, epsilonDigits)
|
||||
return [x0, x1, x2]
|
||||
else:
|
||||
x = pow(sqrt(R2_Q3)+abs(R), 1/3.0)
|
||||
x = x + Q/x
|
||||
x = pow(sqrt(R2_Q3) + abs(R), 1 / 3.0)
|
||||
x = x + Q / x
|
||||
if R >= 0.0:
|
||||
x = -x
|
||||
x = round(x - a1/3.0, epsilonDigits)
|
||||
x = round(x - a1 / 3.0, epsilonDigits)
|
||||
return [x]
|
||||
|
||||
|
||||
@ -711,6 +748,7 @@ def solveCubic(a, b, c, d):
|
||||
# Conversion routines for points to parameters and vice versa
|
||||
#
|
||||
|
||||
|
||||
def calcQuadraticParameters(pt1, pt2, pt3):
|
||||
x2, y2 = pt2
|
||||
x3, y3 = pt3
|
||||
@ -1176,7 +1214,9 @@ def printSegments(segments):
|
||||
for segment in segments:
|
||||
print(_segmentrepr(segment))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
import sys
|
||||
import doctest
|
||||
|
||||
sys.exit(doctest.testmod().failed)
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user