how to find the sharp turn in a 2d-line (curve)?

Hi all,

I am trying to find the sharp turn in a 2d-line (curve). Line is constructed with two vectors, X and Y. In following link you can find a sample line with realized point at which there is sharp turn (red solid point).

https://dl.dropbox.com/u/2190947/line.pdf

I appreciate if you could help me out with this.

Thanks, Payam

Answer 1

[Roger Stafford]

One measure of a "sharp turn" is the amount of curvature between three adjacent points on your curve. Let (x1,y1), (x2,y2), and (x3,y3) be three such adjacent points. By a well-known formula, the curvature of a circle drawn through them is simply four times the area of the triangle formed by the three points divided by the product of its three sides. Using the coordinates of the points this is given by:

 K = 2*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1)) ./ ...
   sqrt(((x2-x1).^2+(y2-y1).^2)*((x3-x1).^2+(y3-y1).^2)*((x3-x2).^2+(y3-y2).^2));

Roger Stafford

[Image Analyst]

Or maybe not so well known, since this is asked here fairly often. I've added Roger's clever solution to the FAQ. Here's a demo based on Roger's formula:

% This demo plots a soft star and then uses ROger Stafford's formula
% to find and mark the locations on the star
% that have high curvature.

% Initialization & clean up stuff.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format longg;
format compact;
fontSize = 20;

%=====================================================
% First make a shape with sharp turns or cusps.
% Demo macro to draw a rounded star (like a splat).

% Select the inner and outer radius.
outerRadius = 44  % You can change this
innerRadius = 19  % You can change this
% Select the number of lobes around the circle.
numberOfLobes = 8;  % You can change this

period = 2 * pi / numberOfLobes;
meanRadius = (outerRadius + innerRadius)/2
amplitude = (outerRadius - innerRadius)/2
t = (0:.005:1)*2*pi; % Independent parameter.
variableRadius = amplitude * cos(2*pi*t/period) + meanRadius;
subplot(2,2,1);
plot(variableRadius, 'LineWidth', 2);
grid on;
ylim([0 outerRadius]);
title('VariableRadius', 'FontSize', fontSize);
period = 2*pi;  % Need to change this now.
xStar = variableRadius .* cos(2*pi*t/period); 
yStar = variableRadius .* sin(2*pi*t/period); 
subplot(2,2,2);
plot(t, xStar, 'LineWidth', 2);
grid on;
title('x2 vs. t', 'FontSize', fontSize);
subplot(2,2,3);
plot(t, yStar, 'LineWidth', 2);
grid on;
title('y2 vs. t', 'FontSize', fontSize);
subplot(2,2,4);
plot(xStar, yStar,'b', 'LineWidth', 2) 
title('x2 vs y2', 'FontSize', fontSize);
axis square;
% Maximize window.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]); % Maximize figure.
set(gcf,'name','Image Analysis Demo','numbertitle','off') 
% OK - all of the above code was just to get some demo data
% that we can use to find the high radius of curvature locations on.

%=====================================================
% Now run along the (x2, y2) soft star curve 
% and find the radius of curvature at each location.
numberOfPoints = length(xStar);
curvature = zeros(1, numberOfPoints);
for t = 1 : numberOfPoints
  if t == 1
    index1 = numberOfPoints;
    index2 = t;
    index3 = t + 1;
  elseif t >= numberOfPoints
    index1 = t-1;
    index2 = t;
    index3 = 1;
  else
    index1 = t-1;
    index2 = t;
    index3 = t + 1;
  end
  % Get the 3 points.
  x1 = xStar(index1);
  y1 = yStar(index1);
  x2 = xStar(index2);
  y2 = yStar(index2);
  x3 = xStar(index3);
  y3 = yStar(index3);
  % Now call Roger's formula:
  % http://www.mathworks.com/matlabcentral/answers/57194#answer_69185
  curvature(t) = 2*abs((x2-x1).*(y3-y1)-(x3-x1).*(y2-y1)) ./ ...
  sqrt(((x2-x1).^2+(y2-y1).^2)*((x3-x1).^2+(y3-y1).^2)*((x3-x2).^2+(y3-y2).^2));
end
% Plot curvature.
figure;
subplot(2, 1, 1);
plot(curvature, 'b-', 'LineWidth', 2);
grid on;
xlim([1 numberOfPoints]); % Set limits for the x axis.
title('Radius of Curvature', 'FontSize', fontSize);

% Find high curvature points - 
% indexes where the curvature is greater than 0.3
highCurvatureIndexes = find(curvature > 0.3);

% Plot soft star again so we can plot 
% high curvature points over it.
subplot(2, 1, 2);
plot(xStar, yStar,'b', 'LineWidth', 2) 
grid on;
axis square;

% Mark high curvature points on the star.
hold on;
plot(xStar(highCurvatureIndexes), yStar(highCurvatureIndexes), ...
  'rd', 'MarkerSize', 15, 'LineWidth', 2);
title('Soft Star with High Curvature Locations Marked', 'FontSize', fontSize);

% Maximize window.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]); % Maximize figure.
set(gcf,'name','Image Analysis Demo','numbertitle','off');
helpdlg('Done with demo');

[Maradona Rodrigues]

Hi I tried using the above code and resulted in some erronous results

my 2d line is list = [0 0; 4 0.5; 8 6; 6 25; 3 7; 1 1] , with the biggest curvature being at the point [6,25]. However i didnt get the second lowest at that point?

[Image Analyst]

Maybe there are not enough points to make a determination. If you have a recent version of MATLAB you might also try findchangepts().

Answer 2

[Jan]

"Sharp" is relative. There is always a zoom level, which let a curve look smooth.

If you do not have a curve defined by a function, but a piecewise defined line, you are looking for neighboring elements with and included angle above a certain limit. But when such a piece has a length of 1e-200, while the others have a length of 1.0, can this have a "sharp turn"?!

But let's imagine, that you can control this fundamental problem by inventing some meaningful thresholds. Then this determines the angle between two lines:

angle = atan2(norm(cross(N1, N2)), dot(N1, N2))

Answer 3

[Image Analyst]

Well for that example, just do

yAtTurn = min(y);
xAtTurn = find(y == yAtTurn);

If you need something more general, flexible, and robust, then you need to say how other curves might look different than the one example you supplied.

[tafteh]

Thanks for your quick response. I updated the sample image as noted in the following post.

[Image Analyst]

What is your data? Is it an image or a set of (x,y) coordinates?

[tafteh]

Its a set of (x,y) coordinates.

[Image Analyst]

Perfect - then my code will work just fine. I put it under Roger's answer since it was really his idea for the algorithm - I just packaged it into a nice demo. Go ahead and mark it as the official "Answer".

[Alessandro]

Hi Image Analyst... your code works fine for me but i have a bw image instead of a set of (x,y) points... so that, i've not an ordered set of (x,y)... how can i figure out? (my purpose is still find inflection points)

[Image Analyst]

You can use bwboundaries() to get a list of (x,y) points.

Roger's Answer is in the FAQ http://matlab.wikia.com/wiki/FAQ#How_do_I_find_.22kinks.22_in_a_curve.3F

Answer 4

[tafteh]

Hi again, I have updated the image to make it more clearer regarding my problem explained in the first post. Please find the link below:

https://dl.dropbox.com/u/2190947/line2.pdf

Thanks,