How do i detect/describe a curvature from a data-set of coordinates

I have a data-set of coordinates that makes up a 2d body. It's a rectangular shape with one side curved or kinda circular. I would like to describe this curved side with a polyfit line so that I can be able to get the radius. How can I achieve this? Thank you for your very much!

[KALYAN ACHARJYA]

Mathematical language is far better communication way to clear your doubt.

Answer 1

[Star Strider]

The polyfit function will not give you the centre and radius of your arc.

This will:

t = linspace(-pi/4, pi/4, 10);                                  % Create Data
r = 2.5;                                                        % Create Data
x = r*cos(t)+1;                                                 % Create Data
y = r*sin(t)+2;                                                 % Create Data

fcn = @(b) norm(((x(:)-b(1)).^2 + (y(:)-b(2)).^2) - b(3).^2);   % Objective Function
B = fminsearch(fcn, rand(3,1));                                 % Estimate Parameters
Cx = B(1)                                                       % Center X-Coordinate
Cy = B(2)                                                       % Center Y-Coordinate
R  = B(3)                                                       % Radius
% B = fsolve(fcn, rand(3,1))

figure
plot(x, y, 'pg')
hold on
plot(Cx, Cy, 'pb')
plot([Cx x(1)], [Cy y(1)], '-r')
hold off
grid
legend('Data', 'Centre', 'Radius')
axis equal

Experiment with it to get the result you want.

[Conrade Muyambo]

Thanks for the response. But my arc is rough not uniform. How can i describe it with a 2nd order polyfit? Here is my figure

[Star Strider]

My pleasure.

If you want parameters from it such as the centre and radius, polyfit will not provide them. I seriously doubt that polyfit can approximate that anyway.

It might be best to isolate each ‘corner’ of your rectangular shape and estimate the radius and centre of each of the four arcs.

[Conrade Muyambo]

Oh i see. Then I will work with that method. How can one seperate the sides? Any codes will be greatly appreciated!

[Star Strider]

Your data have x and y coordinates, so the easiest way would be to use the find function, for example

NE = find(x>3.2 & y> 2.7);
SE = find(x>3.2 & y<1.2);
SW = find(x<1.2 & y<1);
NW = find(x< 1.2 & y>2.7);

These are approximate, and will give you the linear indices into your vectors. You then use those to define the x and y vectors for each of the four times you use my parameter estimation code.

So for instance the first one would define:

xv = x(NE);
yv = y(NE);

and so for the rest. Then use ‘xv’ and ‘yv’ for the ‘x’ and ‘y’ in my code.

I do not have your data so I cannot test my code with them. However these should work.

[Conrade Muyambo]

Thanks. I have used a similar approach as shown in my code below to find the edges or vertex points of a regular rectangular shape. But now my shape is now irregular and I would like to know the points where the arc starts and ends. respectively. I have also attached my coordinates.

ymax = max(y) ;% finding the biggest or maximum value of y
ymv1 = find(y==ymax) ;% locating its position in the array or vector (position vector)
y(ymv1)       ;         % getting the actual value of y at that position
x(ymv1)        ;        % getting the actual value of x at that position
xol=min(x(ymv1)) ;      % only considering the smallest value
yol=ymax ;

[Star Strider]

Try this:

D = load('test2-coordinates.txt');
x = D(:,1);
y = D(:,2);

[xl,xh] = bounds(x);                                                    % X-Coordinate Limits
xr = xh-xl;                                                             % X-Coordinate Range
[yl,yh] = bounds(y);                                                    % Y-Coordinate Limits
yr = yh-yl;                                                             % Y-Coordinate Range

NE = (x>(xh-0.15*xr)) & (y>(yh-0.15*yr));                               % ‘Northeast’ Corner
SE = (x>(xh-0.15*xr)) & (y<(yl+0.15*yr));                               % ‘Southeast’ Corner
SW = (x<(xl+0.15*xr)) & (y<(yl+0.15*yr));                               % ‘Southwest’ Corner
NW = (x<(xl+0.15*xr)) & (y>(yh-0.15*yr));                               % ‘Northwest’ Corner

Lvc = {find(NE),find(SE),find(SW),find(NW)};                            % Vector Cell Array

fcn = @(b,x,y) norm(((x(:)-b(1)).^2 + (y(:)-b(2)).^2) - b(3).^2);       % Objective Function
for k = 1:numel(Lvc)
    xv = x(Lvc{k});
    yv = y(Lvc{k});
    B(:,k) = fminsearch(@(b)fcn(b,xv,yv), rand(3,1));                   % Estimate Parameters
end

figure
plot(x, y)
hold on
for k = 1:size(B,2)
    plot(B(1,k),B(2,k),'xr')
    plot(x(Lvc{k}([1,end])), y(Lvc{k}([1,end])), 'r+')
    text(B(1,k),B(2,k), sprintf('R = %.3f', abs(B(3,k))), 'VerticalAlignment','bottom', 'HorizontalAlignment','center')
end
hold off
grid

The find calls are necessary in order to identify the ends of the arcs. Otherwise, logical indexing alone would suffice in my code. The ends of the arcs are defined as:

xend = x(Lvc{k}([1,end]));
yend = y(Lvc{k}([1,end]));

for each arc, with ‘k’ defined as a for loop index.

My code then produces this figure:

The red ‘x’ markers are the centres of each arc, the red ‘+’ markers are the ends of the arcs (as my code identifies them).

Experiment to get the result you want.

[Conrade Muyambo]

Thank you very much. It is working.. But if I want to know the start and ending poins of the arcs using this xend = x(Lvc{k}([1,end])); I get two points. Am I not supposed to get 8 values also for the yend = y(Lvc{k}([1,end]));

[Star Strider]

My pleasure.

You get the beginning and end indices by addressing them in a for loop:

for k = 1:size(B,2)
    arcsx(:,k) = x(Lvc{k}([1,end]));
    arcsy(:,k) = y(Lvc{k}([1,end]));
end

Each column of each matrix gives the indices of the starting and ending points for the arc.

For example, the matrices are:

arcsx =
    3.6250    3.3500    0.7833    1.0833
    3.3250    3.6417    1.1167    0.7167

and:

arcsy =
    2.9500    0.7000    0.9417    3.3333
    3.3417    0.9750    0.5750    2.9583

so the start and end points for ‘NE’ are:

       x         y
    3.6250    2.9500
    3.3250    3.3417 

for ‘SE’:

       x         y
    3.3500    0.7000
    3.6417    0.9750 

and so for the rest of the arcs.

The columns are in clockwise order of the corners: ‘NE’, ‘SE’, ‘SW’, ‘NW’, in columns 1 through 4 respectively. There is no straightforward way of presenting these so it is also easy to recover them and work with them.

If my Answer helped you solve your problem, please Accept it!

[Conrade Muyambo]

Thank you. That helped a lot. But is it now possible, with the start point of one arc and the end point of the second arc to fit a curve to describe the side?

[Star Strider]

As always, my pleasure.

No.

The arc parameters describe the arcs only.

I have no idea how to describe the sides, since that is not what you requested initially.

See my oiriginal Answer (the ‘Create Data’ assignments) to understand how to describe the arcs.

Answer 2

[Conrade Muyambo]

It served the purpose sir! That was only a question. Thank you

[Star Strider]

As always, my pleasure!

You can likely adapt my code to estimate the sides that are not part of the arcs, although since they are more linear that nonlinear, a linear or ‘slightly nonlinear’ regression could work.

Answer 3

[Conrade Muyambo]

Okay I understand, i will work on it. I will post a different topic later, but not related to this regarding the same shape,.