swift - Bezier and b-spline arc-length algorithm giving me problems -- SOLVED -
i'm having bit of problem calculating arc-length of bezier , b-spline curves. i've been banging head against several days, , think i'm there, can't seem right. i'm developing in swift, think syntax clear plenty knows c/c++ able read it. if not, please allow me know , i'll seek translate c/c++.
i've checked implementations against several sources on , on again, and, far algorithms go, seem correct, although i'm not sure b-spline algorithm. tutorials utilize degree, , utilize order, of curve in calculations, , confused. in addition, in using gauss-legendre quadrature, understand i'm supposed sum integration of spans, i'm not sure i'm understanding how correctly. understand, should integrating on each knot span. correct?
thanks in advance assistance can give me.
when calculate length of bezier curve next command polygon, 28.2842712474619, while 3d software (cinema 4d , maya) tells me length should 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 b-spline off. algorithm produces 5.6062782185353, while should 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 mathematician, i'm struggling math, think have gist of it.
the vector class pretty straight-forwared, i've overloaded operators convenience/legibility makes code quite lengthy, i'm not posting here. i'm not including gauss-legendre weights , abscissae. can download source , xcode project here (53k).
here's 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 { allow gl = gausslegendre() gl.order = 3 // plenty quadratic polynomial allow xprime = gl.integrate(a:0.0, b:1.0, closure:{ (t:double)->double in homecoming self.dx(attime:t) }) allow yprime = gl.integrate(a:0.0, b:1.0, closure:{ (t:double)->double in homecoming self.dy(attime:t) }) homecoming sqrt(xprime*xprime + yprime*yprime) } // vectorize this, correctness > efficiency // derivative of x-component func dx(attime t:double)->double { allow tc = (1.0-t) allow r0 = (3.0 * tc*tc) * (c1.x - c0.x) allow r1 = (6.0 * tc*t) * (c2.x - c1.x) allow r2 = (3.0 * t*t) * (c3.x - c2.x) homecoming r0 + r1 + r2 } // derivative of y-component func dy(attime t:double)->double { allow tc = (1.0-t) allow r0 = (3.0 * tc*tc) * (c1.y - c0.y) allow r1 = (6.0 * tc*t) * (c2.y - c1.y) allow r2 = (3.0 * t*t) * (c3.y - c2.y) homecoming r0 + r1 + r2 } } here 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 homecoming length of curve using gauss-legendre numerical integration func cachespanlengths() { spanlengths = [double]() totallength = 0.0 allow gl = gausslegendre() gl.order = o-1 // derivative should quadratic, o-2 suffice? // doing right? piece-wise integration? in o-1 ..< knots.count-o { allow t0 = knots[i] allow t1 = knots[i+1] allow xprime = gl.integrate(a:t0, b:t1, closure:self.dx) allow yprime = gl.integrate(a:t0, b:t1, closure:self.dy) allow spanlength = sqrt(xprime*xprime + yprime*yprime) spanlengths.append(spanlength) totallength += spanlength } } // 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 } homecoming r } // method calculate , store knot vector func calcknots() { // number of knots in knot vector = number of command points + order (i.e. grade + 1) allow knotcount = cp.count + o knots = [double]() // open b-spline ends incident on first , lastly command points, // first o knots same , lastly o knots same, o order // of curve. var k = 0 in 0 ..< o { knots.append(0.0) } in o ..< cp.count { k++ knots.append(double(k)) } k++ in cp.count ..< knotcount { knots.append(double(k)) } } // vectorize this, correctness > efficiency // derivative of x-component func dx(t:double)->double { var p = double(0.0) allow n = o in 0 ..< cp.count-1 { allow u0 = knots[i + n + 1] allow u1 = knots[i + 1] allow fn = double(n) / (u0 - u1) allow thepoint = (cp[i+1].x - cp[i].x) * fn allow b = basis(i+1, n-1, double(t)) p += thepoint * b } homecoming double(p) } // derivative of y-component func dy(t:double)->double { var p = double(0.0) allow n = o in 0 ..< cp.count-1 { allow u0 = knots[i + n + 1] allow u1 = knots[i + 1] allow fn = double(n) / (u0 - u1) allow thepoint = (cp[i+1].y - cp[i].y) * fn allow b = basis(i+1, n-1, double(t)) p += thepoint * b } homecoming double(p) } } and here 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 allow wgts = gl_weights[order-2] allow absc = gl_abscissae[order-2] in 0..<order { allow a0 = absc[i] allow w0 = wgts[i] result += w0 * f(0.5 * (_b + _a + a0 * (_b - _a))) } homecoming 0.5 * (_b - _a) * result } } and main logic:
let beziercontrolpoints = [ vector(-10.0, -10.0), vector(0.0, -10.0), vector(0.0, 10.0), vector(10.0, 10.0) ] allow 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 mike answer!
i verified correctly remapping numerical integration interval a..b -1..1 purposes of legendre-gauss quadrature. math here (apologies real mathematicians out there, it's best long-forgotten calculus).
i've increased order of legendre-gauss quadrature 5 32 mike suggested.
then after lot of floundering around in mathematica, came , re-read mike's code , discovered code not equivalent his.
i taking square root of sums of squared integrals of derivative components:
when should have been taking integral of magnitudes of derivative vectors:
in terms of code, in bezier class, instead of this:
// wrong func curvelength()->double { allow gl = gausslegendre() gl.order = 3 // plenty quadratic polynomial allow xprime = gl.integrate(a:0.0, b:1.0, closure:{ (t:double)->double in homecoming self.dx(attime:t) }) allow yprime = gl.integrate(a:0.0, b:1.0, closure:{ (t:double)->double in homecoming self.dy(attime:t) }) homecoming sqrt(xprime*xprime + yprime*yprime) } i should have written this:
// right func curvelength()->double { allow gl = gausslegendre() gl.order = 32 homecoming = gl.integrate(a:0.0, b:1.0, closure:{ (t:double)->double in allow x = self.dx(attime:t) allow y = self.dy(attime:t) homecoming sqrt(x*x + y*y) }) } my code calculates arc length as: 3.59835872777095 mathematica: 3.598358727834686
so, result pretty close. interestingly, there discrepancy between plot in mathematica of test bezier curve, , same rendered cinema 4d, explain why arc lengths calculated mathematica , cinema 4d different well. think trust mathematica more correct, though.
in b-spline class, instead of this:
// wrong func cachespanlengths() { spanlengths = [double]() totallength = 0.0 allow gl = gausslegendre() gl.order = o-1 // derivative should quadratic, o-2 suffice? // doing right? piece-wise integration? in o-1 ..< knots.count-o { allow t0 = knots[i] allow t1 = knots[i+1] allow xprime = gl.integrate(a:t0, b:t1, closure:self.dx) allow yprime = gl.integrate(a:t0, b:t1, closure:self.dy) allow spanlength = sqrt(xprime*xprime + yprime*yprime) spanlengths.append(spanlength) totallength += spanlength } } i should have written this:
// right func cachespanlengths() { spanlengths = [double]() totallength = 0.0 allow gl = gausslegendre() gl.order = 32 // doing right? piece-wise integration? in o-1 ..< knots.count-o { allow t0 = knots[i] allow t1 = knots[i+1] allow spanlength = gl.integrate(a:t0, b:t1, closure:{ (t:double)->double in allow x = self.dx(attime:t) allow y = self.dy(attime:t) homecoming sqrt(x*x + y*y) }) spanlengths.append(spanlength) totallength += spanlength } } unfortunately, b-spline math not straight-forward, , haven't been able test in mathematica bezier math, i'm not exclusively sure code working, above changes. post update when verify it.
update 2: problem solved
eureka, discovered off-by 1 error in code calculate b-spline derivative.
instead of
// derivative of x-component func dx(t:double)->double { var p = double(0.0) allow n = o // wrong (should 1 less) in 0 ..< cp.count-1 { allow u0 = knots[i + n + 1] allow u1 = knots[i + 1] allow fn = double(n) / (u0 - u1) allow thepoint = (cp[i+1].x - cp[i].x) * fn allow b = basis(i+1, n-1, double(t)) p += thepoint * b } homecoming double(p) } // derivative of y-component func dy(t:double)->double { var p = double(0.0) allow n = o // wrong (should 1 less_ in 0 ..< cp.count-1 { allow u0 = knots[i + n + 1] allow u1 = knots[i + 1] allow fn = double(n) / (u0 - u1) allow thepoint = (cp[i+1].y - cp[i].y) * fn allow b = basis(i+1, n-1, double(t)) p += thepoint * b } homecoming double(p) } i should have written
// derivative of x-component func dx(t:double)->double { var p = double(0.0) allow n = o-1 // right in 0 ..< cp.count-1 { allow u0 = knots[i + n + 1] allow u1 = knots[i + 1] allow fn = double(n) / (u0 - u1) allow thepoint = (cp[i+1].x - cp[i].x) * fn allow b = basis(i+1, n-1, double(t)) p += thepoint * b } homecoming double(p) } // derivative of y-component func dy(t:double)->double { var p = double(0.0) allow n = o-1 // right in 0 ..< cp.count-1 { allow u0 = knots[i + n + 1] allow u1 = knots[i + 1] allow fn = double(n) / (u0 - u1) allow thepoint = (cp[i+1].y - cp[i].y) * fn allow b = basis(i+1, n-1, double(t)) p += thepoint * b } homecoming double(p) } my code calculates length of b-spline curve 6.87309971722132. mathematica: 6.87309884638438.
it's not scientifically precise, plenty me.
the legendre-gauss procedure defined interval [-1,1], whereas beziers , b-splines defined on [0,1], that's simple conversion , @ to the lowest degree while you're trying create sure code right thing, easy bake in instead of supplying dynamic interval (as say, accuracy on efficiency. 1 time works, can worry optimising)
so, given weights w , abscissae a (both of same length n), you'd do:
z = 0.5 in 1..n w = w[i] = a[i] t = z * + z sum += w * arcfn(t, xpoints, ypoints) homecoming z * sum with pseudo-code assuming list indexing 1. arcfn defined as:
arcfn(t, xpoints, ypoints): x = derive(xpoints, t) y = derive(ypoints, t) c = x*x + y*y homecoming sqrt(c) but part looks right already.
your derivatives right too, main question is: "are using plenty slices in legendre-gauss quadrature?". code suggests you're using 5 slices, isn't plenty result. using http://pomax.github.io/bezierinfo/legendre-gauss.html term data, want set n of 16 or higher (for cubic bezier curves, 24 safe, although still underperformant curves cusps or lots of inflections).
i can recommend taking "unit test" approach here: test bezier , bspline code (separately) known base of operations , derivative values. check out? 1 problem ruled out. on lg code: if perform legendre-gauss on parametric function straight line using:
fx(t) = t fy(t) = t fx'(t) = 1 fy'(t) = 1 over interval t=[0,1], know length should exactly square root of 2, , derivatives simplest possible. if work, 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]; know length should exactly 1. lg implementation yield right value? problem ruled out. if doesn't, check weights , abscissae. match ones linked page (generated verifiably right mathematica program, pretty much guaranteed correct)? using plenty slices? bump number 10, 16, 24, 32; increasing number of slices show stabilising summation, adding more slices doesn't alter digits before 2nd, 3rd, 4th, 5th, etc decimal point increment count.
are curves you're testing known problematic curves? plot them, have cusps or lots of inflections? that's going problem lg, seek simpler curves see if values those, @ least, correct.
finally, check types: using highest precision possible datatype? 32 bit floats going run mysteriously disappearing fpu , wonderful rounding errors @ values need utilize when doing lg reasonable number of slices.
swift bezier numerical-integration
