CODEWITHSUNDEEP
matlab
Aidan K
Aidan K

Reputation: 55

Calculate heading angle from x and y information

I have data that records the x and y positions of an animal in a 2D assay over time stored in a matlab matrix. I can plot these co-ordinates over time, and extract the velocity information and plot this using cline.

The problem I am having at the moment is calculating the heading angle. It should be a trivial trigonometry question, but I am drawing a blank on the best way to start.

The data is stored in a matrix xy representing x and y co-ordinates:

    796.995391705069    151.755760368664
    794.490825688073    150.036697247706
    788.098591549296    145.854460093897
    786.617021276596    144.327659574468
    781.125000000000    140.093750000000
    779.297872340426    138.072340425532
    775.294642857143    133.879464285714

What I would like to be able to do is know the angle of the line drawn from (796.995, 151.755) to (794.490, 150.036), and so on. My research suggests atan2 will be the appropriate function, but I am unsure how to call it correctly to give useful information.

    difx = xy(1,1) - xy(2,1);
    dify = xy(1,2) - xy(2,2);
    angle = atan2(dify,difx);
    angle = angle*180/pi % convert to degrees

The result is 34.4646. Is this correct?

If it is correct, how do I get the value to be in the range 0-360?

Upvotes: 2

Views: 19304

Answers (2)

Floris
Floris

Reputation: 46415

You can use the diff function to get all the differences at once:

dxy = diff(xy); % will contain [xy(2,1)-xy(1,1) xy(2,2)-xy(1,2); ...

Then you compute the angle using the atan2 function:

a = atan2(dxy(:,2), dxy(:,1));

You convert to degrees with

aDeg = 180 * a / pi;

And finally take the angle modulo 360 to get it between 0 and 360:

aDeg = mod(aDeg, 360);

So - you pretty much got it right, yes. Except that you have calculated the heading from point 2 to point 1, and I suspect you want to start at 1 and move towards 2. That would give you a negative number - or modulo 360, an angle of about 325 degrees.

Also, using the diff function gets you the entire array of headings all at once which is a slight improvement over your code. [rc mi]= EDIT the problem of "phase wrapping" - when the heading goes from 359 to 0 - is quite a common problem. If you are interested in knowing when a large change happens, you can try the following trick (using aDeg from above - angle in degrees).

dDeg1 = diff(aDeg);  % the change in angle
dDeg2 = diff(mod(aDeg + 90, 360)); % we moved the phase wrap point by 180 degrees
dDeg12 = [dDeg1(:) dDeg2(:)]';
[rc mi]= min(abs(dDeg12));
indx = sub2ind(size(dDeg12), mi, 1:size(dDeg12, 2));
result = dDeg12(ii);

What I did there: one of the variables (dDeg or dDeg2) does not see the phase wrap, and the min function finds out which one (it will have a smaller absolute difference). The sub2ind looks up that number (it is either positive or negative - but it's the smaller one of the two), and that is the value that ends up in result.

Upvotes: 3

Steve
Steve

Reputation: 4097

You can verify the angle by plotting a little line that starts at the first point and end in the direction of the heading. If the angle is correct, it will point in the direction of the next point in xy. Everything depends on where yo define 0 degrees at (straight up, say) from and whether positive degrees is rotation counterclockwise (I do) or clockwise. In MATLAB you can get the numbers between 0 and 360 but using modulo---or you can just add 180 to your results but this will change the definition of where the 0 degree mark is.

I made the following script that is a bit complex but shows how to calculate the heading/angle for all points in vector format and then displays them.

xy =[ 796.995391705069    151.755760368664
    794.490825688073    150.036697247706
    788.098591549296    145.854460093897
    786.617021276596    144.327659574468
    781.125000000000    140.093750000000
    779.297872340426    138.072340425532
    775.294642857143    133.879464285714];
% t = linspace(0,3/2*pi, 14)';
% xy = [sin(t), cos(t)];

% calculate the angle:
myDiff = diff(xy);
myAngle = mod(atan2(myDiff(:,1), myDiff(:,2))*180/pi, 360);

% Plot the original Data:
figure(1);
clf;
subplot(1,3,1);
plot(xy(:,1), xy(:,2), '-bx', 'markersize', 12);
hold all
axis equal;grid on;
title('Original Data');

% Plot the calculated angle:
subplot(1,3,2);
plot(myAngle);
axis tight; grid on;
title('Heading');

% Now plot the result with little lines pointing int he heading:
subplot(1,3,3);
plot(xy(:,1), xy(:,2), '-bx', 'markersize', 12);
hold all

% Just for visualization:
vectorLength = max(.8, norm(xy(1,:)- xy(2,:)));

for ind = 1:length(xy)-1
    startPoint = xy(ind,:)';
    endPoint = startPoint + vectorLength*[sind(myAngle(ind)); cosd(myAngle(ind))];
    myLine = [startPoint, endPoint];
    plot(myLine(1,:), myLine(2, :), ':r ', 'linewidth', 2)
end

axis equal;grid on;
title('Original Data with Heading Drawn On');

For example, if you use my test data

t = linspace(0,3/2*pi, 14)';
xy = [sin(t), cos(t)];

You get the following:

example plot 1

and if you do yours you get

example plot 2

Note how the little red line starts at the original data point and moves in the direction of the next point---just like the original blue line connecting the points.

Also note that the use of diff in the code to difference all the points properly at once. This is faster and avoids any problems with the direction--looks like in your case it's swapped.

Upvotes: 1

Related Questions