🔒 Locked

Calculating the length of a Bezier curve...

Started by iNsAn1tY Apr 12, 2005 at 5:41 PM 7 replies 50.5k views

Original Post

Apr 12, 2005 at 5:41 PM

Evening all. I'm doing a bit of work with 3D cubic Bezier curves, and I've run into a bit of a problem. I need to know (ideally) the exact length of the curve I'm working with. I know that I can get an approximation of the length by creating straight lines along the length of the curve, then adding together their lengths, but I was wondering if there was a simpler solution in the form of an equation that I can solve for the curve. I can't seem to find much on Google about it, so any input is greatly appreciated...

averisk

Apr 12, 2005 at 6:28 PM

I've considered this problem before... solving it analytically appears to be impossible, since it involves an integral that usually doesn't integrate in general

http://www.cs.uiowa.edu/~kearney/22c196Spring03/ArcLengthParameterization.ppt
<-- Here's something you might find useful (a description of a numeric approach)

Charles B

Apr 12, 2005 at 7:20 PM

The generic formula for a parametric curve length is

Sum[0,1]( sqrt((dx/dt)^2+(dy/dt)^2+(dz/dt)^2)*dt )

BTW, to "verify", for a classical function curve y=f(x), just replace t by x, discard z and you get Sum[xa,xb]( sqrt(1+(dy/dx)^2)*dx )

Since your Bezier spline is of degree 3, dx/dt is of degree 2. (dx/dt)^2 is of degree 4, etc ... All in all you end up with the integration of the square root of a degree 4 polynomial in t.

Then there should be enough litterature on the www for you concerning numerical or formal integration techniques. But this can easilly be solved formally :

Go there http://torte.cs.berkeley.edu:8010/tilu and you'll get the formal answer you need. Select Mathematica as language and type (indefinite integration : finding a primitive):

sqrt(a*x^4+b*x^3+c*x^2+d*x+e)

You'll get (control by differentiating the result, it's easy):

             5          4          2          3sqrt (1/5 a x  + 1/4 b x  + 1/2 d x  + 1/3 c x  + e x)

And in fact replace x by 1 since you just integrated between 0 and 1 and this makes 0 at 0.

length = sqrt(1/5a + 1/4b + 1/2d + 1/3c + e)

I hope you understood how to get a, b, c, d and e from your spline.

With a quick SSE or 3DNow sqrt or the "Carmack" sqrt, you can have both precision and speed. Google with Gamedev search engine to know more about quick sqrts.

Charles B

Apr 12, 2005 at 7:28 PM

Quote:
Original post by averisk
I've considered this problem before... solving it analytically appears to be impossible, since it involves an integral that usually doesn't integrate in general

What I thought, till I realized I had already coded that in a commercial game a long time ago. And after putting the equations on paper ...

The non analytical problem with spline is the min distance issue. This would require the roots of a 5th degree polynomial in t. Still it's very useful to know how far from the spline an object is and find the nearest point, the arc length coordinate on a track. But I remember I could solve this issue with mixed analysis and quick numerical approximations. Other analytically impossible things are defining "parallel" splines, or determining the distance between two splines.

Dave Eberly

Apr 12, 2005 at 9:41 PM

Quote:
Original post by Charles B
All in all you end up with the integration of the square root of a degree 4 polynomial in t.

This is correct.

Quote:
Original post by Charles B
But this can easilly be solved formally :

Go there http://torte.cs.berkeley.edu:8010/tilu and you'll get the formal answer you need. Select Mathematica as language and type (indefinite integration : finding a primitive):

sqrt(a*x^4+b*x^3+c*x^2+d*x+e)

You'll get (control by differentiating the result, it's easy):
             5          4          2          3sqrt (1/5 a x  + 1/4 b x  + 1/2 d x  + 1/3 c x  + e x)

This is incorrect. You cannot formally switch the order of integration and evaluation of sqrt. The derivative of your proposed antiderivative is

  (1/2)*(a*x^4+b*x^3+c*x^2+d*x+e)/sqrt(a*x^5/5+b*x^4/4+c*x^3/3+d*x^2/2+e*x)

which is not equal to what you started with.

The integral of the square root of a fourth-degree polynomial does not have a closed form solution in terms of a finite number of elementary functions. In computer terms, this means you must rely on a numerical integrator to estimate the value of the integral. Look in a text book on numerical methods (or in Numerical Recipes) for various integrators.

averisk

Apr 12, 2005 at 9:54 PM

Wow I didn't realize it was that easy

I'm glad I saw this thread :)

Edit: Damn, apparently I was right

John Schultz

Apr 12, 2005 at 9:56 PM

Quote:
Original post by iNsAn1tY
Evening all.

I'm doing a bit of work with 3D cubic Bezier curves, and I've run into a bit of a problem. I need to know (ideally) the exact length of the curve I'm working with. I know that I can get an approximation of the length by creating straight lines along the length of the curve, then adding together their lengths, but I was wondering if there was a simpler solution in the form of an equation that I can solve for the curve. I can't seem to find much on Google about it, so any input is greatly appreciated...

Hi,

I haven't looked at this code in few years (perhaps 5), but it might be the simplest method (may need some tweaks/changes; I already see some optimizations involving length/sqrt (sqrtInverse(), removal of divides)):

// ^ = dot productflt bezLength3(const Vec3 b[4]) {  Vec3 p0 = (b[0] - b[1]);  Vec3 p1 = (b[2] - b[1]);  Vec3 p2;  Vec3 p3 = (b[3] - b[2]);  flt l0 = p0.length();  flt l1 = p1.length();  flt l3 = p3.length();  if (l0 > 0.f) p0 /= l0;  if (l1 > 0.f) p1 /= l1;  if (l3 > 0.f) p3 /= l3;  p2 = -p1;  flt a = fastAbs(p0 ^ p1) + fastAbs(p2 ^ p3);  if ((a > 1.98f) || ((l0+l1+l3) < ((4.f - a)*8.f))) {    return l0+l1+l3;  } // if  Vec3 bl[4],br[4];  bl[0] = b[0];  bl[1] = (b[0]+b[1])*.5f;  Vec3 mid = (b[1]+b[2])*.5f;  bl[2] = (bl[1]+mid)*.5f;  br[3] = b[3];  br[2] = (b[2]+b[3])*.5f;  br[1] = (br[2]+mid)*.5f;  br[0] = (br[1]+bl[2])*.5f;  bl[3] = br[0];  return bezLength3(bl) + bezLength3(br);} // bezLength3

Charles B

Apr 13, 2005 at 3:59 AM

Quote:
Original post by Dave Eberly
This is incorrect. You cannot formally switch the order of integration and evaluation of sqrt.

Damn you are right. Been faulty and misleaded. Still for sure, this was not my idea. But it's certainly what this buggy program did ...

Quote:
The derivative of your proposed antiderivative is
  (1/2)*(a*x^4+b*x^3+c*x^2+d*x+e)/sqrt(a*x^5/5+b*x^4/4+c*x^3/3+d*x^2/2+e*x)
which is not equal to what you started with.

The integral of the square root of a fourth-degree polynomial does not have a closed form solution in terms of a finite number of elementary functions. In computer terms, this means you must rely on a numerical integrator to estimate the value of the integral. Look in a text book on numerical methods (or in Numerical Recipes) for various integrators.

... I wrote that quickly, and did only pseudo "verify" mentally. It's because I did put to much confidence in this site from berkeley.edu (???). I got it from links that seemed to be serious. And else I don't have my math books where I am currently. Do you know a serious www link to do quick formal computations ?

Charles B

Apr 13, 2005 at 4:39 AM

OK got it. In mathematica, style, type Sqrt[x] not sqrt(x). Probably this was the problem in Tilu giving me this wrong answer. Anyway I found the link to the Wolfram integrator ( http://integrals.wolfram.com/ ) again, and this confirms what was my first idea.

It should be possible to slice the Bezier into quadratic subparts (instead of straight line segments) that should approximate the curve very well locally. Thus more precision for less slices. This leads to dx/dt 1st order thus ds/dt square root of a 2nd order polynomial :

Sqrt[a*x*x+b*x+c]

has a formal solution (I hope I am right this time). It involves sqrts (which counts for nothing with SIMD) and a log. But all these functions can be hardcoded. And I suppose it will be faster o compute 10 quadratic lengths than 1000 line segments for the same precision. One could possibly add criteria on the curvature to know how to slice for constant precision and/or add some recursion for adaptative subdivision.

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.