Bezier and b-spline arc-length algorithm giving me problems
I'm having a bit of a problem calculating the arc-length of my bezier and b-spline curves. I've been banging my head against this for several days, and I think I'm almost there, but can't seem to get it exactly right. I'm developing in Swift, but I think its syntax is clear enough that anyone who knows C/C++ would be able to read it. If not, please let me know and I'll try to translate it into C/C++.
I've checked my implementations against several sources over and over again, and, as far as the algorithms go, they seem to be correct, although I'm not so sure about the B-spline algorithm. Some tutorials use the degree, and some use the order, of the curve in their calculations, and I get really confused. In addition, in using the Gauss-Legendre quadrature, I understand that I'm supposed to sum the integration of the spans, but I'm not sure I'm understanding how to do that correctly. From what I understand, I should be integrating over each knot span. Is that correct?
When I calculate the length of a Bezier curve with the following control polygon, I get 28.2842712474619, while 3D software (Cinema 4D and Maya) tells me the length should be 30.871.
let bezierControlPoints = [
Vector(-10.0, -10.0),
Vector(0.0, -10.0),
Vector(0.0, 10.0),
Vector(10.0, 10.0)
]
The length of the b-spline is similarly off. My algorithm produces 5.6062782185353, while it should be 7.437.
let splineControlPoints = [
Vector(-2.0, -1.0),
Vector(-1.0, 1.0),
Vector(-0.25, 1.0),
Vector(0.25, -1.0),
Vector(1.0, -1.0),
Vector(2.0, 1.0)
]
I'm not a mathematician, so I'm struggling with the math, but I think I have the gist of it.
The Vector class is pretty straight-forwared, but I've overloaded some operators for convenience/legibility which makes the code quite lengthy, so I'm not posting it here. I'm also not including the Gauss-Legendre weights and abscissae. You can download the source and Xcode project from here (53K).
Here's my bezier curve class:
class Bezier
{
var c0:Vector
var c1:Vector
var c2:Vector
var c3:Vector
init(ic0 _ic0:Vector, ic1 _ic1:Vector, ic2 _ic2:Vector, ic3 _ic3:Vector) {
c0 = _ic0
c1 = _ic1
c2 = _ic2
c3 = _ic3
}
// Calculate curve length using Gauss-Legendre quadrature
func curveLength()->Double {
let gl = GaussLegendre()
gl.order = 3 // Good enough for a quadratic polynomial
let xprime = gl.integrate(a:0.0, b:1.0, closure:{ (t:Double)->Double in return self.dx(atTime:t) })
let yprime = gl.integrate(a:0.0, b:1.0, closure:{ (t:Double)->Double in return self.dy(atTime:t) })
return sqrt(xprime*xprime + yprime*yprime)
}
// I could vectorize this, but correctness > efficiency
// The derivative of the x-component
func dx(atTime t:Double)->Double {
let tc = (1.0-t)
let r0 = (3.0 * tc*tc) * (c1.x - c0.x)
let r1 = (6.0 * tc*t) * (c2.x - c1.x)
let r2 = (3.0 * t*t) * (c3.x - c2.x)
return r0 + r1 + r2
}
// The derivative of the y-component
func dy(atTime t:Double)->Double {
let tc = (1.0-t)
let r0 = (3.0 * tc*tc) * (c1.y - c0.y)
let r1 = (6.0 * tc*t) * (c2.y - c1.y)
let r2 = (3.0 * t*t) * (c3.y - c2.y)
return r0 + r1 + r2
}
}
Here is my b-spline class:
class BSpline
{
var spanLengths:[Double]! = nil
var totalLength:Double = 0.0
var cp:[Vector]
var knots:[Double]! = nil
var o:Int = 4
init(controlPoints:[Vector]) {
cp = controlPoints
calcKnots()
}
// Method to return length of the curve using Gauss-Legendre numerical integration
func cacheSpanLengths() {
spanLengths = [Double]()
totalLength = 0.0
let gl = GaussLegendre()
gl.order = o-1 // The derivative should be quadratic, so o-2 would suffice?
// Am I doing this right? Piece-wise integration?
for i in o-1 ..< knots.count-o {
let t0 = knots[i]
let t1 = knots[i+1]
let xprime = gl.integrate(a:t0, b:t1, closure:self.dx)
let yprime = gl.integrate(a:t0, b:t1, closure:self.dy)
let spanLength = sqrt(xprime*xprime + yprime*yprime)
spanLengths.append(spanLength)
totalLength += spanLength
}
}
// The b-spline basis function
func basis(i:Int, _ k:Int, _ x:Double)->Double {
var r:Double = 0.0
switch k {
case 0:
if (knots[i] <= x) && (x <= knots[i+1]) {
r = 1.0
} else {
r = 0.0
}
default:
var n0 = x - knots[i]
var d0 = knots[i+k]-knots[i]
var b0 = basis(i,k-1,x)
var n1 = knots[i+k+1] - x
var d1 = knots[i+k+1]-knots[i+1]
var b1 = basis(i+1,k-1,x)
var left = Double(0.0)
var right = Double(0.0)
if b0 != 0 && d0 != 0 { left = n0 * b0 / d0 }
if b1 != 0 && d1 != 0 { right = n1 * b1 / d1 }
r = left + right
}
return r
}
// Method to calculate and store the knot vector
func calcKnots() {
// The number of knots in the knot vector = number of control points + order (i.e. degree + 1)
let knotCount = cp.count + o
knots = [Double]()
// For an open b-spline where the ends are incident on the first and last control points,
// the first o knots are the same and the last o knots are the same, where o is the order
// of the curve.
var k = 0
for i in 0 ..< o {
knots.append(0.0)
}
for i in o ..< cp.count {
k++
knots.append(Double(k))
}
k++
for i in cp.count ..< knotCount {
knots.append(Double(k))
}
}
// I could vectorize this, but correctness > efficiency
// Derivative of the x-component
func dx(t:Double)->Double {
var p = Double(0.0)
let n = o
for i in 0 ..< cp.count-1 {
let u0 = knots[i + n + 1]
let u1 = knots[i + 1]
let fn = Double(n) / (u0 - u1)
let thePoint = (cp[i+1].x - cp[i].x) * fn
let b = basis(i+1, n-1, Double(t))
p += thePoint * b
}
return Double(p)
}
// Derivative of the y-component
func dy(t:Double)->Double {
var p = Double(0.0)
let n = o
for i in 0 ..< cp.count-1 {
let u0 = knots[i + n + 1]
let u1 = knots[i + 1]
let fn = Double(n) / (u0 - u1)
let thePoint = (cp[i+1].y - cp[i].y) * fn
let b = basis(i+1, n-1, Double(t))
p += thePoint * b
}
return Double(p)
}
}
And here is my Gauss-Legendre implementation:
class GaussLegendre
{
var order:Int = 5
init() {
}
// Numerical integration of arbitrary function
func integrate(a _a:Double, b _b:Double, closure f:(Double)->Double)->Double {
var result = 0.0
let wgts = gl_weights[order-2]
let absc = gl_abscissae[order-2]
for i in 0..<order {
let a0 = absc[i]
let w0 = wgts[i]
result += w0 * f(0.5 * (_b + _a + a0 * (_b - _a)))
}
return 0.5 * (_b - _a) * result
}
}
And my main logic:
let bezierControlPoints = [
Vector(-10.0, -10.0),
Vector(0.0, -10.0),
Vector(0.0, 10.0),
Vector(10.0, 10.0)
]
let splineControlPoints = [
Vector(-2.0, -1.0),
Vector(-1.0, 1.0),
Vector(-0.25, 1.0),
Vector(0.25, -1.0),
Vector(1.0, -1.0),
Vector(2.0, 1.0)
]
var bezier = Bezier(controlPoints:bezierControlPoints)
println("Bezier curve length: \(bezier.curveLength())\n")
var spline:BSpline = BSpline(controlPoints:splineControlPoints)
spline.cacheSpanLengths()
println("B-Spline curve length: \(spline.totalLength)\n")
UPDATE: PROBLEM (PARTIALLY) SOLVED
Thanks to Mike for his answer!
I verified that I am correctly remapping the numerical integration from the interval a..b to -1..1 for the purposes of Legendre-Gauss quadrature. The math is here (apologies to any real mathematicians out there, it's the best I could do with my long-forgotten calculus).
I've increased the order of the Legendre-Gauss quadrature from 5 to 32 as Mike suggested.
Then after a lot of floundering around in Mathematica, I came back and re-read Mike's code and discovered that my code was NOT equivalent to his.
I was taking the square root of the sums of the squared integrals of the derivative components:
when I should have been taking the integral of the magnitudes of the derivative vectors:
In terms of code, in my Bezier class, instead of this:
// INCORRECT
func curveLength()->Double {
let gl = GaussLegendre()
gl.order = 3 // Good enough for a quadratic polynomial
let xprime = gl.integrate(a:0.0, b:1.0, closure:{ (t:Double)->Double in return self.dx(atTime:t) })
let yprime = gl.integrate(a:0.0, b:1.0, closure:{ (t:Double)->Double in return self.dy(atTime:t) })
return sqrt(xprime*xprime + yprime*yprime)
}
I should have written this:
// CORRECT
func curveLength()->Double {
let gl = GaussLegendre()
gl.order = 32
return = gl.integrate(a:0.0, b:1.0, closure:{ (t:Double)->Double in
let x = self.dx(atTime:t)
let y = self.dy(atTime:t)
return sqrt(x*x + y*y)
})
}
My code calculates the arc length as: 3.59835872777095
Mathematica: 3.598358727834686
So, my result is pretty close. Interestingly, there is a discrepancy between a plot in Mathematica of my test Bezier curve, and the same rendered by Cinema 4D, which would explain why the arc lengths calculated by Mathematica and Cinema 4D are different as well. I think I trust Mathematica to be more correct, though.
In my B-Spline class, instead of this:
// INCORRECT
func cacheSpanLengths() {
spanLengths = [Double]()
totalLength = 0.0
let gl = GaussLegendre()
gl.order = o-1 // The derivative should be quadratic, so o-2 would suffice?
// Am I doing this right? Piece-wise integration?
for i in o-1 ..< knots.count-o {
let t0 = knots[i]
let t1 = knots[i+1]
let xprime = gl.integrate(a:t0, b:t1, closure:self.dx)
let yprime = gl.integrate(a:t0, b:t1, closure:self.dy)
let spanLength = sqrt(xprime*xprime + yprime*yprime)
spanLengths.append(spanLength)
totalLength += spanLength
}
}
I should have written this:
// CORRECT
func cacheSpanLengths() {
spanLengths = [Double]()
totalLength = 0.0
let gl = GaussLegendre()
gl.order = 32
// Am I doing this right? Piece-wise integration?
for i in o-1 ..< knots.count-o {
let t0 = knots[i]
let t1 = knots[i+1]
let spanLength = gl.integrate(a:t0, b:t1, closure:{ (t:Double)->Double in
let x = self.dx(atTime:t)
let y = self.dy(atTime:t)
return sqrt(x*x + y*y)
})
spanLengths.append(spanLength)
totalLength += spanLength
}
}
Unfortunately, the B-Spline math is not as straight-forward, and I haven't been able to test it in Mathematica as easily as the Bezier math, so I'm not entirely sure my code is working, even with the above changes. I will post another update when I verify it.
UPDATE 2: PROBLEM SOLVED
Eureka, I discovered an off-by one error in my code to calculate the B-Spline derivative.
Instead of
// Derivative of the x-component
func dx(t:Double)->Double {
var p = Double(0.0)
let n = o // INCORRECT (should be one less)
for i in 0 ..< cp.count-1 {
let u0 = knots[i + n + 1]
let u1 = knots[i + 1]
let fn = Double(n) / (u0 - u1)
let thePoint = (cp[i+1].x - cp[i].x) * fn
let b = basis(i+1, n-1, Double(t))
p += thePoint * b
}
return Double(p)
}
// Derivative of the y-component
func dy(t:Double)->Double {
var p = Double(0.0)
let n = o // INCORRECT (should be one less_
for i in 0 ..< cp.count-1 {
let u0 = knots[i + n + 1]
let u1 = knots[i + 1]
let fn = Double(n) / (u0 - u1)
let thePoint = (cp[i+1].y - cp[i].y) * fn
let b = basis(i+1, n-1, Double(t))
p += thePoint * b
}
return Double(p)
}
I should have written
// Derivative of the x-component
func dx(t:Double)->Double {
var p = Double(0.0)
let n = o-1 // CORRECT
for i in 0 ..< cp.count-1 {
let u0 = knots[i + n + 1]
let u1 = knots[i + 1]
let fn = Double(n) / (u0 - u1)
let thePoint = (cp[i+1].x - cp[i].x) * fn
let b = basis(i+1, n-1, Double(t))
p += thePoint * b
}
return Double(p)
}
// Derivative of the y-component
func dy(t:Double)->Double {
var p = Double(0.0)
let n = o-1 // CORRECT
for i in 0 ..< cp.count-1 {
let u0 = knots[i + n + 1]
let u1 = knots[i + 1]
let fn = Double(n) / (u0 - u1)
let thePoint = (cp[i+1].y - cp[i].y) * fn
let b = basis(i+1, n-1, Double(t))
p += thePoint * b
}
return Double(p)
}
My code now calculates the length of the B-Spline curve as 6.87309971722132. Mathematica: 6.87309884638438.
It's probably not scientifically precise, but good enough for me.
The Legendre-Gauss procedure is specifically defined for the interval [-1,1], whereas Beziers and B-Splines are defined over [0,1], so that's a simple conversion and at least while you're trying to make sure your code does the right thing, easy to bake in instead of supplying a dynamic interval (as you say, accuracy over efficiency. Once it works, we can worry about optimising)
So, given weights W and abscissae A (both of same length n), you'd do:
z = 0.5
for i in 1..n
w = W[i]
a = A[i]
t = z * a + z
sum += w * arcfn(t, xpoints, ypoints)
return z * sum
with the pseudo-code assuming list indexing from 1. The arcfn would be defined as:
arcfn(t, xpoints, ypoints):
x = derive(xpoints, t)
y = derive(ypoints, t)
c = x*x + y*y
return sqrt(c)
But that part looks right already.
Your derivatives look correct too, so the main question is: "are you using enough slices in your Legendre-Gauss quadrature?". Your code suggests you're using only 5 slices, which isn't nearly enough to get a good result. Using http://pomax.github.io/bezierinfo/legendre-gauss.html as term data, you generally want a set for n of 16 or higher (for cubic Bezier curves, 24 is generally safe, although still underperformant for curves with cusps or lots of inflections).
I can recommend taking the "unit test" approach here: test your bezier and bspline code (separately) for known base and derivative values. Do those check out? One problem ruled out. On to your LG code: if you perform Legendre-Gauss on a parametric function for a straight line using:
fx(t) = t
fy(t) = t
fx'(t) = 1
fy'(t) = 1
over interval t=[0,1], we know the length should be exactly the square root of 2, and the derivatives are the simplest possible. If those work, do a non-linear test using:
fx(t) = sin(t)
fy(t) = cos(t)
fx'(t) = cos(t)
fy'(t) = -sin(t)
over interval t=[0,1]; we know the length should be exactly 1. Does your LG implementation yield the correct value? Another problem ruled out. If it doesn't, check your weights and abscissae. Do they match the ones from the linked page (generated with a verifiably correct Mathematica program, so pretty much guaranteed to be correct)? Are you using enough slices? Bump the number up to 10, 16, 24, 32; increasing the number of slices will show a stabilising summation, where adding more slices doesn't change digits before the 2nd, 3rd, 4th, 5th, etc decimal point as you increase the count.
Are the curves you're testing with known to be problematic curves? Plot them, do they have cusps or lots of inflections? That's going to be a problem for LG, try simpler curves to see if the values you get back for those, at least, are correct.
Finally, check your types: Are you using the highest precision possible datatype? 32 bit floats are going to run into mysteriously disappearing FPU and wonderful rounding errors at the values we need to use when doing LG with a reasonable number of slices.
- iOS Alternate Icon not changing
- Picker still interactable when disabled
- Property wrappers and SwiftUI environment: how can a property wrapper access the environment of its enclosing object?
- Why use @ObservableObject when you might as well use @StateObject?
- Live update relative time in SwiftUI
- Using the Apple Music API to create playlist?
- Error when compiling AsyncStack<> with Swift 6
- Why does my snapshot() in SwiftUI not capture downloaded images from WebImage (SDWebImageSwiftUI)?
- How to avoid sharing state on WindowGroup / multiscreen macOS with SwiftUI?
- How can I change the UISearchBar search text color?
- Swift Combine Pipeline to compare draft model
- Multiple function calls, but single execution of async Task
- Data in HTTPBody with a PUT method fails, while it works with a POST?
- How to use a custom activity indicator view while using API call in my project
- Make sheet the exact size of the content inside
- Send post variables to API
- Setting properties in init method using closures
- How can my "didFinishLaunchingWithOptions" wait until call API finish
- EXC_BAD_INSTRUCTION (code=EXC_i386_INVOP, subcode=0x0) crash
- API authentication with swift
- Swift Best Practices - how / when to check for internet connection and 404s when working with APIs
- Reading Websites in iOS
- I cannot play music in the background (Xcode 9, Swift 4)
- Variable was never mutated...or Constant being used before initialized
- JSON API transfer data in Swift
- Having trouble decoding JSON file from GitHub API swift 5
- Swift 5 Parse multi level JSON
- Decoding JSON error "but found a dictionary instead"
- Swift weather data from Openweathermap API?
- How to add placeholder text to TextEditor in SwiftUI?