围绕拟合线旋转轴 MATLAB答案

【问题标题】：Rotating Axes Around Line of Fit MATLAB围绕拟合线旋转轴 MATLAB
【发布时间】：2013-03-11 17:56:53
【问题描述】：

我目前对以下问题感到沮丧：

我有轨迹数据（即：经度和纬度数据），我对其进行插值以找到线性拟合（在 matlab 中使用 polyfit 和 polyval）。

我想要做的是以 x 轴（经度轴）最终位于最佳拟合线上的方式旋转轴，因此我的数据现在应该位于此（旋转）轴。

我尝试从拟合线的斜率（一级多项式y=mx+q 的公式中的 m）评估旋转矩阵

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

然后将我的原始数据乘以这个矩阵......无济于事！

我不断获得一个图，其中我的数据位于中间，而不是在我期望的 x 轴上。

我错过了什么？

感谢您的帮助！

最好的问候，

冬眠

【问题讨论】：

使用这个旋转矩阵，你已经假设旋转轴在原点。这是真的还是您还需要翻译数据？另外，我无法从您写的内容中理解出了什么问题，您能否提供一些示例数据\您使用的代码？
实际上，正如 Azim 的 Pursuit 让我注意到的那样，我应该使用的矩阵必须用 atan(m)（或与之等效的线性代数）而不是 m 来构建，并且应该转置获得“别名”类型的旋转。现在你让我注意到了，这是真的：在旋转问题中，我忽略了y=mx+q 等式中的+q...我还应该添加翻译吗？如何？感谢您的帮助！

标签： matlab rotational-matrices cartesian-coordinates

【解决方案1】：

有几点：

如果你有一个线性函数 y=mx+b，那条线的角度是atan(m)，而不是m。这些对于小的 m', but very different for largem` 大致相同。
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')

【讨论】：

谢谢！真的很有帮助！
我终于去了办公室并尝试了代码：我知道你在那里做了什么，这是我的错误！您的线性矩阵等效于将 m 三角解释为角度的正切（因此 cos(angle)=(1/sqrt(1+tan^2(angle))) 等）并且您使用了所谓的“别名”旋转（轴围绕向量的旋转，它是规范的“不在场证明”旋转矩阵的转置......我错误地使用的那个）：我解释正确了吗？
我将旋转矩阵的两个版本称为“body”旋转（在框架中旋转物体）和“frame”旋转（旋转框架）。老实说，我不确定这是哪一个；我已经这样做了多年，我总是需要检查结果。通常实际做的是选择一个点（在这种情况下，（b，a）），决定我想要那个点去哪里（例如（b，a）->（1,0）），写下（希望如此）正确的矩阵，然后检查它。我发现的最有效的方法是计划编写大量测试代码。

【解决方案2】：

您在旋转矩阵中寻找的角度是线与水平线的角度。这可以作为斜率的反正切得到，因为：

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

所以t = atan(m) 并注意到您要将线旋转回水平，将旋转矩阵定义为：

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

现在您可以使用R 旋转您的点

【讨论】：