Hello guys,

I am in a desperate need of help.

Situation: I have a set of data (1D graph) for which I want to calculate curvatures.

In order to calculate curvatures I'm using Catmullrom spline interpolation
f(x)' = 0.5*(-p0 + p2) + (2*p0-5*p1+4*p2-p3) * t + 1.5 * (-p0 + 3*p1-3*p2 + p3)*t*t;
f(x)" = (2*p0-5*p1+4*p2-p3) + 3 * (-p0 + 3*p1-3*p2 + p3)*t;

1D (graph) curvature formula is: y" / (1+y'^2)^3/2

Sample data: 100, 90, 100, 1, 70, 100

What I expected to get is the highest curvature © value for 1 (sample data), however, the result is:

p0: 110, p1: 100, p2: 90, p3: 100, t: 0 x': -10 , x": -20 , c: 0.0197037
p0: 100, p1: 90, p2: 100, p3: 1, t: 1.49012e-008 x': 2.22027e-006 , x": 149 , c: 149 <--- wasn't expected
p0: 90, p1: 100, p2: 1, p3: 70, t: 2.98023e-008 x': -44.5 , x": -386 , c: 0.00437701
p0: 100, p1: 1, p2: 70, p3: 100, t: 1.19209e-007 x': -15 , x": 375 , c: 0.110375 <--- expected a huge spike here
p0: 1, p1: 70, p2: 100, p3: 130, t: 5.96046e-008 x': 49.5 , x": -78 , c: 0.000642707
p0: 100, p1: 1, p2: 70, p3: 100, t: 1 x': 49.5 , x": -246 , c: 0.002027

From the table above we can see that t is almost 0, thus
f(x)' may be rewritten as 0.5*(-p0 + p2) <--- relies only on two adjacent control points!
f(x)" would be then (2*p0-5*p1+4*p2-p3)

So we modify our set of data to manipulate f(x)'.
Data set: 100, 90, 90, 1, 90, 100

p0: 110, p1: 100, p2: 90, p3: 90, tt: 0 x': -10 , x": -10 , c: 0.00985185
p0: 100, p1: 90, p2: 90, p3: 1, tt: 1.49012e-008 x': -5 , x": 109 , c: 0.82218 <---- ok
p0: 90, p1: 90, p2: 1, p3: 90, tt: 2.98023e-008 x': -44.5 , x": -356 , c: 0.00403683
p0: 90, p1: 1, p2: 90, p3: 100, tt: 1.19209e-007 x': 5.1856e-005 , x": 435 , c: 435 <--- works as expected
p0: 1, p1: 90, p2: 100, p3: 110, tt: 5.96046e-008 x': 49.5 , x": -158 , c: 0.00130189
p0: 90, p1: 1, p2: 90, p3: 100, tt: 1 x': 49.5 , x": -336 , c: 0.00276858


My final question: is it possible to rely on curvature calculation using catmullrom interpolated spline. If yes, what is wrong with my approach? Currently it seems that the correctness of the curvature calculation depends solely on my set of data and it is absolutely unreliable.

This problem makes me dizzy for the third day now. Don't know what to do...
Would be very grateful if you could help me! Thank you very much in advance!

Ginga

PS. I was encouraged to post here by this related thread: http://www.gamedev.net/topic/518279-mesaurement-of-curvature-of-camull-rom-spline-segment/

Your derivatives seem to be incorrect:

q(t) = 0.5 * ( (2 * P1) +
(-P0 + P2) * t +
(2*P0 - 5*P1 + 4*P2 - P3) * t^2 +
(-P0 + 3*P1 - 3*P2 + P3) * t^3)

q'(t) = 0.5 * ( (-P0 + P2) +
2 * (2 * P0 - 5*P1 + 4*P2 - P3) * t +
3 * (-P0 + 3*P1 - 3*P2 + P3) * t^2)

Second derivative should be easy to derive from this.

If we open the brackets of your provided derivative, we will get the same as mine:

q'(t) = 0.5 * ( (-P0 + P2) + 2 * (2 * P0 - 5*P1 + 4*P2 - P3) * t +3 * (-P0 + 3*P1 - 3*P2 + P3) * t^2)=
0.5*(-P0 + P2) + (2 * P0 - 5*P1 + 4*P2 - P3) * t + 1.5 * (-P0 + 3*P1 - 3*P2 + P3) * t^2)


But thank you for the quick answer!

You're right, didn't recognise it in that form and didn't give it a second thought. Not sure what you're trying to accomplish, but plotting the graphs might help you see where the problem lies.

I want to express my data set as a variation. Curvature would suit this purpose perfectly. In order to calculate curvatures I must get f' and f'', hence I need a formula to get derivatives. Catmullrom seems to be he simplest solution.

Recently I found that Catmullrom's 2nd derivative (curvature. 1st is slope) is not continuous (linear). I wonder if this may be the case for my problem.

I am not a very strong mathematician.

Thought to experiment with NURBS, but have no idea how to get 2nd deriv...

Would a tool like Wolfram Alpha help you to find the derivative that you're looking for?

http://www.wolframalpha.com/input/?i=second+derivative+of+x^2

Recently I found that Catmullrom's 2nd derivative (curvature. 1st is slope) is not continuous (linear). I wonder if this may be the case for my problem.

You can use cubic splines that have continuous second derivative, and one usually makes the second derivative at the ends be 0 (look up "natural cubic splines").

I still don't quite understand the importance of the curvature, though. Perhaps you can elaborate.

It is possible to calculate discrete analogous of the curvature without computing any spline curve. Depending on your needs, it may works for you.

You may know that the curvature at a point of a curve is one over the radius of the osculating circle at that point. A simple discrete approximation of the curvature can thus be found by computing the radius of the circle intersecting each set of three consecutive points in the data set. The actual formula for this can be easily found with google.

Another discrete analogous of the curvature can be found by thinking at the curvature as the rate of change of the angle between the tangent and the x-axis. If a and b are the length of the segments between three consecutive points and theta is the angle between the vector from the first point to the second with the vector from the second point to the third, then the curvature can be defined as the quantity 2*theta / (a + b).

I have experimented a bit with your example and probably understood why the curvature is so much different from what you expected. My previous methods (and probably any method similar to what has been presented in this thread) should probably have the same behaviour. When you trace a smooth curve between three points, the angle theta is smoothed in all the length of the curve and the curvature at the centre point is thus reduced based on how long the segments between the points are. This is probably correct, even if it isn't what you would expect/want.

A simple solution to this problem, based on the last method I have presented in my last post, is to assume a + b = 2 (that is we assume the distances between adjacent points to be of length one) and thus define the curvature as the angle theta. Note that this isn't a correct approximation of the curvature since it doesn't converge to the correct one. However it can be probably useful. In formulas you have that the angle is given by atan2( h*(2*y1 - y0 - y2), h^2 + (y2 - y1)*(y1 - y0) ) (this is based on the general formula atan2( perp_dot(u,v), dot(u,v) )). y0, y1 and y2 are the adjacent values in the array and h is the distance along the x-axis between adjacent samples.