Radius of Curvature Array

The code generates radius of curvatures for three points(coordinates x1,y1 x2,y2 and x3,y3. However, i have a set of coordinates (77) in all, and I am trying to rewrite the code to loop through the coordinates and generate a radius of curvature for every three points within the set of 77, and store the result in another variable say R. I have not been able to do this and I need urgent help: See the radius of curvature code below that generates when radius of curvature for 3 points.

function lambda=mycurvature(xinput, yinput)
global xstar ystar;

xstar=xinput;ystar=yinput;
alpha=atan2(ystar(2)-ystar(1),xstar(2)-xstar(1))-atan2(ystar(3)-ystar(2),xstar(3)-xstar(2));

if abs(alpha)<0.002;                       %set the curvature to 0 if the points are aligned
R=inf; lambda=0;
%plot(xstar,ystar,'-r',xstar,ystar,'*b','linewidth',2,'markersize',10);
else
initialguess = [1; 1; 1];                % Make a starting guess at the solution
options=optimset('Display','none');      % Option to display output
res = fsolve(@userfun,initialguess,options);    % Call optimizer
x0=res(1);y0=res(2);R=res(3);     %x-coord of center; y-coord of center; radius of curvature
lambda=1/R;           %curvature

%t=0:0.01:2*pi;xx=x0+R*cos(t);yy=y0+R*sin(t);
%plot(xx,yy,'-r',xstar,ystar,'*b','linewidth',2,'markersize',10);
end

%str=sprintf('radius=%6.3f, curvature=%6.3f',R,lambda);
%title(str, 'fontsize', 14);xlabel('x'),ylabel('y');
%axis equal;


function F = userfun(vars)
global xstar ystar;
F = [...
(xstar(1)-vars(1))^2+(ystar(1)-vars(2))^2-(vars(3))^2;  ...
(xstar(2)-vars(1))^2+(ystar(2)-vars(2))^2-(vars(3))^2;  ...
(xstar(3)-vars(1))^2+(ystar(3)-vars(2))^2-(vars(3))^2;  ...  
];

xinput are the array of x coordinates for the three points, yinput are the array of y coordinates for the 3 points.

What i am doing that is not working is

x = [array of my x-coordinates of my 77 points]; y =[array of my y-coordinates of my 77 points]
for i = 1:3:length(x)-1
xinput = x(i), x(i+1), x(i+2)
yinput = y(i), y(i+1), y(i+2)

and then removed the function at the first line of code for the radius of curvature code above. This did not work. What am I doing wrong? Please help.

Answer 1

[darova]

Try not to use global variables if its possible

function [x0, y0, lambda] = mycurvature(xinput, yinput)
alpha=atan2(yinput(2)-yinput(1),xinput(2)-xinput(1))-atan2(yinput(3)-yinput(2),xinput(3)-xinput(2));
if abs(alpha)<0.002;                       %set the curvature to 0 if the points are aligned
    R=inf; lambda=0;
    %plot(xinput,yinput,'-r',xinput,yinput,'*b','linewidth',2,'markersize',10);
else
initialguess = [1; 1; 1];                % Make a starting guess at the solution
options=optimset('Display','none');      % Option to display output
res = fsolve(@(vars)userfun(vars,xinput,yinput),initialguess,options);    % Call optimizer
x0=res(1);y0=res(2);R=res(3);     %x-coord of center; y-coord of center; radius of curvature
lambda=1/R;           %curvature
%t=0:0.01:2*pi;xx=x0+R*cos(t);yy=y0+R*sin(t);
%plot(xx,yy,'-r',xinput,yinput,'*b','linewidth',2,'markersize',10);
end
%str=sprintf('radius=%6.3f, curvature=%6.3f',R,lambda);
%title(str, 'fontsize', 14);xlabel('x'),ylabel('y');
%axis equal;
function F = userfun(vars,xinput,yinput)
F = [...
(xinput(1)-vars(1))^2+(yinput(1)-vars(2))^2-(vars(3))^2;  ...
(xinput(2)-vars(1))^2+(yinput(2)-vars(2))^2-(vars(3))^2;  ...
(xinput(3)-vars(1))^2+(yinput(3)-vars(2))^2-(vars(3))^2;  ...  
];

Main code:

clc , clear
a = 2;
b = 2;
t = linspace(0,pi/2,30);
x = a*cos(t);
y = b*sin(t);
plot(x,y,'.-r')
hold on
for i = 2:length(x)-1
    xinput = [x(i-1) x(i) x(i+1)];
    yinput = [y(i-1) y(i) y(i+1)];
    [xc, yc, lam] = mycurvature(xinput,yinput);
    plot([xc x(i)],[yc y(i)],'.-b')
end
hold off
grid on

[darova]

Your code didn't work because of (i think)

xinput = x(i), x(i+1), x(i+2)
yinput = y(i), y(i+1), y(i+2)

Not because of global variables

[Jide Williams]

Thanks darova, I really appreciate your response. However, running the code the way you instructed still did not return an array of radius of curvature. Look at, what do you think?

[darova]

You have tou create then those arrays for curvate

lambda = zeros(length(xinput)-2,1);
for i = 1:length(x)-2
    xinput = [x(i) x(i+1) x(i+2)];
    yinput = [y(i) y(i+1) y(i+2)];
    [xc, yc, lambda(i)] = mycurvature(xinput,yinput);
end
disp(lambda)

[Jide Williams]

Thanks darova for the effort, well, maybe I should give you my set of x,y's and see for yourself.

x = [0.2 0.4 0.6 0.8 1 1.1 1.2 1.4 1.7 1.9 2.2 2.6 2.9 3.3 4 4.3 4.71 4.99 5.295 5.4 5.8 6.2 6.6 6.9 7.19 7.5 7.59 7.5 7.4 7.3 7.1 6.7 6.4 6.0 5.6 5.3 5.1 4.8 4.6 4.5 4.5 4.6 4.7 4.9 5.1 5.2 5.4 5.65 6 6.4 6.7 7.0 7.2 7.6 7.9 8.2 8.51 8.9 9.1 9.4 9.6 9.9 10.2 10.4 10.8 11.30 11.6 12.1 12.4 12.8 13 13.4 13.6 14 14.2 14.7 15];

y = [8 8 8 8 8 3.4 3 2.5 2.2 1.9 1.6 1.3 1.1 0.8 0.8 0.8 0.8 0.8 0.87 0.87 1.0 1.2 1.4 1.6 1.9 2.3 2.7 3.3 3.7 3.9 4.2 4.4 4.5 4.5 4.4 4.2 4 3.6 3.3 2.9 2.5 2.0 1.7 1.4 1.2 1.09 1 0.8 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 1.1 1.3 1.6 1.8 2.1 2.4 2.6 2.8 2.8 2.7 2.4 2.29 2.1 2 1.7 1.6 1.3 1.2 1.2 1.2];

with these values, see if you can generate the radius of curvature array. Be mindful that 3 points will have thesame radius of curvature.

Thanks

[Jide Williams]

I found a way out of it. Thanks Darova!!!!

[Divya]

Y I am getting a errors

[Divya]

please reply

Answer 2

[Jide Williams]

the code isn't wrong xinput = x(i), x(i+1), x(i+2); and y(i), y(i+1), y(i+2); it just that its not looping all through the 77 points. It should loop through all the 77 points pick three points at a time and finding the Radius of curvature for three points at a time.