Android:SensorManager.getRotationMatrix 和 SensorManager.getOrientation() 的算法

Android: Algorithms for SensorManager.getRotationMatrix and SensorManager.getOrientation()

要在 Android 中获取欧拉角(例如俯仰角、横滚角、方位角)的方位,需要执行以下操作:

  1. SensorManager.getRotationMatrix(float[] R, float[] I, float[] 重力, float[] 地磁);
  2. SensorManager.getOrientation(float[] R, float[] 方向);

在第一个中,我意识到它使用了一种TRIAD算法;旋转矩阵(R[])由重力组成,地磁X重力,重力X(地磁X重力)---X是叉积。
请参阅以下代码:

    float Ax = gravity[0];
    float Ay = gravity[1];
    float Az = gravity[2];
    final float Ex = geomagnetic[0];
    final float Ey = geomagnetic[1];
    final float Ez = geomagnetic[2];
    float Hx = Ey*Az - Ez*Ay;
    float Hy = Ez*Ax - Ex*Az;
    float Hz = Ex*Ay - Ey*Ax;
    final float normH = (float)Math.sqrt(Hx*Hx + Hy*Hy + Hz*Hz);
    if (normH < 0.1f) {
        // device is close to free fall (or in space?), or close to
        // magnetic north pole. Typical values are  > 100.
        return false;
    }
    final float invH = 1.0f / normH;
    Hx *= invH;
    Hy *= invH;
    Hz *= invH;
    final float invA = 1.0f / (float)Math.sqrt(Ax*Ax + Ay*Ay + Az*Az);
    Ax *= invA;
    Ay *= invA;
    Az *= invA;
    final float Mx = Ay*Hz - Az*Hy;
    final float My = Az*Hx - Ax*Hz;
    final float Mz = Ax*Hy - Ay*Hx;
    if (R != null) {
        if (R.length == 9) {
            R[0] = Hx;     R[1] = Hy;     R[2] = Hz;
            R[3] = Mx;     R[4] = My;     R[5] = Mz;
            R[6] = Ax;     R[7] = Ay;     R[8] = Az;
        } else if (R.length == 16) {
            R[0]  = Hx;    R[1]  = Hy;    R[2]  = Hz;   R[3]  = 0;
            R[4]  = Mx;    R[5]  = My;    R[6]  = Mz;   R[7]  = 0;
            R[8]  = Ax;    R[9]  = Ay;    R[10] = Az;   R[11] = 0;
            R[12] = 0;     R[13] = 0;     R[14] = 0;    R[15] = 1;
        }
    }

但是,我无法理解SensorManager.getOrientation()。

 azimuth = (float)Math.atan2(R[1], R[4]);
 pitch = (float)Math.asin(-R[7]);
 roll = (float)Math.atan2(-R[6], R[8]);

获得欧拉角的确切算法是什么?

让我试着解释一下: getRotationMatrix 根据重力和磁力矢量组成旋转矩阵。

我们这里的主要目标是构建NED frame

我们假设重力指向地球中心,磁铁指向北极。但在实际情况下,这些向量是不垂直的,这就是为什么我们首先计算与 E 和 A 正交且属于切平面的向量 H 的原因。 H 是叉积 (E x A) 并且与 E 和 A 正交。

float Hx = Ey*Az - Ez*Ay;
float Hy = Ez*Ax - Ex*Az;
float Hz = Ex*Ay - Ey*Ax;
final float normH = (float)Math.sqrt(Hx*Hx + Hy*Hy + Hz*Hz);

归一化加速度和H向量(因为这些向量将构成ENU坐标系的基础)

final float invH = 1.0f / normH;
    Hx *= invH;
    Hy *= invH;
    Hz *= invH;
final float invA = 1.0f / (float)Math.sqrt(Ax*Ax + Ay*Ay + Az*Az);
    Ax *= invA;
    Ay *= invA;
    Az *= invA;

找到最后一个基向量 (M) 作为 H 和 A 的叉积:

double Mx = Ay * Hz - Az * Hy;
double My = Az * Hx - Ax * Hz;
double Mz = Ax * Hy - Ay * Hx;

体坐标系中任意向量(a)的坐标通过NED坐标表示为a = Ra' R - Transformation matrix矩阵,其列为新基向量在old basis中的坐标

但是NED坐标系中的坐标是 计算为 a' = T^(-1) * a。对于正交变换矩阵逆等于转置矩阵。因此我们有:

R[0] = Hx;     R[1] = Hy;     R[2] = Hz;
R[3] = Mx;     R[4] = My;     R[5] = Mz;
R[6] = Ax;     R[7] = Ay;     R[8] = Az;

一旦我们有了旋转矩阵,我们就可以将其转换为欧拉角表示。转换公式取决于您使用的 convention。 你的公式

azimuth = (float)Math.atan2(R[1], R[4]);
pitch = (float)Math.asin(-R[7]);
roll = (float)Math.atan2(-R[6], R[8]);

对于具有约定 Y-X-Z 的 Tiat Bryan 角度是正确的。为了更好地理解从旋转矩阵到欧拉角的转换,我建议研究 Gregory G. Slabaugh - "Computing Euler angles from a rotation matrix"

的一篇文章

我将给出 getOrientation 的几何解释,并解释如何仅使用 RotationMatrix 计算各种旋转。如果您了解 getOrientation 的作用,那么就没有必要使用它,如果您不了解,那么使用它会给您带来各种麻烦。

具体来说,它将回答许多关于在 Stack Exchange 上 posted 的问题,例如

  • 在人像模式下为平板电脑编写指南针。为什么调用 getRotation 后需要加上 90 度才能得到正确答案。
  • 为什么要先调用RemapCoordinateSystem再调用getOrientation获取后置摄像头的方向(设备的负方向z-axis)
  • 为什么在调用RemapCoordinateSystem然后getOrientation正确获取设备的方向后z-axis,当设备是平的时候结果就没有意义了。
  • 如何独立于设备位置计算围绕任何设备轴的旋转,即平坦或不平坦。

首先我需要更详细地解释一下RotationMatrix以及它能做什么,然后再解释getOrientation

这里要考虑 2 个坐标系。
一种是世界坐标系,x-axis 指向东方,y-axis 指向北方,z-axis 指向天空。
另一个是设备坐标系,x-axis 短边在 phone 上(平板电脑上的长边),y-axis 长边在 phone 上(短边)在平板电脑上)和 z-axis 是正交于指向您的屏幕的向量。

从数学上讲,我们的对象是 3 维实向量 space,我们为此向量使用以下基数 space。

世界基础W = {E, N, SKY} 其中
E是位于东方向的单位向量
N是位于北方向的单位向量
SKY是天空方向的单位向量

设备基础D = {X, Y, Z} 其中
X是位于phone短边方向的单位向量(tablet为长边)
Y是位于phone长边方向的单位向量(平板电脑为短边)
Z 是一个单位向量,位于与指向您的屏幕正交的方向上

对于那些忘记什么是基的人来说,这意味着3 space中的每个向量v都可以写成

以世界为基础
v = a_1 E + a_2 N + a_3 SKY 其中 a_1, a_2, a_3 是实数
通常写为 (a_1, a_2, a_3)_W 或只是 (a_1, a_2, a_3)。这 3 个元组被称为 v 相对于世界基础的坐标,或者只是 v 的坐标,隐含的基础被理解为世界基础

在设备基础上
v = b_1 X + b_2 Y + b_3 Z 其中 b_1, b_2, b_3 是实数
通常写成 (b_1, b_2, b_3)_W 或只是 (b_1, b_2, b_3)。这 3 个元组称为 v 相对于 Device 基础的坐标,或者只是 v 的坐标,其基础隐式理解为 Device 基础。传感器返回的值以设备为基础,例如加速度计返回为值 [0]、值 [1] 和值 [2],写在设备基础中将是
acc = 值[0] X + 值[1] Y + 值[2] Z

特别是
X = 1 X + 0 Y + 0 Z
Y = 0 X + 1 Y + 0 Z
Z = 0 X + 0 Y + 1 Z

XYZ相对于设备的坐标basis 分别是 (1, 0, 0), (0, 1, 0) 和 (0, 0, 1)。您将在稍后的 getOrientation 中看到如何使用这些向量来解释计算。

请注意,World basis 是固定的,但 Device basis 会随着 phone 位置的变化而变化。即单位向量X,Y,Z随着设备位置的变化而变化。它们以相同的方式定义,但是当位置改变时它们是不同的向量。你仍然把它们写成 X, Y, Z 但它们是不同的 X, Y, Z.
因此,如果 phone 静止不动,作用在 phone 上的唯一力是重力,因此加速度计矢量理论上是位于世界天空轴上的矢量,即在世界基础上坐标为 (0, 0,克)。在设备基础上,对于某些 a_1、a_2 和 a_3,它是 (a_1、a_2、a_3)。如果设备静止在不同的方向,加速度计仍然是相同的,即世界基础上的 (0, 0, g) 但现在它是 (b_1,b_2、b_3) 在设备基础中,其中至少一个 a 与 b 不同。

getRotationMatrix 的参数包括 gravitygeomagneticgravity 假设是真实的 gravity 即它在世界基础上的坐标是 (0, 0, k) 并且地磁假设位于世界 N-Sky 平面上它的坐标是(0, a, b) 以世界为基础。

根据这些假设,让我们看看 Rotation Matrix 是如何在 getRotationMatrix 中计算的。它在这个方法中尝试的是获得一个矩阵 M 使得给定任何在 Device 中具有坐标 (a_1, a_2, a_3) 的向量根据产品 M (a_1, a_2, a_3)_T(转置)给出坐标 (b_1, b_2, b_3) 以世界为基础。从数学上讲,getRotationMatrix 计算基矩阵的变化。

gravityg的传入参数应该是onSensorChangedTYPE_GRAVITY或[=的低通滤波器得到的值50=] 和 geomagnetic m 应该是 TYPE_MAGNETIC_FIELD

的值

g = a_1 X + a_2 Y + a_3 Z
m = b_1 X + b_2 Y + b_3 Z

不失一般性假设gm已经归一化,即g[=的范数385=] 和 m 等于 1。 我们现在要做的是在Device basis中写上World basis{E, N, SKY}

由于g假设为重力,即世界基础中的(0, 0, 1) 这正是向量SKY.

天空 = g = a_1 X + a_2 Y + a_3 Z
SKY在Device basis中的坐标为(a_1,a_2,a_3)

现在因为 m 被假定位于世界 N-SKY 平面 mg 是与世界 N-SKY 平面正交的单位向量,指向右侧,因此这是 E。

E = m x g = b_1 X + b_2 Y + b_3 Z

也就是E相对于Device basis的坐标是(b_1, b_2, b_3) b的是mg

的叉积得到的值

最后SKYE的叉积是一个正交于E-SKY平面的向量,是N

N = 天空 x E = c_1 X + c_2 Y + c_3 Z

我们有

E = b_1 X + b_2 Y + b_3 Z
N = c_1 X + c_2 Y + c_3 Z
天空 = a_1 X + a_2 Y + a_3 Z

因此给定World basis中的任何坐标我们都可以找到Device basis中的坐标。例如 (1, 2, 3) 是一个向量 v

v = E + 2 N + 3天空

写在设备坐标中将通过替换并乘以

v = b_1 X + b_2 Y + b_3 Z + 2 (c_1 X + c_2 Y + c_3 Z) + 3 (a_1 X + a_2 Y + a_3 Z)
v = (b_1 + 2 c_1 + 3 a_1) X + (b_2 + 2 c_2 + 3 a_2) Y + (b_3 + 2 c_3 + 3 a_3) Z

但我们真正感兴趣的是在设备基础上给定任何坐标的逆向找到世界基础上的坐标。好吧,对于那些不记得线性代数的人来说,在上面我们有 3 个未知的 3 个方程,因此我们应该能够求解这个方程 XYZ 根据 ENSKY

从数学上讲,将上面的a,b,c分列得到的矩阵就是基矩阵从World基到Device基的变化,所以我们需要求它的逆矩阵。但是这个矩阵是一个正交矩阵,因此它的逆矩阵就是它的转置。

写在代码里

a[0]  a[1]  a[2]
a[3]  a[4]  a[5]
a[6]  a[7]  a[8]

给定设备基础中的任何坐标(a_1、a_2、a_3)我们可以找到坐标(b_1、b_2、b_3) 在世界坐标系中,取上面矩阵的乘积和 (a_1, a_2, a_3) 的转置。

特别是

a[0]  a[1]  a[2]      1       a[0]
a[3]  a[4]  a[5]  x   0   =   a[3]
a[6]  a[7]  a[8]      0       a[6]

因此的坐标 在世界基础上是 (a[0], a[3], a[6])

a[0]  a[1]  a[2]      0       a[1]
a[3]  a[4]  a[5]  x   1   =   a[4]
a[6]  a[7]  a[8]      0       a[7]

因此Y在World basis中的坐标为(a[1], a[4], a[7])

a[0]  a[1]  a[2]      0       a[2]
a[3]  a[4]  a[5]  x   0   =   a[5]
a[6]  a[7]  a[8]      1       a[8]

因此Z在World basis中的坐标为(a[2],a[5],a[8])

现在让我们看看 getOrientation 做了什么。

它的代码是

azimuth = (float)Math.atan2(R[1], R[4]);
pitch = (float)Math.asin(-R[7]);
roll = (float)Math.atan2(-R[6], R[8]);

文档说

values[0]: azimuth, rotation around the -Z axis, i.e. the opposite direction of Z axis.
values[1]: pitch, rotation around the -X axis, i.e the opposite direction of X axis.
values[2]: roll, rotation around the Y axis.

不考虑设备的位置就照字面理解上面的文档是错误的。仅当设备是扁平的时,文档才适用于 azimuth

在我们继续之前,请注意所有计算都必须使用相对于同一基础的坐标来完成。例如,如果你想找到 2 个向量之间的角度,这 2 个向量的坐标必须相对于相同的基础。

现在当你计算旋转时,你有点隐含地理解旋转是从固定位置的偏差。也就是说,如果答案是 90 度,那么它与什么成 90 度?在 azimuth 的情况下,它与磁北成 90 度,即您计算与 N 矢量的偏差。如果你不能拼出这个隐含的固定位置,你就会遇到很多麻烦。例如,音高的固定矢量是什么?卷?设备方向(portrait-landscape)?

让我们看一下每个计算,看看它真正计算的是什么。

对于azimuth官方计算是

azimuth = (float)Math.atan2(R[1], R[4]); 

让我们回到Y在Device坐标中的写法。如前所述,Y 的坐标在 Device 基础中是 (0, 1, 0) 而 Y 在 World 基础中的坐标是 ( a[1], a[4], a[7])

Y投影到World XY平面的World basis坐标为(a[1], a[4])。这个投影向量和N向量的夹角是(画出投影向量的图,不懂就N

Math.atan2(R[1], R[4]);

这正是azimuth的计算结果。

因此getOrientation计算设备y-axis投射到世界XY-plane与世界北轴之间的角度。因此,如果设备垂直放置,则此计算没有意义,因为投影的 y-coordinate 始终为 0。从几何上讲,如果您指向天空,则计算正北方向没有意义。此外,在这种情况下,围绕 -Z 轴旋转意味着远离纵向旋转,因此文档不正确。

对于Pitch官方计算是

pitch = (float)Math.asin(-R[7]);

Y向量投影到世界N-SKY平面的世界基础坐标是(a[4],a[7])和此投影向量与 Y 向量到世界 E-N 平面的投影之间的角度与 Z[ 的投影之间的角度相同=385=]矢量进入Y-Z平面,重力矢量为(自己再画图)

Math.asin(-R[7]) 

因此,间距是设备 y-axis 到世界 N-SKY 平面的投影与 Y 到世界的投影之间的角度 E-N飞机

同样,Roll 是设备 x-axis 到世界 X-SKY 平面的投影与 X 的投影之间的角度进入世界E-N位面。

例如我之前指出的。

  • 在人像模式下为平板电脑编写指南针。为什么调用 getRotation 后需要加 90 度才能得到正确答案。

在这种情况下,我们想要计算 x-axis 到世界 XY 平面的投影之间的角度,因此需要加上 90 度,因为 x 轴和 y 轴之间的角度始终为 90 度.人们可以改为获取 X 向量的世界基础坐标,即 (a[0], a[3], a[6]),然后投影到世界 XY 平面中以获得 (a[0] , a[3]) 并进行计算 Math.atan2(R[0], R[3]).

  • 为什么要先调用RemapCoordinateSystem再调用getOrientation获取后置摄像头的方向(设备的负方向z-axis)。

在这种情况下,您要计算 z 轴的方向,而不是 getOrientation 计算的 y 轴。因此,您必须调用 remapCoordinateSystem 将 Z 轴映射到 Y 轴。在几何上,您所做的是将 z-axis 旋转到 y-axix,现在 y-axis 变为 -z-axis。因此,您需要取反旋转矩阵中的 y-column 以返回原始定义坐标系。您可以只获取 -Z 的世界基础坐标,即 (-a[2], -a[5], -a[8]) 然后投影到世界 XY 平面得到 (-a[2], -a[5]) 并进行计算。

注意: 我第一次使用传感器的 classes 是在 2011 年,在知道 getRotationMatrixgetOrientation 做,我创建了一个库来仅使用旋转矩阵进行计算。我不认为 remapCoordinateSystemgetOrientation 是必需的,但直到我写下这个答案我才意识到它们不好用,因为它会给出正确的 azymuth 但随后给出错误pitch 的值。 remapCoordinateSystem的作用是让getOrientation中的azymuth按照你要计算的方向进行正确的计算。但是 pitch 值不正确,您必须为其添加 -90 度以获得后置摄像头方向情况下的正确值以及其他情况下的任何值。

如果我要写这个 class 我会创建

float getPitch(float[] rotMatrix)   
float getRoll(float[] rotMatrix) 
float getOrientation(float[] rotMatrix, int axis_you_want_to_calculate) 

那么就不会出现混淆,所有的答案都是正确的。

  • 为什么在调用RemapCoordinateSystem然后getOrientation正确获取设备z-axis的方向后,当设备是平的时结果就不再有意义了。

您正在尝试计算正北的 -z 轴方向,现在它指向天空,因此计算不再有意义。

  • 如何独立于设备位置(即平坦或不平坦)计算围绕任何设备轴的旋转。

明天我会post回答在Get Euler Yaw angle of an Android Device处绕z轴旋转的情况。

注意:如果您只是想编写一个指南针应用程序并且设备始终是平的,那么 TYPE_ORIENTATION 仍然不错。

注2:本人画画很烂,无法post任何图片说明。我什至无法使用 Gimp 绘制直线(我按照说明按住 shift 键,但我的线条参差不齐)。如果有人擅长绘画并愿意提供帮助,请发表评论,我会说明我想到的图片插图以及将它们放在哪里。