非等间距离散点的曲率求解（Matlab）

环境：MATLAB

一、离散点曲率计算原理

参考：计算离散点的曲率（附Python, MATLAB代码） - 知乎

 二、二维离散点曲率matlab求解

1.三点法：

function [kappa,norm_k] = Solve2Dcurvature(x,y)
    %二维曲率标量形式kappa=(x"y'-x'y")/(((x')^2+(y')^2).^(3/2))，
    %三点构成的曲线参数设为x=a(1)+a(2)*t+a(3)*t^2，y=b(1)+b(2)*t+b(3)*t^2,
    %其中在t=0时对应中间那个点（x(2),y(2)），t=-t_a时对应第一个点（x(1),y(1)），t=t_b时对应第三个点（x(3),y(3)）
    x = reshape(x,3,1);
    y = reshape(y,3,1);
    t_a = norm([x(2)-x(1),y(2)-y(1)]);
    t_b = norm([x(3)-x(2),y(3)-y(2)]);
    M =[[1, -t_a, t_a^2];
        [1, 0,    0    ];
        [1,  t_b, t_b^2]];
    a = M\x;%M的逆乘以x
    b = M\y;
    %求中间点的曲率kappa
    kappa  = 2.*(a(3)*b(2)-b(3)*a(2)) / (a(2)^2.+b(2)^2.)^(1.5);
    norm_k =  [b(2),-a(2)]/sqrt(a(2)^2.+b(2)^2.);
    
%{
    % n=i（i=1，2，3）是取第i个点为最终的曲率值
    xd1 = a(2)+2*a(3)*[-t_a,0,t_b];%x一阶导
    xd2 = 2*a(3);%x二阶导
    yd1 = b(2)+2*b(3)*[-t_a,0,t_b];
    yd2 = 2*b(3);
    n=2;
    Curvature = (xd1(n)*yd2 - xd2*yd1(n))./((xd1(n).^2+yd1(n).^2).^(3/2));
%}
end

点多时（速度V、加速度A可得时） 

function [kappa] = VA2Curvature(v,a)
    %input:二维平面速度v:2X1;加速度a:2X1
    %二维曲率标量形式kappa=(x"y'-x'y")/(((x')^2+(y')^2).^(3/2))，
    xd1 = v(1,1);%x一阶导
    xd2 = a(1,1);%x二阶导
    yd1 = v(2,1);
    yd2 = a(2,1);
    n=1;
    kappa = (xd1*yd2 - xd2*yd1)./((xd1.^2+yd1.^2).^(3/2));%kappa = (x"y'-x'y")/(((x')^2+(y')^2).^(3/2))
end

2.三点法效果：

         取相隔点数为d的三个点计算中间点的曲率：

d = 5;%选取用于计算曲率的曲线上的三个点
for i=d+1:column-d
    xi=[pos(1,i-d);pos(1,i);pos(1,i+d)];
    yi=[pos(2,i-d);pos(2,i);pos(2,i+d)];
    [cur(i),norm_cur(:,i)]=Solve2Dcurvature(xi,yi);
end

离散点轨迹图：

 对应曲率：

 

 

 PS:曲率计算效果受d值（即三点选取策略）影响很大。

三、三维离散点曲率求解

高斯曲率：参考 (22条消息) 主曲率，主方向，高斯曲率与平均曲率公式与matlab代码_zb12138的博客-CSDN博客_matlab 高斯曲率

持续更新ing...