Is it possible to get signed curvature from splines of a non-single-valued curve?

Question

Jon 2016-2-13

2
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/267933-is-it-possible-to-get-signed-curvature-from-splines-of-a-non-single-valued-curve

编辑： Jon 2016-2-18

testcl.mat

I am trying to compute the curvature of a meandering river centerline. I have been using a method of finding the angle of the centerline (as we move along it) with respect to the horizontal using forward differencing (i.e. theta = arctan(y2-y1/x2-x1), then forward differencing again to find curvature (i.e. kappa = (theta2-theta1)/ds). This is only one way to compute curvature with discrete differencing; I have used two others. The resulting curvature signals using discrete differencing like this are very noisy, and require smoothing (something I'd like to avoid if possible).

However, the problem is that the curvature signal that results from these methods is sensitive to the discretization of the centerline. I.e., if the centerline is sampled more coarsely, then the average of the curvature signal will decrease. I would like a way to obtain the curvature signal that is (more?) independent of the sampling density. First question, is this possible? My idea is to use splines, but I can't recover a curvature signal from them that is comparable to the discrete methods.

The centerline is just an array of (x,y) coordinates in space. My approach so far has been to fit splines to both the x and y vectors; just the x-vectors is shown below:

t = 1:length(cl);
xs = cl(:,1);
ys = cl(:,2);
% The spline that interpolates (x,y)
sxs = spline(t,xs);
% first derivative, pp form
s1xs = sxs;
s1xs.order = sxs.order-1; 
s1xs.coefs = bsxfun(@times, sxs.coefs(:,1:end-1), s1xs.order:-1:1);
% second derivative, pp form
s2xs = s1xs;
s2xs.order = s1xs.order-1; 
s2xs.coefs = bsxfun(@times, s1xs.coefs(:,1:end-1), s2xs.order:-1:1);
% Evaluate the curvature
A1xs = ppval(s1xs,t); % first derivative at x
A2xs = ppval(s2xs,t); % second derivative at x

then curvature is computed as

kappa = A2xs./(1 + A1xs.^2).^(3/2);

I can do the same for the y-vector, but then I am not sure how to "combine" the curvatures of the x and y vectors into a "total" curvature. Can I simply do a sum-of-squares, or is it more complicated? Or am I approaching this problem incorrectly? Finally, I would also like to retain the sign of curvature. That is, when the concave curvatures are positive and convex are negative. A sum-of-squares approach would erase the sign of the curvature...

Any thoughts are appreciated! I am attaching a .mat file of a sample centerline.