有几件事:
如果你有一个线性函数 y=mx+b,那条线的角度是atan(m)
,而不是m
。这些对于小m', but very different for large
m` 大致相同。
2+ 阶 polyfit 的线性分量不同于 1 阶 polyfit 的线性分量。您需要对数据进行两次拟合,一次在您的工作级别,一次使用一阶拟合。
给定一个斜率m
,有比使用三角函数更好的计算旋转矩阵的方法(例如cos(atan(m))
)。在执行几何时,我总是尽量避免使用三角函数,并用线性代数运算代替它们。这通常更快,并且导致奇点问题更少。请参阅下面的代码。
这种方法会导致一些轨迹出现问题。例如,考虑一个北/南轨迹。但这是一个更长的讨论。
使用所描述的方法,加上上面的注释,这里是一些实现这个的示例代码:
%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')