No BSD License  

Highlights from
curvature to cartesian

image thumbnail
from curvature to cartesian by skouliki katelouzos
with this function you can convert a (curvature- curve length) relation to Cartesian coordinates.

[xo yo]=curv2cart(a,b,c,d,e,l,p)
function [xo yo]=curv2cart(a,b,c,d,e,l,p)
%This function it is used to transform a (curvature-curve length) relation 
%to cartesian coordinates.k=f(s)->x=x(t),y=y(t)
%usage:[xo yo]=curv2cart(a,b,c,d,e,l,p)
%'a','b','c','d' and 'e' are part of the polynomial
%relation : k=a*s^4+b*s^3+c*s^2+d*s+e;
%l: length of the curve in points
%p: precision.Distance between two points.
%
%EXAMPLES:
%
%curv2cart(0,0,0,0,c,10,0.001);         %plots a circle of radius c
%curv2cart(0,0,0,0,c,10,0.1);           %plots a badly calculated circle due to low
%precision
%
%Interesting and reason for existance of this function, are the polynomial
%spirals:
%
%curv2cart(0,0,0,1,0,10,0.001);         %k=s
%curv2cart(0,0,1,0,0,10,0.001);         %k=s^2
%curv2cart(0,0,1,0,-2.19,10,0.001);     %k=s^2-2.19
%curv2cart(0,0,1,0,-4,10,0.001);        %k=s^2-4
%curv2cart(0,0,1,0,1,10,0.001);         %k=s^2+1
%curv2cart(5,0,-18,0,5,4,0.001);        %k=5s^4-18s^2+5
%
%curv2cart, creates some cool pseudofractals:
%
%curv2cart(0,0,2,0,-1,50,0.01);

s=0:p:l;
s2=-s;
x1=zeros(1,length(s));
y1=zeros(1,length(s));
x2=zeros(1,length(s));
y2=zeros(1,length(s));
k=a*s.^4+b*s.^3+c*s.^2+d*s+e;
k2=a*s2.^4+b*s2.^3+c*s2.^2+d*s2+e;
x1(2)=s(2);
y1(2)=0;
x2(2)=s2(2);
y2(2)=0;

for i=3:(length(s))
    dth=k(i)*s(2);
    cth=(x1(i-1)-x1(i-2))/norm([x1(i-1) y1(i-1)]-[x1(i-2) y1(i-2)]);
    sth=(y1(i-1)-y1(i-2))/norm([x1(i-1) y1(i-1)]-[x1(i-2) y1(i-2)]);
    if sth>=0
        scorr=1;
    else scorr=-1;
    end
    newth=dth+scorr*acos(cth);
    xx=cos(newth);
    yy=sin(newth);
    x1(i)=(xx*s(2))+x1(i-1);
    y1(i)=(yy*s(2))+y1(i-1);
end

for i=3:(length(s))
    dth=k2(i)*s2(2);
    cth=(x2(i-1)-x2(i-2))/norm([x2(i-1) y2(i-1)]-[x2(i-2) y2(i-2)]);
    sth=(y2(i-1)-y2(i-2))/norm([x2(i-1) y2(i-1)]-[x2(i-2) y2(i-2)]);
    if sth>=0
        scorr=1;
    else scorr=-1;
    end
    newth=dth+scorr*acos(cth);
    xx=cos(newth);
    yy=sin(newth);
    x2(i)=(xx*s(2))+x2(i-1);
    y2(i)=(yy*s(2))+y2(i-1);
end
xo=[fliplr(x2) x1];
yo=[fliplr(y2) y1];
plot(xo,yo);

Contact us at files@mathworks.com