2

我有数据记录了动物在 2D 测定中随时间存储在 matlab 矩阵中的 x 和 y 位置。我可以随时间绘制这些坐标,并提取速度信息并使用 cline 绘制它。

我现在遇到的问题是计算航向角。这应该是一个微不足道的三角问题,但我在最好的开始方式上画了一个空白。

数据存储在表示 x 和 y 坐标的矩阵 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

我想做的是知道从 (796.995, 151.755) 到 (794.490, 150.036) 绘制的线的角度,等等。我的研究表明 atan2 将是合适的函数,但我不确定如何正确调用它以提供有用的信息。

    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

结果是 34.4646。这个对吗?

如果正确,如何使值在 0-360 范围内?

4

2 回答 2

3

您可以使用该diff函数一次获取所有差异:

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

然后使用以下atan2函数计算角度:

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

你转换成度数

aDeg = 180 * a / pi;

最后取角度模 360 得到 0 到 360 之间:

aDeg = mod(aDeg, 360);

所以 - 你几乎是对的,是的。除了你已经计算了从点 2 到点 1 的航向,我怀疑你想从 1 开始向 2 移动。这会给你一个负数 - 或模 360,一个大约 325 度的角度。

此外,使用该diff函数可以一次获得整个标题数组,这对您的代码略有改进。[rc mi]= 编辑“相位环绕”的问题 - 当标题从 359 变为 0 时 - 是一个很常见的问题。如果您想知道何时发生大的变化,您可以尝试以下技巧(aDeg从上方使用 - 以度为单位的角度)。

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);

我在那里做了什么:其中一个变量(dDegdDeg2)没有看到相位环绕,并且该min函数找出哪个变量(它将具有较小的绝对差)。sub2ind查找该数字(它是正数或负数 - 但它是两者中较小的一个),这就是最终在result.

于 2013-09-27T21:21:31.170 回答
1

您可以通过绘制一条从第一个点开始并在航向方向结束的小线来验证角度。如果角度正确,它将指向 xy 中下一个点的方向。一切都取决于你在哪里定义 0 度(比如说),以及正度是逆时针旋转(我愿意)还是顺时针旋转。在 MATLAB 中,您可以获得 0 到 360 之间的数字,但使用模数 - 或者您可以将 180 添加到结果中,但这会改变 0 度标记所在位置的定义。

我制作了以下有点复杂的脚本,但显示了如何计算矢量格式的所有点的航向/角度,然后显示它们。

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');

例如,如果您使用我的测试数据

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

您会得到以下信息:

示例图 1

如果你做你的,你会得到

示例图 2

注意小红线是如何从原始数据点开始向下一个点移动的——就像连接这些点的原始蓝线一样。

另请注意,diff在代码中使用 可以一次正确地区分所有点。这更快并且避免了方向上的任何问题——在你的情况下看起来它被交换了。

于 2013-09-27T21:37:05.813 回答