Rotating Axes Around Line of Fit MATLAB

Question

I'm currently frustrated by the following problem:

I've got trajectory data (i.e.: Longitude and Latitude data) which I interpolate to find a linear fitting (using polyfit and polyval in matlab).

What I'd like to do is to rotate the axes in a way that the x-axis (the Longitude one) ends up lying on the best-fit line, and therefore my data should now lie on this (rotated) axis.

What I've tried is to evaluate the rotation matrix from the slope of the line-of-fit (m in the formula for a first grade polynomial y=mx+q) as

[cos(m) -sin(m);sin(m) cos(m)]

and then multiply my original data by this matrix...to no avail!

I keep obtaining a plot where my data lay in the middle and not on the x-axis where I expect them to be.

What am I missing?

Thank you for any help!

Best Regards,

Wintermute

Answer 1

A couple of things:

1. If you have a linear function y=mx+b, the angle of that line is atan(m), not m. These are approximately the same for small m', but very different for largem`.

2. The linear component of a 2+ order polyfit is different than the linear component of a 1st order polyfit. You'll need to fit the data twice, once at your working level, and once with a first order fit.

3. Given a slope m, there are better ways of computing the rotation matrix than using trig functions (e.g. cos(atan(m))). I always try to avoid trig functions when performing geometry and replace them with linear algebra operations. This is usually faster, and leads to fewer problems with singularities. See code below.

4. This method is going to lead to problems for some trajectories. For example, consider a north/south trajectory. But that is a longer discussion.

Using the method described, plus the notes above, here is some sample code which implements this:

%Setup some sample data
long = linspace(1.12020, 1.2023, 1000);
lat  = sin ( (long-min(long)) / (max(long)-min(long))*2*pi  )*0.0001 + linspace(.2, .31, 1000);

%Perform polynomial fit
p = polyfit(long, lat, 4);

%Perform linear fit to identify rotation
pLinear = polyfit(long, lat, 1);
m = pLinear(1);  %Assign a common variable for slope
angle = atan(m);

%Setup and apply rotation
%    Compute rotation metrix using trig functions
rotationMatrix = [cos(angle) sin(angle); -sin(angle) cos(angle)];

%    Compute same rotation metrix without trig
a = sqrt(m^2/(1+m^2));  %a, b are the solution to the system:    
b = sqrt(1/(1+m^2));    %    {a^2+b^2 = 1}, {m=a/b}
%                       %That is, the point (b,a) is on the unit
%                       circle, on a line with slope m
rotationMatrix = [b a; -a b];  %This matrix rotates the point (b,a) to (1,0)

%    Generally you rotate data after removing the mean value
longLatRotated = rotationMatrix * [long(:)-mean(long) lat(:)-mean(lat)]';

%Plot to confirm
figure(2937623)
clf

subplot(211)
hold on
plot(long, lat, '.')
plot(long, polyval(p, long), 'k-')
axis tight
title('Initial data')
xlabel('Longitude')
ylabel('Latitude')

subplot(212)
hold on;
plot(longLatRotated(1,:), longLatRotated(2,:),'.b-');
axis tight
title('Rotated data')
xlabel('Rotated x axis')
ylabel('Rotated y axis')

Answer 2

The angle you are looking for in the rotation matrix is the angle of the line makes to the horizontal. This can be found as the arc-tangent of the slope since:

tan(\theta) = Opposite/Adjacent = Rise/Run = slope

so t = atan(m) and noting that you want to rotate the line back to horizontal, define the rotation matrix as:

R = [cos(-t) sin(-t)
     sin(-t) cos(-t)]

Now you can rotate your points with R