我就废话不多说了，大家还是直接看代码吧~
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
|
x0 = linspace(0.1,2,100);%x0,y0验证函数离散点，可以非等间隔y0 = 1./x0;h1 = abs(diff(x0)) ;h = [h1 h1(end)];ht = h;yapp1 = gradient(y0)./ht; %matlab数值近似yapp2 = del2(y0)./ht; %matlab数值近似k2 = abs(yapp2)./(1+yapp1.^2).^(3/2);figureplot(k2)title('曲率曲线')[~,maxflag] = max(k2);%曲率最大位置x_max = x0(maxflag);y_max = y0(maxflag);%画出图像 标注曲率最大点figureplot(x0,y0,'.-');hold on;plot(x_max,y_max,'rp')title('标注最大曲率点')xlabel('log10((norm(b*xk-l)))')ylabel('log10((norm(xk)))') |
补充：matlab 插值+计算离散点曲率
思路：点足够密的话直接用 diff、gradient 求曲率，稀疏的话先插值再算曲率。
公式：
点密的情况 输入曲线坐标（1-2）求一、二阶导数（4-9）通过公式求得曲率（10）
|
1
2
3
4
5
6
7
8
9
10
11
|
x = 0:0.01:7;y = cos(x*0.5*pi); h1 = abs(diff(x));h = [h1 h1(end)];ht = h; y1 = gradient(y)./ht;y2 = gradient(y1)./ht;curv = abs(y2)./sqrt((1+y1.^2).^3); plot(x,y,'-',x,curv,'--r);legend('raw data, 'curvature','location',"best");grid on |
图像与下文理论值图像相同
点稀疏的情况
1、输入散点坐标（1-2）
2、用样条曲线（b-spline）等方法插值得到拟合曲线（3-4）
3、diff、gradient 函数求拟合曲线的一、二阶导数（6-11）
4、通过公式求得曲率（12）
例：余弦函数取 8 个点，用 b-spline 插值
|
1
2
3
4
5
6
7
8
9
10
11
12
13
|
x = 0:1:7;y = cos(x*0.5*pi);xx = 0:0.01:7;yy = spline(x,y,xx); h1 = abs(diff(xx));h = [h1 h1(end)];ht = h; yy1 = gradient(yy)./ht;yy2 = gradient(yy1)./ht;curv = abs(yy2)./sqrt((1+yy1.^2).^3); plot(xx,yy,'-',xx,curv,'--r',x,y,'o-');legend('b-spline', 'curvature','raw data','location',"best");grid on |
补充用法
求最大曲率并在图中标出
|
1
2
3
|
[max_val,max_ind]=max(curv);hold onplot(xx(max_ind),yy(max_ind),'*r'); |
与理论值（余弦函数曲线）对比
曲率对比
几种插值方法对比
列举四种方法，分别为：分段线性插值、三次样条曲线（b-spline）插值、三次 hermite 插值（pchip）、修正 akima 分段三次 hermite 插值（akima）
case 1： 三维螺线
三维螺线散点
插值
俯视
侧视
case 2：二维梯形波
二维梯形波
case 3：三维不规则折线
三维不规则折线（不等间距）
对比可得：
case 1：b-spline>akima>pchip>linear
case 2：linear>pchip>akima>b-spline
case 3：linear≈pchip≈akima>b-spline
故在插值的时候需要选择适合的计算方法
以上为个人经验，希望能给大家一个参考，也希望大家多多支持服务器之家。如有错误或未考虑完全的地方，望不吝赐教。
原文链接：https://blog.csdn.net/qq_38991255/article/details/85127468
