Curve curvature display

Use

load matlab

x = locusx(1:5:end)';
y = locusy(1:5:end)';
[u,v] = PJ(x,y);
figure(1)
quiver(x(2:end-1),y(2:end-1),u,v,1)
hold on
plot(x,y)
axis equal

Subfunctions

function [kappa,norm_k] = PJcurvature(x,y)
    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;
    b = M\y;

    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.);
    
    
end

function [u,v] = PJ(x,y)
kappa_arr = zeros(length(x)-2,1);
norm_arr = zeros(length(x)-2,2);
    for i=2:(length(x)-1)
        
        
        [kappa,norm_l] = PJcurvature(x(i-1:i+1),y(i-1:i+1));
        kappa_arr(i-1) = kappa;
        norm_arr(i-1,:) = norm_l;
    end
    u = kappa_arr.* norm_arr(:,1);
    v = kappa_arr.* norm_arr(:,2);
end