曲线的曲率特征

 找了很多计算曲率的例子，都是x轴等间距的例子，我自己写了一个非等间距的曲率计算方法，将一个sin函数偏移并且旋转一定角度，分别测量他们的曲率和与二阶导的和

clc;
clear;

X=(0:0.001:2*pi);
gradient=10;
x0=sin(X');
x0=[X',x0];

%旋转点
alpha=pi/6;
x1=[(x0(:,1)+1)*cos(alpha)-x0(:,2)*sin(alpha),(x0(:,2)+1)*cos(alpha)+x0(:,1)*sin(alpha)];
%添加白噪声
%x1(:,2)=awgn(x1(:,2),80);
size0=size(x0,1);
height=size0(1);

%求得曲率累积
dist_cur0=StandardErrorPoint(gradient,x0);
dist_cur1=StandardErrorPoint(gradient,x1);

%求得二阶梯度累积
dist_gra0=StandardError(0,10,x0(:,2));
dist_gra1=StandardError(0,10,x1(:,2));

%显示
label0=strcat( '\leftarrow ',num2str(dist_cur0),'+e3');
label1=strcat( '\leftarrow ',num2str(dist_cur1),'+e3');
plot(x0(:,1),x0(:,2),'-r');
text(x0(height,1),x0(height,2),label0);
text(x1(height,1),x1(height,2),label1);
hold on
plot(x1(:,1),x1(:,2),'-g');

 下面是计算曲率的和的函数

function [curvature]=StandardErrorPoint1(gradient,x0)
    size_x=size(x0,1);    
    div=zeros(size_x,1);
    %一阶导
    for i=gradient+1:size_x-gradient
        div(i)=(x0(i+gradient,2)-x0(i,2))/(x0(i+gradient,1)-x0(i,1));
    end
    for i=1:gradient
       div(i)=div(gradient+1); 
    end
    for i=size_x-gradient+1:size_x
       div(i)=div(size_x-gradient); 
    end  
    %二阶导
    div2=zeros(size_x,1);
    for i=gradient+1:size_x-gradient
        div2(i)=(div(i+gradient)-div(i))/(x0(i+gradient,1)-x0(i,1));
    end
    for i=1:gradient
       div2(i)=div2(gradient+1); 
    end
    for i=size_x-gradient+1:size_x
       div2(i)=div2(size_x-gradient); 
    end
    %曲率
    cuv=zeros(size_x,1);
    for i=gradient+1:size_x-gradient
        d10=(x0(i-gradient,2)-x0(i,2))/(x0(i-gradient,1)-x0(i,1));
        d11=(x0(i+gradient,2)-x0(i,2))/(x0(i+gradient,1)-x0(i,1));
        d2=(d11-d10)/(x0(i+gradient,1)-x0(i,1));
        cuv(i)=abs(d2)/(1+d11*d11).^(3/2);
        %cuv(i)=abs(div2(i))/(1+div(i)*div(i)).^(3/2);
    end
    for i=1:gradient
       cuv(i)=cuv(gradient+1); 
    end
    for i=size_x-gradient+1:size_x
       cuv(i)=cuv(size_x-gradient); 
    end
    curvature=sum(cuv);
%    figure(3600)
%     plot(div,'-r'); 
%     hold on
%     plot(div2,'-g');
%     hold on
%     plot(cuv,'-b');
%     hold on
%     offset=fix(x0(1,1)*100);
%     text(1000+offset,div(1000+offset),['div' num2str(x0(1,1))]);
%     text(2000+offset,div2(2000+offset),['div2' num2str(x0(1,1))]);
%     text(2000+offset,cuv(2000+offset),['cuv' num2str(x0(1,1))]);

 下面这个是计算的二阶导的和，有点问题，有空再搞搞

%统计二阶导
function dist = StandardError(crop,gradient, y)
    size0=size(y);
    y_crop=y(crop+1:size0(1)-crop);
    size_crop=size(y_crop);
    size_crop=size_crop(1);
    y_gradient=zeros(size_crop,1);
    for i=gradient:size_crop-gradient
        y_gradient(i)=abs(y_crop(i-fix(gradient)+1)+y_crop(i+fix(gradient))-2*y_crop(i));
        %y_gradient(i)=y_crop(i-fix(gradient)+1)+y_crop(i+fix(gradient))-2*y_crop(i)/(size_crop-2*gradient+1);
    end
    %平滑梯度
    filter_num=1;   
    for i=1:filter_num
        %y_gradient=Gaussianfilter(100,10,y_gradient);                     %窗体大小和sigma没影响，
        %y_gradient=Medfilter(300,y_gradient);                             %修改滤波次数,经过检测，变化不大；窗体大小变大，梯度表差越小
        %y_gradient=Meanfilter(500,y_gradient);                             %窗体变化，影响不大
    end
%     figure(100);
%     plot(y_gradient);
    dist=sum(y_gradient);
end
        

 实验效果，