views:

272

answers:

2

Using scipy's interpolate.splprep function get a parametric spline on parameter u, but the domain of u is not the line integral of the spline, it is a piecewise linear connection of the input coordinates. I've tried integrate.splint, but that just gives the individual integrals over u. Obviously, I can numerically integrate a bunch of Cartesian differential distances, but I was wondering if there was closed-form method for getting the length of a spline or spline segment (using scipy or numpy) that I was overlooking.

Edit: I am looking for a closed-form solution or a very fast way to converge to a machine-precision answer. I have all but given up on the numerical root-finding methods and am now primarily after a closed-form answer. If anyone has any experience integrating elliptical functions or can point me to a good resource (other than Wolfram), That would be great.

I'm going to try Maxima to try to get the indefinite integral of what I believe is the function for one segment of the spline: I cross-posted this on MathOverflow

+1  A: 

Because both x & y are cubic parametric functions, there isn't a closed solution in terms of simple functions. Numerical integration is the way to go. Either integrating the arc length expression or simply adding line segment lengths - depends on the accuracy you are after and how much effort you want to exert.

An accurate and fast "Adding length of line segments" method:

Using recurvise subdivision (a form of de Casteljeau's algorithm) to generate points, can give you a highly accurate representation with minimal number of points. Only subdivide subdivisions if they fail to meet a criteria. Usually the criteria is based on the length joining the control points (the hull or cage). For cubic, usually comparing closeness of P0P1+P1P2+P2P3 to P0P3, where P0, P1, P2 & P3 are the control points that define your bezier.

You can find some Delphi code here: link text

It should be relatively easy to convert to Python. It will generate the points. The code already calculates the length of the segments in order to test the criteria. You can simply accumulate those length values along the way.

symmetry
This solution is nice for approximations, but I should have stated in my question that I am after an exact or a machine-precision answer. It is similar to EOL's romberg solution, but i would have to iterate in pure python.
bpowah
As the general analytic solution involves elliptic integrals you will be numerically integrating no matter what approach you take.The method I have outlined doesn't directly integrate the arc length expression but it is a numeric integration as well. All that matters is that the method converges to a desired accuracy in a desired time.
symmetry
To be a bit more precise about those lengths I was referring to that you would accumulate along the way...have a look here:http://steve.hollasch.net/cgindex/curves/cbezarclen.htmlUse the average of the hull length (L1) and the chord length (L0) to approximate the arc-length of the individual subsegments.
symmetry
+1  A: 

You can integrate the function sqrt(x'(u)**2+y'(u)**2) over u, where you calculate the derivatives x' and y' of your coordinates with scipy.interpolate.splev. The integration can be done with one of the routines from scipy.integrate (quad is precise [Clenshaw-Curtis], romberg is generally faster). This should be more precise, and probably faster than adding up lots of small distances (which is equivalent to integrating with the rectangle rule).

EOL
I like this solution as you are able to set your desired accuracy, but unfortunately I am after machine precision-type accuracy. I implemented your suggestion, but QUADPACK (as awesome as it is) doesn't meet my performance requirements. I am currently seeking a closed-form solution and have cross-posted on MathOverflow. BTW: Clenshaw-Curtis is only exact/precise for polynomials. The length of this spline is an elliptical function.
bpowah