diff --git a/BezierKit/BezierKit.xcodeproj/project.pbxproj b/BezierKit/BezierKit.xcodeproj/project.pbxproj index ae556ac..19f1556 100644 --- a/BezierKit/BezierKit.xcodeproj/project.pbxproj +++ b/BezierKit/BezierKit.xcodeproj/project.pbxproj @@ -49,6 +49,10 @@ FD84360622B0091500AA90EF /* PathComponent+WindingCount.swift in Sources */ = {isa = PBXBuildFile; fileRef = FD84360422B0091500AA90EF /* PathComponent+WindingCount.swift */; }; FDA727591ED5035300011871 /* CubicCurveTests.swift in Sources */ = {isa = PBXBuildFile; fileRef = FDA727581ED5035300011871 /* CubicCurveTests.swift */; }; FDA7275A1ED5035300011871 /* CubicCurveTests.swift in Sources */ = {isa = PBXBuildFile; fileRef = FDA727581ED5035300011871 /* CubicCurveTests.swift */; }; + FDB6011B25BB9B3700BAB067 /* BezierCurve+Polynomial.swift in Sources */ = {isa = PBXBuildFile; fileRef = FDB6011A25BB9B3700BAB067 /* BezierCurve+Polynomial.swift */; }; + FDB6011C25BB9B3700BAB067 /* BezierCurve+Polynomial.swift in Sources */ = {isa = PBXBuildFile; fileRef = FDB6011A25BB9B3700BAB067 /* BezierCurve+Polynomial.swift */; }; + FDB6012225BBA06600BAB067 /* BezierCurve+PolynomialTests.swift in Sources */ = {isa = PBXBuildFile; fileRef = FDB6012125BBA06600BAB067 /* BezierCurve+PolynomialTests.swift */; }; + FDB6012325BBA06600BAB067 /* BezierCurve+PolynomialTests.swift in Sources */ = {isa = PBXBuildFile; fileRef = FDB6012125BBA06600BAB067 /* BezierCurve+PolynomialTests.swift */; }; FDB6B4021EAFD6DF00001C61 /* BezierCurve.swift in Sources */ = {isa = PBXBuildFile; fileRef = FDB6B3F71EAFD6DF00001C61 /* BezierCurve.swift */; }; FDB6B4031EAFD6DF00001C61 /* BezierCurve.swift in Sources */ = {isa = PBXBuildFile; fileRef = FDB6B3F71EAFD6DF00001C61 /* BezierCurve.swift */; }; FDB6B4041EAFD6DF00001C61 /* CubicCurve.swift in Sources */ = {isa = PBXBuildFile; fileRef = FDB6B3F81EAFD6DF00001C61 /* CubicCurve.swift */; }; @@ -170,6 +174,8 @@ FD80FF8922B1CEE00018C592 /* LockTests.swift */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.swift; path = LockTests.swift; sourceTree = ""; }; FD84360422B0091500AA90EF /* PathComponent+WindingCount.swift */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.swift; path = "PathComponent+WindingCount.swift"; sourceTree = ""; }; FDA727581ED5035300011871 /* CubicCurveTests.swift */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.swift; path = CubicCurveTests.swift; sourceTree = ""; }; + FDB6011A25BB9B3700BAB067 /* BezierCurve+Polynomial.swift */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.swift; path = "BezierCurve+Polynomial.swift"; sourceTree = ""; }; + FDB6012125BBA06600BAB067 /* BezierCurve+PolynomialTests.swift */ = {isa = PBXFileReference; lastKnownFileType = sourcecode.swift; path = "BezierCurve+PolynomialTests.swift"; sourceTree = ""; }; FDB6B3F71EAFD6DF00001C61 /* BezierCurve.swift */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.swift; path = BezierCurve.swift; sourceTree = ""; }; FDB6B3F81EAFD6DF00001C61 /* CubicCurve.swift */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.swift; path = CubicCurve.swift; sourceTree = ""; }; FDB6B3F91EAFD6DF00001C61 /* Draw.swift */ = {isa = PBXFileReference; fileEncoding = 4; lastKnownFileType = sourcecode.swift; path = Draw.swift; sourceTree = ""; }; @@ -279,6 +285,7 @@ FD0F550A1DC43FFB0084CDCD /* Info.plist */, FD0F55081DC43FFB0084CDCD /* BezierKitTestHelpers.swift */, FDF0664D1FFA0C9900123308 /* BezierCurveTests.swift */, + FDB6012125BBA06600BAB067 /* BezierCurve+PolynomialTests.swift */, FDE6CD8E1EC8F9BD00FAB479 /* LineSegmentTests.swift */, FDA727581ED5035300011871 /* CubicCurveTests.swift */, FD40244F2110CF5100FA723C /* QuadraticCurveTests.swift */, @@ -326,6 +333,7 @@ FD149EB92135CBFF009E791D /* AugmentedGraph.swift */, FDC859622119274A00AF7642 /* BoundingBoxHierarchy.swift */, FDB6B3F71EAFD6DF00001C61 /* BezierCurve.swift */, + FDB6011A25BB9B3700BAB067 /* BezierCurve+Polynomial.swift */, FD5CF14B22400FCA00FE15A6 /* BezierCurve+Intersection.swift */, FDB6B3F81EAFD6DF00001C61 /* CubicCurve.swift */, FDB6B3FC1EAFD6DF00001C61 /* CGPoint+Overloads.swift */, @@ -657,6 +665,7 @@ FDB6B40F1EAFD6DF00001C61 /* PathComponent.swift in Sources */, FDB6B4151EAFD6DF00001C61 /* Utils.swift in Sources */, FDCE99A6223C404E00597989 /* Path+Data.swift in Sources */, + FDB6011C25BB9B3700BAB067 /* BezierCurve+Polynomial.swift in Sources */, FDC2EB4B2298735C007768FC /* Lock.swift in Sources */, FDB6B4071EAFD6DF00001C61 /* Draw.swift in Sources */, FDB6B4131EAFD6DF00001C61 /* Types.swift in Sources */, @@ -683,6 +692,7 @@ FDB6B40E1EAFD6DF00001C61 /* PathComponent.swift in Sources */, FDB6B4141EAFD6DF00001C61 /* Utils.swift in Sources */, FDCE99A5223C404E00597989 /* Path+Data.swift in Sources */, + FDB6011B25BB9B3700BAB067 /* BezierCurve+Polynomial.swift in Sources */, FDC2EB4A2298735C007768FC /* Lock.swift in Sources */, FDB6B4061EAFD6DF00001C61 /* Draw.swift in Sources */, FDB6B4121EAFD6DF00001C61 /* Types.swift in Sources */, @@ -708,6 +718,7 @@ FD12F1E22288CF6900404CE1 /* UtilsTests.swift in Sources */, FDA727591ED5035300011871 /* CubicCurveTests.swift in Sources */, FDB9D7331EB28CEB00413F0E /* BezierKit_iOSTests.swift in Sources */, + FDB6012225BBA06600BAB067 /* BezierCurve+PolynomialTests.swift in Sources */, FD80FF8A22B1CEE00018C592 /* LockTests.swift in Sources */, FDE6CD8F1EC8F9BD00FAB479 /* LineSegmentTests.swift in Sources */, FD4024502110CF5100FA723C /* QuadraticCurveTests.swift in Sources */, @@ -737,6 +748,7 @@ FD12F1E32288CF6900404CE1 /* UtilsTests.swift in Sources */, FDA7275A1ED5035300011871 /* CubicCurveTests.swift in Sources */, FDB9D7421EB28D1900413F0E /* BezierKit_MacTests.swift in Sources */, + FDB6012325BBA06600BAB067 /* BezierCurve+PolynomialTests.swift in Sources */, FD80FF8B22B1CEE00018C592 /* LockTests.swift in Sources */, FDE6CD901EC8F9BD00FAB479 /* LineSegmentTests.swift in Sources */, FD4024512110CF5100FA723C /* QuadraticCurveTests.swift in Sources */, diff --git a/BezierKit/BezierKitTests/BezierCurve+PolynomialTests.swift b/BezierKit/BezierKitTests/BezierCurve+PolynomialTests.swift new file mode 100644 index 0000000..70ee666 --- /dev/null +++ b/BezierKit/BezierKitTests/BezierCurve+PolynomialTests.swift @@ -0,0 +1,78 @@ +// +// BezierCurve+PolynomialTests.swift +// BezierKit +// +// Created by Holmes Futrell on 1/22/21. +// Copyright © 2021 Holmes Futrell. All rights reserved. +// + +import BezierKit +import XCTest + +class BezierCurve_PolynomialTests: XCTestCase { + + func testPolynomialLineSegment() { + let lineSegment = LineSegment(p0: CGPoint(x: 3, y: 4), p1: CGPoint(x: 5, y: 6)) + XCTAssertEqual(lineSegment.xPolynomial, BernsteinPolynomial1(b0: 3, b1: 5)) + XCTAssertEqual(lineSegment.yPolynomial, BernsteinPolynomial1(b0: 4, b1: 6)) + } + + func testPolynomialQuadratic() { + let quadratic = QuadraticCurve(p0: CGPoint(x: 1, y: 0), + p1: CGPoint(x: 2, y: -2), + p2: CGPoint(x: 3, y: -1)) + XCTAssertEqual(quadratic.xPolynomial, BernsteinPolynomial2(b0: 1, b1: 2, b2: 3)) + XCTAssertEqual(quadratic.yPolynomial, BernsteinPolynomial2(b0: 0, b1: -2, b2: -1)) + } + + func testPolynomialCubic() { + let cubic = CubicCurve(p0: CGPoint(x: 1, y: 0), + p1: CGPoint(x: 2, y: 2), + p2: CGPoint(x: 3, y: 1), + p3: CGPoint(x: 4, y: -1)) + XCTAssertEqual(cubic.xPolynomial, BernsteinPolynomial3(b0: 1, b1: 2, b2: 3, b3: 4)) + XCTAssertEqual(cubic.yPolynomial, BernsteinPolynomial3(b0: 0, b1: 2, b2: 1, b3: -1)) + } + + func testExtremaLine() { + let l1 = LineSegment(p0: CGPoint(x: 1.0, y: 2.0), p1: CGPoint(x: 4.0, y: 6.0)) + let (x1, y1, all1) = l1.extrema() + XCTAssertTrue(x1.isEmpty) + XCTAssertTrue(y1.isEmpty) + XCTAssertTrue(all1.isEmpty) + + let l2 = LineSegment(p0: CGPoint(x: 1.0, y: 2.0), p1: CGPoint(x: 4.0, y: 2.0)) + let (x2, y2, all2) = l2.extrema() + XCTAssertTrue(x2.isEmpty) + XCTAssertTrue(y2.isEmpty) + XCTAssertTrue(all2.isEmpty) + } + + func testExtremaQuadratic() { + let f: [CGFloat] = [4, -2, 1] // f(t) = 4t^2 - 2t + 1, which has a local minimum at t = 0.25 + let g: [CGFloat] = [1, -4, 4] // g(t) = t^2 -4t + 4, which has a local minimum at t = 2 (outside parameter range) + let q = BezierKitTestHelpers.quadraticCurveFromPolynomials(f, g) + let (x, y, all) = q.extrema() + XCTAssertEqual(all.count, 1) + XCTAssertEqual(all[0], 0.25) + XCTAssertEqual(x.count, 1) + XCTAssertEqual(x[0], 0.25) + XCTAssertTrue(y.isEmpty) + } + + func testExtremaCubic() { + let f: [CGFloat] = [1, -1, 0, 0] // f(t) = t^3 - t^2, which has two local minimum at t=0, t=2/3 and an inflection point t=1/3 + let g: [CGFloat] = [0, 3, -2, 0] // g(t) = 3t^2 - 2t, which has a local minimum at t=1/3 + let c = BezierKitTestHelpers.cubicCurveFromPolynomials(f, g) + let (x, y, all) = c.extrema() + XCTAssertEqual(all.count, 3) + XCTAssertEqual(all[0], 0.0) + XCTAssertEqual(all[1], 1.0 / 3.0) + XCTAssertEqual(all[2], 2.0 / 3.0) + XCTAssertEqual(x[0], 0.0) + XCTAssertEqual(x[1], 1.0 / 3.0) + XCTAssertEqual(x[2], 2.0 / 3.0) + XCTAssertEqual(y.count, 1) + XCTAssertEqual(y[0], 1.0 / 3.0) + } +} diff --git a/BezierKit/BezierKitTests/BezierKitTestHelpers.swift b/BezierKit/BezierKitTests/BezierKitTestHelpers.swift index e4faaa4..c126e2d 100644 --- a/BezierKit/BezierKitTests/BezierKitTestHelpers.swift +++ b/BezierKit/BezierKitTests/BezierKitTestHelpers.swift @@ -76,6 +76,14 @@ class BezierKitTestHelpers { return sum } + static func quadraticCurveFromPolynomials(_ f: [CGFloat], _ g: [CGFloat]) -> QuadraticCurve { + precondition(f.count == 3 && g.count == 3) + let curve = QuadraticCurve(p0: CGPoint(x: f[2], y: g[2]), + p1: CGPoint(x: 0.5 * f[1] + f[2], y: 0.5 * g[1] + g[2]), + p2: CGPoint(x: f[0] + f[1] + f[2], y: g[0] + g[1] + g[2])) + return curve + } + static func cubicCurveFromPolynomials(_ f: [CGFloat], _ g: [CGFloat]) -> CubicCurve { precondition(f.count == 4 && g.count == 4) // create a cubic bezier curve from two polynomials @@ -89,12 +97,7 @@ class BezierKitTestHelpers { let b = r / 3.0 + a let c = q / 3.0 + 2.0 * b - a let d = p + a - 3.0 * b + 3.0 * c - // check that it worked - let curve = CubicCurve(p0: a, p1: b, p2: c, p3: d) - for t: CGFloat in stride(from: 0, through: 1, by: 0.1) { - assert(distance(curve.point(at: t), CGPoint(x: evaluatePolynomial(f, at: t), y: evaluatePolynomial(g, at: t))) < 0.001, "internal error! failed to fit polynomial!") - } - return curve + return CubicCurve(p0: a, p1: b, p2: c, p3: d) } static func isSatisfactoryReduceResult(_ result: [Subcurve], for curve: A) -> Bool { diff --git a/BezierKit/BezierKitTests/CubicCurveTests.swift b/BezierKit/BezierKitTests/CubicCurveTests.swift index 203851c..bc1c7dd 100644 --- a/BezierKit/BezierKitTests/CubicCurveTests.swift +++ b/BezierKit/BezierKitTests/CubicCurveTests.swift @@ -250,22 +250,6 @@ class CubicCurveTests: XCTestCase { XCTAssertEqual(c1.length(), 5.0, accuracy: epsilon) } - func testExtrema() { - let f: [CGFloat] = [1, -1, 0, 0] // f(t) = t^3 - t^2, which has two local minimum at t=0, t=2/3 and an inflection point t=1/3 - let g: [CGFloat] = [0, 3, -2, 0] // g(t) = 3t^2 - 2t, which has a local minimum at t=1/3 - let c = BezierKitTestHelpers.cubicCurveFromPolynomials(f, g) - let (x, y, all) = c.extrema() - XCTAssertEqual(all.count, 3) - XCTAssertEqual(all[0], 0.0) - XCTAssertEqual(all[1], 1.0 / 3.0) - XCTAssertEqual(all[2], 2.0 / 3.0) - XCTAssertEqual(x[0], 0.0) - XCTAssertEqual(x[1], 1.0 / 3.0) - XCTAssertEqual(x[2], 2.0 / 3.0) - XCTAssertEqual(y.count, 1) - XCTAssertEqual(y[0], 1.0 / 3.0) - } - func testProject() { let epsilon: CGFloat = 1.0e-5 // test a cubic diff --git a/BezierKit/BezierKitTests/LineSegmentTests.swift b/BezierKit/BezierKitTests/LineSegmentTests.swift index a9d53cc..79fff7c 100644 --- a/BezierKit/BezierKitTests/LineSegmentTests.swift +++ b/BezierKit/BezierKitTests/LineSegmentTests.swift @@ -100,14 +100,6 @@ class LineSegmentTests: XCTestCase { XCTAssertEqual(l.length(), 5.0) } - func testExtrema() { - let l = LineSegment(p0: CGPoint(x: 1.0, y: 2.0), p1: CGPoint(x: 4.0, y: 6.0)) - let (x, y, all) = l.extrema() - XCTAssertTrue(x.isEmpty) - XCTAssertTrue(y.isEmpty) - XCTAssertTrue(all.isEmpty) - } - func testProject() { let l = LineSegment(p0: CGPoint(x: 1.0, y: 2.0), p1: CGPoint(x: 5.0, y: 6.0)) let p1 = l.project(CGPoint(x: 0.0, y: 0.0)) // should project to p0 diff --git a/BezierKit/BezierKitTests/PolynomialTests.swift b/BezierKit/BezierKitTests/PolynomialTests.swift index 0ff51be..1681921 100644 --- a/BezierKit/BezierKitTests/PolynomialTests.swift +++ b/BezierKit/BezierKitTests/PolynomialTests.swift @@ -11,45 +11,45 @@ import XCTest class PolynomialTests: XCTestCase { - let accuracy = 1.0e-5 + let accuracy: CGFloat = 1.0e-5 func testEvaluation() { let point = BernsteinPolynomial0(b0: 3.0) XCTAssertEqual(point.reduce(a1: 1, a2: 2), 0) - XCTAssertEqual(point.f(0), 3) - XCTAssertEqual(point.f(0.5), 3) - XCTAssertEqual(point.f(1), 3) + XCTAssertEqual(point.value(at: 0), 3) + XCTAssertEqual(point.value(at: 0.5), 3) + XCTAssertEqual(point.value(at: 1), 3) XCTAssertEqual(point.derivative, BernsteinPolynomial0(b0: 0.0)) - XCTAssertEqual(point.analyticalRoots(between: 0, and: 1), []) + XCTAssertEqual(point.distinctAnalyticalRoots(between: 0, and: 1), []) XCTAssertEqual(point.coefficients, [3.0]) let line = BernsteinPolynomial1(b0: 2.0, b1: 4.0) - XCTAssertEqual(line.f(0), 2) - XCTAssertEqual(line.f(0.5), 3) - XCTAssertEqual(line.f(1), 4) + XCTAssertEqual(line.value(at: 0), 2) + XCTAssertEqual(line.value(at: 0.5), 3) + XCTAssertEqual(line.value(at: 1), 4) XCTAssertEqual(line.derivative, BernsteinPolynomial0(b0: 2)) - XCTAssertEqual(line.analyticalRoots(between: -2, and: 1), [-1]) - XCTAssertEqual(line.analyticalRoots(between: 0, and: 1), []) + XCTAssertEqual(line.distinctAnalyticalRoots(between: -2, and: 1), [-1]) + XCTAssertEqual(line.distinctAnalyticalRoots(between: 0, and: 1), []) XCTAssertEqual(line.coefficients, [2, 4]) let quad = BernsteinPolynomial2(b0: -1, b1: 1.0, b2: 0.0) - XCTAssertEqual(quad.f(0), -1) - XCTAssertEqual(quad.f(0.5), 0.25) - XCTAssertEqual(quad.f(1), 0) + XCTAssertEqual(quad.value(at: 0), -1) + XCTAssertEqual(quad.value(at: 0.5), 0.25) + XCTAssertEqual(quad.value(at: 1), 0) XCTAssertEqual(quad.derivative, BernsteinPolynomial1(b0: 4, b1: -2)) XCTAssertEqual(quad.coefficients, [-1, 1, 0]) } func testDegree1() { let polynomial = BernsteinPolynomial1(b0: -3, b1: 2) - let roots = findRoots(of: polynomial, between: -1, and: 1) + let roots = findDistinctRoots(of: polynomial, between: -1, and: 1) XCTAssertEqual(roots.count, 1) - XCTAssertEqual(roots[0], 0.6, accuracy: accuracy) + XCTAssertEqual(roots[0], CGFloat(0.6), accuracy: accuracy) } func testDegree2() { let polynomial = BernsteinPolynomial2(b0: -5, b1: -6, b2: -4) - let roots = findRoots(of: polynomial, between: -10, and: 10) + let roots = findDistinctRoots(of: polynomial, between: -10, and: 10) XCTAssertEqual(roots[0], -1, accuracy: accuracy) XCTAssertEqual(roots[1], 1.0 + 2.0 / 3.0, accuracy: accuracy) } @@ -57,27 +57,37 @@ class PolynomialTests: XCTestCase { func testDegree3() { // x^3 - 6x^2 + 11x - 6 let polynomial = BernsteinPolynomial3(b0: -6, b1: -7.0 / 3.0, b2: -2.0 / 3.0, b3: 0) - XCTAssertEqual(polynomial.coefficients, [-6, -7.0 / 3.0, -2.0 / 3.0, 0.0]) - let roots = findRoots(of: polynomial, between: 0, and: 4) + XCTAssertEqual(polynomial.coefficients, [-6, CGFloat(-7.0 / 3.0), CGFloat(-2.0 / 3.0), 0.0]) + let roots = findDistinctRoots(of: polynomial, between: 0, and: 4) XCTAssertEqual(roots[0], 1, accuracy: accuracy) XCTAssertEqual(roots[1], 2, accuracy: accuracy) XCTAssertEqual(roots[2], 3, accuracy: accuracy) } - func testDegree3RepeatedRoot() { + func testDegree3RepeatedRoot1() { // x^3 - 4x^2 + 5x - 2 // repeated root at x = 1 let polynomial = BernsteinPolynomial3(b0: -2, b1: -1.0 / 3.0, b2: 0, b3: 0) - let roots = findRoots(of: polynomial, between: -1, and: 3) + let roots = findDistinctRoots(of: polynomial, between: -1, and: 3) XCTAssertEqual(roots[0], 1, accuracy: accuracy) XCTAssertEqual(roots[1], 2, accuracy: accuracy) } +// func testDegree3RootExactlyZero() { +// // root is exactly t = 0 (at the start of unit interval), +// // so may be accidentally discarded due to numerical precision +// let polynomial = BernsteinPolynomial3(b0: 0, b1: 96, b2: -24, b3: -36) +// let roots = findRoots(of: polynomial, between: 0, and: 1) +// XCTAssertEqual(roots.count, 2) +// XCTAssertEqual(roots[0], 0.0) +// XCTAssertEqual(roots[1], 2.0 / 3.0, accuracy: accuracy) +// } + func testDegree4() { // x^4 - 2.44x^2 + 1.44 - let polynomial = BernsteinPolynomial4(b0: 1.44, b1: 1.44, b2: 1.44 - 1.22 / 3, b3: 0.22, b4: 0) - XCTAssertEqual(polynomial.coefficients, [1.44, 1.44, 1.44 - 1.22 / 3, 0.22, 0]) - let roots = findRoots(of: polynomial, between: -2, and: 2) + let polynomial = BernsteinPolynomial4(b0: 1.44, b1: 1.44, b2: CGFloat(1.44 - 1.22 / 3), b3: 0.22, b4: 0) + XCTAssertEqual(polynomial.coefficients, [1.44, 1.44, CGFloat(1.44 - 1.22 / 3), 0.22, 0]) + let roots = findDistinctRoots(of: polynomial, between: -2, and: 2) XCTAssertEqual(roots[0], -1.2, accuracy: accuracy) XCTAssertEqual(roots[1], -1, accuracy: accuracy) XCTAssertEqual(roots[2], 1, accuracy: accuracy) @@ -87,7 +97,7 @@ class PolynomialTests: XCTestCase { func testDegree4RepeatedRoots() { // x^4 - 2x^2 + 1 let polynomial = BernsteinPolynomial4(b0: 1, b1: 1, b2: 2.0 / 3.0, b3: 0, b4: 0) - let roots = findRoots(of: polynomial, between: -2, and: 2) + let roots = findDistinctRoots(of: polynomial, between: -2, and: 2) XCTAssertEqual(roots.count, 2) XCTAssertEqual(roots[0], -1, accuracy: accuracy) XCTAssertEqual(roots[1], 1, accuracy: accuracy) @@ -97,8 +107,7 @@ class PolynomialTests: XCTestCase { // 0.2x^5 - 0.813333x^3 - 8.56x let polynomial = BernsteinPolynomial5(b0: 0, b1: -1.712, b2: -3.424, b3: -5.2173333, b4: -7.1733332, b5: -9.173333) XCTAssertEqual(polynomial.coefficients, [0, -1.712, -3.424, -5.2173333, -7.1733332, -9.173333]) - let roots = findRoots(of: polynomial, between: -4, and: 4) - XCTAssertEqual(polynomial.analyticalRoots(between: -5, and: 5), nil, "shouldn't be possible to solve analytically") + let roots = findDistinctRoots(of: polynomial, between: -4, and: 4) XCTAssertEqual(roots[0], -2.9806382, accuracy: accuracy) XCTAssertEqual(roots[1], 0, accuracy: accuracy) XCTAssertEqual(roots[2], 2.9806382, accuracy: accuracy) @@ -106,7 +115,7 @@ class PolynomialTests: XCTestCase { func testDegree4RealWorldIssue() { let polynomial = BernsteinPolynomial4(b0: 1819945.4373168945, b1: -3353335.8194732666, b2: 3712712.6330566406, b3: -2836657.1703338623, b4: 2483314.5947265625) - let roots = findRoots(of: polynomial, between: 0, and: 1) + let roots = findDistinctRootsInUnitInterval(of: polynomial) XCTAssertEqual(roots.count, 2) XCTAssertEqual(roots[0], 0.15977874432923783, accuracy: 1.0e-5) XCTAssertEqual(roots[1], 0.407811682610126, accuracy: 1.0e-5) diff --git a/BezierKit/BezierKitTests/QuadraticCurveTests.swift b/BezierKit/BezierKitTests/QuadraticCurveTests.swift index c735790..e404bad 100644 --- a/BezierKit/BezierKitTests/QuadraticCurveTests.swift +++ b/BezierKit/BezierKitTests/QuadraticCurveTests.swift @@ -100,8 +100,6 @@ class QuadraticCurveTests: XCTestCase { // func testLength() { // } // -// func testExtrema() { -// } func testProject() { let epsilon: CGFloat = 1.0e-5 diff --git a/BezierKit/Library/BezierCurve+Polynomial.swift b/BezierKit/Library/BezierCurve+Polynomial.swift new file mode 100644 index 0000000..4388b12 --- /dev/null +++ b/BezierKit/Library/BezierCurve+Polynomial.swift @@ -0,0 +1,51 @@ +// +// BezierCurve+Polynomial.swift +// BezierKit +// +// Created by Holmes Futrell on 1/22/21. +// Copyright © 2021 Holmes Futrell. All rights reserved. +// + +import Foundation + +/// a parametric function whose x and y coordinates can be considered as separate polynomial functions +/// eg `f(t) = (xPolynomial(t), yPolynomial(t))` +public protocol ComponentPolynomials { + associatedtype Polynomial: BernsteinPolynomial + var xPolynomial: Polynomial { get } + var yPolynomial: Polynomial { get } +} + +extension LineSegment: ComponentPolynomials { + public var xPolynomial: BernsteinPolynomial1 { return BernsteinPolynomial1(b0: self.p0.x, b1: self.p1.x) } + public var yPolynomial: BernsteinPolynomial1 { return BernsteinPolynomial1(b0: self.p0.y, b1: self.p1.y) } +} + +extension QuadraticCurve: ComponentPolynomials { + public var xPolynomial: BernsteinPolynomial2 { return BernsteinPolynomial2(b0: self.p0.x, b1: self.p1.x, b2: self.p2.x) } + public var yPolynomial: BernsteinPolynomial2 { return BernsteinPolynomial2(b0: self.p0.y, b1: self.p1.y, b2: self.p2.y) } +} + +extension CubicCurve: ComponentPolynomials { + public var xPolynomial: BernsteinPolynomial3 { return BernsteinPolynomial3(b0: self.p0.x, b1: self.p1.x, b2: self.p2.x, b3: self.p3.x) } + public var yPolynomial: BernsteinPolynomial3 { return BernsteinPolynomial3(b0: self.p0.y, b1: self.p1.y, b2: self.p2.y, b3: self.p3.y) } +} + +extension BezierCurve where Self: ComponentPolynomials { + /// default implementation of `extrema` by finding roots of component polynomials + public func extrema() -> (x: [CGFloat], y: [CGFloat], all: [CGFloat]) { + func rootsForPolynomial(_ polynomial: B) -> [CGFloat] { + let firstOrderDerivative = polynomial.derivative + var roots = findDistinctRootsInUnitInterval(of: firstOrderDerivative) + if self.order >= 3 { + let secondOrderDerivative = firstOrderDerivative.derivative + roots += findDistinctRootsInUnitInterval(of: secondOrderDerivative) + } + return roots.sortedAndUniqued() + } + let xRoots = rootsForPolynomial(self.xPolynomial) + let yRoots = rootsForPolynomial(self.yPolynomial) + let allRoots = (xRoots + yRoots).sortedAndUniqued() + return (x: xRoots, y: yRoots, all: allRoots) + } +} diff --git a/BezierKit/Library/BezierCurve.swift b/BezierKit/Library/BezierCurve.swift index 3530b7f..b90efc1 100644 --- a/BezierKit/Library/BezierCurve.swift +++ b/BezierKit/Library/BezierCurve.swift @@ -66,27 +66,6 @@ extension BezierCurve { return Utils.length({(_ t: CGFloat) in self.derivative(at: t)}) } - public func extrema() -> (x: [CGFloat], y: [CGFloat], all: [CGFloat]) { - func sequentialDifference(_ array: [T]) -> [T] where T: FloatingPoint { - return (1.. [CGFloat] { - let values = self.points.map { $0[dimension] } - let firstOrderDiffs = sequentialDifference(values) - var roots = Utils.droots(firstOrderDiffs) - if self.order >= 3 { - let secondOrderDiffs = sequentialDifference(firstOrderDiffs) - roots += Utils.droots(secondOrderDiffs) - } - return roots.filter({$0 >= 0 && $0 <= 1}).sortedAndUniqued() - } - guard self.order > 1 else { return (x: [], y: [], all: []) } - let xRoots = rootsForDimension(0) - let yRoots = rootsForDimension(1) - let allRoots = (xRoots + yRoots).sortedAndUniqued() - return (x: xRoots, y: yRoots, all: allRoots) - } - // MARK: - public func hull(_ t: CGFloat) -> [CGPoint] { return Utils.hull(self.points, t) diff --git a/BezierKit/Library/CubicCurve.swift b/BezierKit/Library/CubicCurve.swift index 8962ccc..46fc49b 100644 --- a/BezierKit/Library/CubicCurve.swift +++ b/BezierKit/Library/CubicCurve.swift @@ -220,14 +220,14 @@ public struct CubicCurve: NonlinearBezierCurve, Equatable { minimumDistanceSquared = lengthSquaredEnd } // the roots represent the values at which the curve and its derivative are perpendicular - // ie, the dot product of q and l is equal to zero - let points = BernsteinPolynomial5(b0: Double(p0.x + p0.y), - b1: Double(p1.x + p1.y), - b2: Double(p2.x + p2.y), - b3: Double(p3.x + p3.y), - b4: Double(p4.x + p4.y), - b5: Double(p5.x + p5.y)) - for t in findRoots(of: points, between: 0, and: 1) { + // ie, the dot product of c and q is equal to zero + let polynomial = BernsteinPolynomial5(b0: p0.x + p0.y, + b1: p1.x + p1.y, + b2: p2.x + p2.y, + b3: p3.x + p3.y, + b4: p4.x + p4.y, + b5: p5.x + p5.y) + for t in findDistinctRootsInUnitInterval(of: polynomial) { guard t > 0.0, t < 1.0 else { break } let point = c.point(at: CGFloat(t)) let distanceSquared = point.lengthSquared diff --git a/BezierKit/Library/Polynomial.swift b/BezierKit/Library/Polynomial.swift index 435ec38..284be21 100644 --- a/BezierKit/Library/Polynomial.swift +++ b/BezierKit/Library/Polynomial.swift @@ -11,51 +11,97 @@ import CoreGraphics #endif import Foundation -internal protocol BernsteinPolynomial: Equatable { - func f(_ x: Double) -> Double +public protocol BernsteinPolynomial: Equatable { + func value(at x: CGFloat) -> CGFloat var order: Int { get } - var coefficients: [Double] { get } -// var last: Double { get } -// var first: Double { get } -// func enumerated(block: (Int, Double) -> Void) - associatedtype Difference: BernsteinPolynomial - func difference(a1: Double, a2: Double) -> Difference - func reduce(a1: Double, a2: Double) -> Double - var derivative: Difference { get } -// init(_ d: Difference, last: Double) -// init(first: Double, _ d: Difference) + var coefficients: [CGFloat] { get } +// var last: CGFloat { get } +// var first: CGFloat { get } +// func enumerated(block: (Int, CGFloat) -> Void) + associatedtype NextLowerOrderPolynomial: BernsteinPolynomial + /// a polynomial of the next lower order where each coefficient `b[i]` is defined by `a1 * b[i] + a2 * b[i+1]` + func difference(a1: CGFloat, a2: CGFloat) -> NextLowerOrderPolynomial + /// reduces the polynomial by repeatedly applying `difference` until left with a constant value + func reduce(a1: CGFloat, a2: CGFloat) -> CGFloat + var derivative: NextLowerOrderPolynomial { get } +// init(_ d: Difference, last: CGFloat) +// init(first: CGFloat, _ d: Difference) // func reversed() -> Self -// func split(to x: Double) -> Self -// func split(from x: Double) -> Self -// func split(from tMin: Double, to tMax: Double) -> Self +// func split(to x: CGFloat) -> Self +// func split(from x: CGFloat) -> Self +// func split(from tMin: CGFloat, to tMax: CGFloat) -> Self } -internal extension BernsteinPolynomial { - func f(_ x: Double) -> Double { +internal protocol AnalyticalRoots { + func distinctAnalyticalRoots(between start: CGFloat, and end: CGFloat) -> [CGFloat] +} + +extension BernsteinPolynomial0: AnalyticalRoots { + internal func distinctAnalyticalRoots(between start: CGFloat, and end: CGFloat) -> [CGFloat] { + return [] + } +} + +extension BernsteinPolynomial1: AnalyticalRoots { + internal func distinctAnalyticalRoots(between start: CGFloat, and end: CGFloat) -> [CGFloat] { + var result: [CGFloat] = [] + Utils.droots(self.b0, self.b1) { + guard $0 >= start, $0 <= end else { return } + result.append($0) + } + return result + } +} + +extension BernsteinPolynomial2: AnalyticalRoots { + internal func distinctAnalyticalRoots(between start: CGFloat, and end: CGFloat) -> [CGFloat] { + var result: [CGFloat] = [] + Utils.droots(self.b0, self.b1, self.b2) { + guard $0 >= start, $0 <= end else { return } + result.append($0) + } + return result + } +} + +extension BernsteinPolynomial3: AnalyticalRoots { + internal func distinctAnalyticalRoots(between start: CGFloat, and end: CGFloat) -> [CGFloat] { + var result: [CGFloat] = [] + Utils.droots(self.b0, self.b1, self.b2, self.b3) { + guard $0 >= start, $0 <= end else { return } + result.append($0) + } + return result + } +} + +public extension BernsteinPolynomial { + func value(at x: CGFloat) -> CGFloat { let oneMinusX = 1.0 - x return self.reduce(a1: oneMinusX, a2: x) } - var derivative: Difference { - let order = Double(self.order) + var derivative: NextLowerOrderPolynomial { + let order = CGFloat(self.order) return self.difference(a1: -order, a2: order) } - func reduce(a1: Double, a2: Double) -> Double { + func reduce(a1: CGFloat, a2: CGFloat) -> CGFloat { return self.difference(a1: a1, a2: a2).reduce(a1: a1, a2: a2) } -// func split(to x: Double) -> Self { +// func split(to x: CGFloat) -> Self { // let oneMinusX = 1.0 - x // let difference = self.difference(a1: oneMinusX, a2: x) // let differenceSplit: Difference = difference.split(to: x) // return Self(first: self.first, differenceSplit) // } -// func split(from x: Double) -> Self { +// func split(from x: CGFloat) -> Self { // let oneMinusX = 1.0 - x // let difference = self.difference(a1: oneMinusX, a2: x) // let differenceSplit: Difference = difference.split(from: x) // return Self(differenceSplit, last: self.last) // } -// func split(from tMin: Double, to tMax: Double) -> Self { +// func split(from tMin: CGFloat, to tMax: CGFloat) -> Self { // guard tMax > tMin else { +// #warning("I think this goes into infinite recursion if tMax = tMin = 0.5") // return self.reversed().split(from: 1.0 - tMin, to: 1.0 - tMax) // } // var clippedPolynomial = self.split(to: tMax) @@ -72,181 +118,181 @@ internal extension BernsteinPolynomial { // } } -internal struct BernsteinPolynomial0: BernsteinPolynomial { -// func enumerated(block: (Int, Double) -> Void) { +public struct BernsteinPolynomial0: BernsteinPolynomial { +// func enumerated(block: (Int, CGFloat) -> Void) { // block(0, b0) // } -// var last: Double { return b0 } -// var first: Double { return b0 } -// init(_ d: BernsteinPolynomial0, last: Double) { self.b0 = last } -// init(first: Double, _ d: BernsteinPolynomial0) { self.b0 = first } +// var last: CGFloat { return b0 } +// var first: CGFloat { return b0 } +// init(_ d: BernsteinPolynomial0, last: CGFloat) { self.b0 = last } +// init(first: CGFloat, _ d: BernsteinPolynomial0) { self.b0 = first } // func reversed() -> BernsteinPolynomial0 { return self } -// func split(to x: Double) -> Self { return self } -// func split(from x: Double) -> Self { return self } - init(b0: Double) { self.b0 = b0 } - var b0: Double - var coefficients: [Double] { return [b0] } - func f(_ x: Double) -> Double { +// func split(to x: CGFloat) -> Self { return self } +// func split(from x: CGFloat) -> Self { return self } + public init(b0: CGFloat) { self.b0 = b0 } + public var b0: CGFloat + public var coefficients: [CGFloat] { return [b0] } + public func value(at x: CGFloat) -> CGFloat { return b0 } - var order: Int { return 0 } - func reduce(a1: Double, a2: Double) -> Double { return 0.0 } - func difference(a1: Double, a2: Double) -> BernsteinPolynomial0 { + public var order: Int { return 0 } + public func reduce(a1: CGFloat, a2: CGFloat) -> CGFloat { return 0.0 } + public func difference(a1: CGFloat, a2: CGFloat) -> BernsteinPolynomial0 { return BernsteinPolynomial0(b0: 0.0) } } -internal struct BernsteinPolynomial1: BernsteinPolynomial { -// func enumerated(block: (Int, Double) -> Void) { +public struct BernsteinPolynomial1: BernsteinPolynomial { +// func enumerated(block: (Int, CGFloat) -> Void) { // block(0, b0) // block(1, b1) // } // -// var last: Double { return b1 } -// var first: Double { return b0 } +// var last: CGFloat { return b1 } +// var first: CGFloat { return b0 } // -// init(_ d: BernsteinPolynomial0, last: Double) { +// init(_ d: BernsteinPolynomial0, last: CGFloat) { // self.b0 = d.b0 // self.b1 = last // } // -// init(first: Double, _ d: BernsteinPolynomial0) { +// init(first: CGFloat, _ d: BernsteinPolynomial0) { // self.b0 = first // self.b1 = d.b0 // } // func reversed() -> BernsteinPolynomial1 { BernsteinPolynomial1(b0: b1, b1: b0) } - init(b0: Double, b1: Double) { + public init(b0: CGFloat, b1: CGFloat) { self.b0 = b0 self.b1 = b1 } - typealias Difference = BernsteinPolynomial0 - var b0, b1: Double - var coefficients: [Double] { return [b0, b1] } - func reduce(a1: Double, a2: Double) -> Double { + public typealias NextLowerOrderPolynomial = BernsteinPolynomial0 + public var b0, b1: CGFloat + public var coefficients: [CGFloat] { return [b0, b1] } + public func reduce(a1: CGFloat, a2: CGFloat) -> CGFloat { return a1 * b0 + a2 * b1 } - func difference(a1: Double, a2: Double) -> BernsteinPolynomial0 { + public func difference(a1: CGFloat, a2: CGFloat) -> BernsteinPolynomial0 { return BernsteinPolynomial0(b0: self.reduce(a1: a1, a2: a2)) } - var order: Int { return 1 } + public var order: Int { return 1 } } -internal struct BernsteinPolynomial2: BernsteinPolynomial { -// func enumerated(block: (Int, Double) -> Void) { +public struct BernsteinPolynomial2: BernsteinPolynomial { +// func enumerated(block: (Int, CGFloat) -> Void) { // block(0, b0) // block(1, b1) // block(2, b2) // } -// var last: Double { return b2 } -// var first: Double { return b0 } -// init(_ d: BernsteinPolynomial1, last: Double) { +// var last: CGFloat { return b2 } +// var first: CGFloat { return b0 } +// init(_ d: BernsteinPolynomial1, last: CGFloat) { // self.b0 = d.b0 // self.b1 = d.b1 // self.b2 = last // } -// init(first: Double, _ d: BernsteinPolynomial1) { +// init(first: CGFloat, _ d: BernsteinPolynomial1) { // self.b0 = first // self.b1 = d.b0 // self.b2 = d.b1 // } - init(b0: Double, b1: Double, b2: Double) { + public init(b0: CGFloat, b1: CGFloat, b2: CGFloat) { self.b0 = b0 self.b1 = b1 self.b2 = b2 } - typealias Difference = BernsteinPolynomial1 - var b0, b1, b2: Double - var coefficients: [Double] { return [b0, b1, b2] } - func difference(a1: Double, a2: Double) -> BernsteinPolynomial1 { + public typealias NextLowerOrderPolynomial = BernsteinPolynomial1 + public var b0, b1, b2: CGFloat + public var coefficients: [CGFloat] { return [b0, b1, b2] } + public func difference(a1: CGFloat, a2: CGFloat) -> BernsteinPolynomial1 { return BernsteinPolynomial1(b0: a1 * b0 + a2 * b1, b1: a1 * b1 + a2 * b2) } - var order: Int { return 2 } + public var order: Int { return 2 } } -internal struct BernsteinPolynomial3: BernsteinPolynomial { -// func enumerated(block: (Int, Double) -> Void) { +public struct BernsteinPolynomial3: BernsteinPolynomial { +// func enumerated(block: (Int, CGFloat) -> Void) { // block(0, b0) // block(1, b1) // block(2, b2) // block(3, b3) // } -// var last: Double { return b3 } -// var first: Double { return b0 } -// init(_ d: BernsteinPolynomial2, last: Double) { +// var last: CGFloat { return b3 } +// var first: CGFloat { return b0 } +// init(_ d: BernsteinPolynomial2, last: CGFloat) { // self.b0 = d.b0 // self.b1 = d.b1 // self.b2 = d.b2 // self.b3 = last // } -// init(first: Double, _ d: BernsteinPolynomial2) { +// init(first: CGFloat, _ d: BernsteinPolynomial2) { // self.b0 = first // self.b1 = d.b0 // self.b2 = d.b1 // self.b3 = d.b2 // } - init(b0: Double, b1: Double, b2: Double, b3: Double) { + public init(b0: CGFloat, b1: CGFloat, b2: CGFloat, b3: CGFloat) { self.b0 = b0 self.b1 = b1 self.b2 = b2 self.b3 = b3 } - typealias Difference = BernsteinPolynomial2 - var b0, b1, b2, b3: Double - var coefficients: [Double] { return [b0, b1, b2, b3] } - func difference(a1: Double, a2: Double) -> BernsteinPolynomial2 { + public typealias NextLowerOrderPolynomial = BernsteinPolynomial2 + public var b0, b1, b2, b3: CGFloat + public var coefficients: [CGFloat] { return [b0, b1, b2, b3] } + public func difference(a1: CGFloat, a2: CGFloat) -> BernsteinPolynomial2 { return BernsteinPolynomial2(b0: a1 * b0 + a2 * b1, b1: a1 * b1 + a2 * b2, b2: a1 * b2 + a2 * b3) } - var order: Int { return 3 } + public var order: Int { return 3 } } -internal struct BernsteinPolynomial4: BernsteinPolynomial { -// func enumerated(block: (Int, Double) -> Void) { +public struct BernsteinPolynomial4: BernsteinPolynomial { +// func enumerated(block: (Int, CGFloat) -> Void) { // block(0, b0) // block(1, b1) // block(2, b2) // block(3, b3) // block(4, b4) // } -// var last: Double { return b4 } -// var first: Double { return b0 } -// init(_ d: BernsteinPolynomial3, last: Double) { +// var last: CGFloat { return b4 } +// var first: CGFloat { return b0 } +// init(_ d: BernsteinPolynomial3, last: CGFloat) { // self.b0 = d.b0 // self.b1 = d.b1 // self.b2 = d.b2 // self.b3 = d.b3 // self.b4 = last // } -// init(first: Double, _ d: BernsteinPolynomial3) { +// init(first: CGFloat, _ d: BernsteinPolynomial3) { // self.b0 = first // self.b1 = d.b0 // self.b2 = d.b1 // self.b3 = d.b2 // self.b4 = d.b3 // } - init(b0: Double, b1: Double, b2: Double, b3: Double, b4: Double) { + public init(b0: CGFloat, b1: CGFloat, b2: CGFloat, b3: CGFloat, b4: CGFloat) { self.b0 = b0 self.b1 = b1 self.b2 = b2 self.b3 = b3 self.b4 = b4 } - typealias Difference = BernsteinPolynomial3 - var b0, b1, b2, b3, b4: Double - var coefficients: [Double] { return [b0, b1, b2, b3, b4] } - func difference(a1: Double, a2: Double) -> BernsteinPolynomial3 { + public typealias NextLowerOrderPolynomial = BernsteinPolynomial3 + public var b0, b1, b2, b3, b4: CGFloat + public var coefficients: [CGFloat] { return [b0, b1, b2, b3, b4] } + public func difference(a1: CGFloat, a2: CGFloat) -> BernsteinPolynomial3 { return BernsteinPolynomial3(b0: a1 * b0 + a2 * b1, b1: a1 * b1 + a2 * b2, b2: a1 * b2 + a2 * b3, b3: a1 * b3 + a2 * b4) } - var order: Int { return 4 } + public var order: Int { return 4 } } -internal struct BernsteinPolynomial5: BernsteinPolynomial { -// func enumerated(block: (Int, Double) -> Void) { +public struct BernsteinPolynomial5: BernsteinPolynomial { +// func enumerated(block: (Int, CGFloat) -> Void) { // block(0, b0) // block(1, b1) // block(2, b2) @@ -254,9 +300,9 @@ internal struct BernsteinPolynomial5: BernsteinPolynomial { // block(4, b4) // block(5, b5) // } -// var last: Double { return b5 } -// var first: Double { return b0 } -// init(_ d: BernsteinPolynomial4, last: Double) { +// var last: CGFloat { return b5 } +// var first: CGFloat { return b0 } +// init(_ d: BernsteinPolynomial4, last: CGFloat) { // self.b0 = d.b0 // self.b1 = d.b1 // self.b2 = d.b2 @@ -264,7 +310,7 @@ internal struct BernsteinPolynomial5: BernsteinPolynomial { // self.b4 = d.b4 // self.b5 = last // } -// init(first: Double, _ d: BernsteinPolynomial4) { +// init(first: CGFloat, _ d: BernsteinPolynomial4) { // self.b0 = first // self.b1 = d.b0 // self.b2 = d.b1 @@ -272,7 +318,7 @@ internal struct BernsteinPolynomial5: BernsteinPolynomial { // self.b4 = d.b3 // self.b5 = d.b4 // } - init(b0: Double, b1: Double, b2: Double, b3: Double, b4: Double, b5: Double) { + public init(b0: CGFloat, b1: CGFloat, b2: CGFloat, b3: CGFloat, b4: CGFloat, b5: CGFloat) { self.b0 = b0 self.b1 = b1 self.b2 = b2 @@ -280,39 +326,26 @@ internal struct BernsteinPolynomial5: BernsteinPolynomial { self.b4 = b4 self.b5 = b5 } - typealias Difference = BernsteinPolynomial4 - var b0, b1, b2, b3, b4, b5: Double - var coefficients: [Double] { return [b0, b1, b2, b3, b4, b5] } - func difference(a1: Double, a2: Double) -> BernsteinPolynomial4 { + public typealias NextLowerOrderPolynomial = BernsteinPolynomial4 + public var b0, b1, b2, b3, b4, b5: CGFloat + public var coefficients: [CGFloat] { return [b0, b1, b2, b3, b4, b5] } + public func difference(a1: CGFloat, a2: CGFloat) -> BernsteinPolynomial4 { return BernsteinPolynomial4(b0: a1 * b0 + a2 * b1, b1: a1 * b1 + a2 * b2, b2: a1 * b2 + a2 * b3, b3: a1 * b3 + a2 * b4, b4: a1 * b4 + a2 * b5) } - var order: Int { return 5 } + public var order: Int { return 5 } } -internal extension BernsteinPolynomial { - func analyticalRoots(between start: Double, and end: Double) -> [Double]? { - let order = self.order - guard order > 0 else { return [] } - guard order < 4 else { return nil } // cannot solve - return Utils.droots(self.coefficients.map { CGFloat($0) }).compactMap { - let t = Double($0) - guard t > start, t < end else { return nil } - return t - } - } -} - -private func newton(polynomial: P, derivative: P.Difference, guess: Double, relaxation: Double = 1) -> Double { +private func newton(polynomial: P, derivative: P.NextLowerOrderPolynomial, guess: CGFloat, relaxation: CGFloat = 1) -> CGFloat { let maxIterations = 20 var x = guess for _ in 0..(polynomial: P, derivative: P.Differe return x } -private func findRootBisection(of polynomial: P, start: Double, end: Double) -> Double { +private func findRootBisection(of polynomial: P, start: CGFloat, end: CGFloat) -> CGFloat { var guess = (start + end) / 2 var low = start var high = end - let lowSign = polynomial.f(low).sign - let highSign = polynomial.f(high).sign + let lowSign = polynomial.value(at: low).sign + let highSign = polynomial.value(at: high).sign assert(lowSign != highSign) let maxIterations = 20 var iterations = 0 while high - low > 1.0e-5 { let midGuess = (low + high) / 2 guess = midGuess - let nextGuessF = polynomial.f(guess) + let nextGuessF = polynomial.value(at: guess) if nextGuessF == 0 { return guess } else if nextGuessF.sign == lowSign { @@ -348,21 +381,25 @@ private func findRootBisection(of polynomial: P, start: return guess } -internal func findRoots(of polynomial: P, between start: Double, and end: Double) -> [Double] { +public func findDistinctRootsInUnitInterval(of polynomial: P) -> [CGFloat] { + return findDistinctRoots(of: polynomial, between: 0, and: 1) +} + +internal func findDistinctRoots(of polynomial: P, between start: CGFloat, and end: CGFloat) -> [CGFloat] { assert(start < end) - if let roots = polynomial.analyticalRoots(between: start, and: end) { - return roots + if let analytical = polynomial as? AnalyticalRoots { + return analytical.distinctAnalyticalRoots(between: start, and: end) } let derivative = polynomial.derivative - let criticalPoints: [Double] = findRoots(of: derivative, between: start, and: end) - let intervals: [Double] = [start] + criticalPoints + [end] - var lastFoundRoot: Double? - let roots = (0.. Double? in + let criticalPoints: [CGFloat] = findDistinctRoots(of: derivative, between: start, and: end) + let intervals: [CGFloat] = [start] + criticalPoints + [end] + var lastFoundRoot: CGFloat? + let roots = (0.. CGFloat? in let start = intervals[i] let end = intervals[i+1] - let fStart = polynomial.f(start) - let fEnd = polynomial.f(end) - let root: Double + let fStart = polynomial.value(at: start) + let fEnd = polynomial.value(at: end) + let root: CGFloat if fStart * fEnd < 0 { // TODO: if a critical point is a root we take this // codepath due to roundoff and converge only linearly to one end of interval @@ -382,7 +419,7 @@ internal func findRoots(of polynomial: P, between start: guard abs(value - guess) < 1.0e-5 else { return nil // did not converge near guess } - guard abs(polynomial.f(value)) < 1.0e-10 else { + guard abs(polynomial.value(at: value)) < 1.0e-10 else { return nil // not actually a root } root = value @@ -398,15 +435,15 @@ internal func findRoots(of polynomial: P, between start: return roots } -//internal func findRoots(of polynomial: P, between start: Double, and end: Double) -> [Double] { +//internal func findRoots(of polynomial: P, between start: CGFloat, and end: CGFloat) -> [CGFloat] { // assert(start < end) // -// var tMin: Double = Double.infinity -// var tMax: Double = -Double.infinity +// var tMin: CGFloat = CGFloat.infinity +// var tMax: CGFloat = -CGFloat.infinity // var intersected = false // -// func x(_ i: Int) -> Double { -// return Double(i) / Double(polynomial.order) +// func x(_ i: Int) -> CGFloat { +// return CGFloat(i) / CGFloat(polynomial.order) // } // // compute the intersections of each pair of lines with the x axis // polynomial.enumerated { i, c1 in @@ -439,7 +476,7 @@ internal func findRoots(of polynomial: P, between start: // assert(tMax >= tMin) // // // find [adjustedStart, adjustedEnd] range represented by [tMin, tMax] in original polynomial -// func adjustedT(_ t: Double) -> Double { +// func adjustedT(_ t: CGFloat) -> CGFloat { // return start * (1.0 - t) + end * t // } // let adjustedStart = adjustedT(tMin) diff --git a/BezierKit/Library/Utils.swift b/BezierKit/Library/Utils.swift index df9b1cf..97d891b 100644 --- a/BezierKit/Library/Utils.swift +++ b/BezierKit/Library/Utils.swift @@ -262,23 +262,6 @@ internal class Utils { callback(p0 / (p0 - p1)) } - static func droots(_ p: [CGFloat]) -> [CGFloat] { - // quadratic roots are easy - var result: [CGFloat] = [] - let callback = { result.append($0) } - switch p.count { - case 4: - droots(p[0], p[1], p[2], p[3], callback: callback) - case 3: - droots(p[0], p[1], p[2], callback: callback) - case 2: - droots(p[0], p[1], callback: callback) - default: - fatalError("unsupported") - } - return result - } - static func lerp(_ r: CGFloat, _ v1: CGPoint, _ v2: CGPoint) -> CGPoint { return v1 + r * (v2 - v1) }