旋转矩阵转欧拉角, 及在3D人脸中的使用

长虹剑

已于 2024-01-22 19:52:55 修改

阅读量2k

点赞数 2

分类专栏：视觉图形学文章标签：旋转矩阵欧拉角

于 2019-08-18 00:24:21 首次发布

本文链接：https://blog.csdn.net/hongmaodaxia/article/details/99699496

版权

视觉图形学专栏收录该内容

13 篇文章 4 订阅

订阅专栏

被欧拉角折磨了几年了… 这次还是记一下吧

文章目录

旋转的基础知识总结
旋转矩阵转欧拉角
3D 人脸最好使用的欧拉角为 ...
通过相机位置算出旋转矩阵的验证
sympy 帮助计算旋转矩阵
多年之后理解，内外旋转【最应该先看这个】
坐标系变换，对应旋转与平移形变
- 理论：
bvh 理解旋转
深入了解旋转和内外旋转
UE FRotator 大坑

旋转的基础知识总结

一般表示旋转有以下几种方式

旋转矩阵 R(似乎对应唯一的旋转, 可以看这里的说明). 唯一很重要, 一般表示旋转都可以转到R, 然后变到其他(本博客主要说明变到欧拉角ZXY)
欧拉角 X, Y, Z (pitch, yaw, roll). 但是它们的组合可以达到12种,因为还可以使用类似 X1 Z2 X3 这种(就是某个方向用两次, 因为每次都是相对旋转) 具体可看 wiki. 因此使用欧拉角还必须明确使用的角度和使用的顺序.
轴角. 不过多介绍属于李代数so3. 这个也是经常用到的, 需要借用罗德里格斯转到 R
四元数. 计算非常方便，坐标系变换通过直接改变虚部数值实现
相机位置, head 朝向(一般假设[0, 1, 0]). (有6D rotation 的感觉)这个是有特殊应用场景的, 比如opengl 中的 lookat 函数. 简单介绍下如何算出Rc2w (相机到世界坐标系的旋转), x=R[:, 2]可以通过相机位置和坐标原点的向量获得, 然后y=R[:, 1] 在z和 head 所在平面上, 可利用这个信息算出 x = R[:, 0], 然后再算出y. 代码如下

# 用 torch 实现的. 在ogl 中 z = cam_at 
def ogl_cam_R(z):
    y_t = z.new(3).zero_()
    #y_t.fill_(0)
    y_t[1]=1
    x = y_t.cross(z)
    y = z.cross(x)
    xyz = [ (e/e.norm()).view(-1, 1) for e in [x, y, z] ]
    R = torch.cat(xyz, 1)
    return R

旋转矩阵转欧拉角

前面说过 R 的形式是由你的欧拉角确定的, 要从 R 恢复欧拉角, 首先要知道你要确定你的欧拉角是怎样的定义.
接下来先做一个非常有必要的规定

旋转矩阵乘在顶点的左边即 $t v = R * v$
给定欧拉角的顺序1->2>3 那么其组成的R 按照 $R=R_3R_2R_1$

(其实我自己之所以感觉绕就是因为各种博客和paper 没有明确给出其内在的顺序, 或者是自己没看到)

可以有不同的定义方式如果你没有搞清楚内在的定义顺序你将很糊涂比如

示例1: http://pages.mtu.edu/~dcclark/dice/Background/eulerangles2.pdf

示例2: https://www1.udel.edu/biology/rosewc/kaap686/notes/matrices_rotations_2.pdf

示例3: https://en.wikipedia.org/wiki/Euler_angles

示例3定义的顺序与我的定义相同
至于怎么从R转欧拉角, 当你定义好顺序之后, 看这些矩阵的表达式你将很容易恢复出那些角度来. 比如示例3中我写为 $R_{zxy}$ , 其要恢复x, y 如下

# R=Y \times X \times Z
def get_yaw_pitch_from_Rzxy(self, R):
    view = R[:, 2]
    yaw = math.atan2(view[0], view[2])
    pitch = math.asin(-view[1])
    # roll = math.atan2(R[1, 0], R[1, 1])
    return yaw, pitch

3D 人脸最好使用的欧拉角为 …

先看完上一个小节, 然后本标题的 … 为 $R_{zxy}$
为什么?

对于人脸一般最主要的就是 yaw, pitch,(yaw更重要) 做3D 和 2D 特征点的同学应该都知道, 因为3D特征点需要根据角度动态选择. 选用 $R_{zxy}$ 可以很快算出这两个角.
核能知识点
最先旋转的角在最后呈现的模型中越看不出来, 越是后面旋转的角在最后的模型中越能看出来.
$R_{zxy}$ 的本质是按照roll -> pitch -> yaw 旋转, 也就是最后的 yaw 可以通过打开模型一眼看出来, 就是基本能知道yaw的正负和大小. 方便调试等.

通过相机位置算出旋转矩阵的验证

这个主要是为了记录一下自己总结的plt的使用, 可以忽略

import matplotlib.pyplot as plt
# 2018-3-4
# o p 可以快捷操作图片
# 终于找到退出程序的方法了
def on_release(event):
    if event.button == 3:  plt.close()
def on_press(event):
    if event.dblclick == True:  exit(0)

def create_fig():
    fig = plt.figure()
    fig.canvas.mpl_connect("button_press_event", on_press)
    fig.canvas.mpl_connect("button_release_event", on_release)
    return fig

# 给定Rw2c

cam_direct = Rw2c[:, 2]
R = ogl_cam_R(cam_direct)
fig = create_fig()
ax = plt.subplot(111, projection='3d')

bg = torch.zeros(1, 3)
# 前面的长一些, 方便能看到
for i in range(3):
    pts = torch.cat((bg, Rw2c[:, i].view(1, -1)*1.1), 0).numpy()
    ax.plot(pts[:, 0], pts[:, 1], pts[:, 2], c='g')

for i in range(3):
    pts = torch.cat((bg, R[:, i].view(1, -1)), 0).numpy()
    ax.plot(pts[:, 0], pts[:, 1], pts[:, 2], c='r')

ax.set_zlabel('Z')  # 坐标轴
ax.set_ylabel('Y')
ax.set_xlabel('X')
plt.show()

# 可以防止阻塞
#plt.draw()
#plt.pause(.001)

sympy 帮助计算旋转矩阵

import sympy

#x, y, z=sympy.Symbol('x y z')
x, y, z=sympy.symbols('x y z')
mx = sympy.Matrix([[1, 0, 0], [0, sympy.cos(x), -sympy.sin(x)], [0, sympy.sin(x), sympy.cos(x)]])
my = sympy.Matrix([[sympy.cos(y), 0, sympy.sin(y)], [0, 1, 0], [-sympy.sin(y), 0, sympy.cos(y)] ] ) 
mz = sympy.Matrix([[sympy.cos(z), -sympy.sin(z), 0], [sympy.sin(z), sympy.cos(z), 0], [0, 0, 1]])

#R=mx*my*mz
R=mz*my*mx

#p( R )
#p( sympy.latex(R) )

my_func=sympy.Function('f') # 可以用这种手段把上面的式子变的更加简单, 比如 sin=s, cos=c

p( my_func(R) )

#p( R.diff(x) )

输出

f(Matrix([
[cos(y)*cos(z), sin(x)*sin(y)*cos(z) - sin(z)*cos(x),  sin(x)*sin(z) + sin(y)*cos(x)*cos(z)],
[sin(z)*cos(y), sin(x)*sin(y)*sin(z) + cos(x)*cos(z), -sin(x)*cos(z) + sin(y)*sin(z)*cos(x)],
[      -sin(y),                        sin(x)*cos(y),                         cos(x)*cos(y)]]))

多年之后理解，内外旋转【最应该先看这个】

Intrinsic and extrinsic rotations 这次主要结合这个

欧拉角分为，内在和外在旋转。 intrinsic rotation
一组内在旋转表示相对于对象空间的旋转，每次旋转后都会发生变化。
一组外在旋转表示相对于固定坐标系（通常是世界坐标系）的旋转
在这里插入图片描述
就记住旋转矩阵按照顺序写的是内，按照顺序依次执行旋转的是外
由此看出，内在旋转满足之前说的，最后的旋转保留了一开始的坐标系值。

scipy 里面的 rotation 是根据大小写区分内外的
在这里插入图片描述

大写就是内，这个和 pytorch3d 里面的一样, 而小写是外，从左到右旋转，其实比较符合直觉。
bvh 是内具体介绍在这里

当角度表示围绕三个不同轴的旋转时，一系列三个基本旋转称为Tail-Bryan 角
在这里插入图片描述
正是因为 bvh 是内，所以转到 smpl 系数的时候需要

eulr_order='XYZ'
rot = sciR.from_euler(eulr_order, angles, degrees=True).as_rotvec()

而在 blender 中, smpl 2 blender

eulr = sciR.from_rotvec( vec ).as_euler('xyz')

旋转要规定好旋转的角度
为什么绕y旋转和绕x,z 旋转不一样，这是因为绕y的时候就顺时针了(要从z绕到x)，所以是负数。

其实规定是很重要的，包括极坐标的使用，我还见过这样定义的极坐标，和以前以及网上认识的都不一样(坐标系，角度的表示，角度的标识位置)
在这里插入图片描述

坐标系变换，对应旋转与平移形变

问题来源：
其实很多地方会遇到
1） unity 是左手坐标系，我们一般是右手，涉及到变换，怎么变到unity 中
2）blender 不会处理自动骨骼坐标系之后内部的坐标系应用。（就是自动了之后保存模型的时候坐标系已经变了）
在这里插入图片描述
这样的坐标系虽然便于制作IK，但是默认的 smpl 坐标系都是标准的 ogl 坐标系 这样会导致同样的系数难以产生符合要求的动画

理论：

结论1：坐标系变换等于顶点变换的逆，当然一般容易知道坐标变换，所以应该是
顶点变换是坐标系变换的逆 $P_2= E^{-1}_{1to2} P_1$ , $P_2 = E_{2to1} P_1$
举个粒子
在这里插入图片描述
对于正交矩阵，逆就是转置，所以 $E_{2to1}$ 很直接，两个P 的值也很容易得到

结论2：对应在 $E_1$ 中的旋转变换 $R_1$ ，对应到 $E_2$ 中为，(注意E是坐标系变换，和点变换时互逆的)
$R_2 P_2 = E_{2to1} R_1 E_{1to2} P_2 ⇒ R_2 = E^TR_1E$
这里写的时候带了P，这个是为了直观明白为什么是上面那种对应关系.
下面是做 smpl 系数的时候用到的

_pose = np.einsum( "njab,jbc->njac", pose, bone_mat )
_pose = np.einsum( "jba,njbc->njac", bone_mat, _pose )
pose = _pose

结论3：对应在 $E_1$ 中的平移变换 $t_1$ ，对应到 $E_2$ 中为，
$t_2 + P_2 = E_{2to1} (t_1 + E_{1to2} P_2) ⇒ t_2 = E^Tt_1$

bvh 理解旋转

bvh quat 中有，这个 to_euler 揭示出bvh 内旋的信息

def to_euler(x, order='zyx'):
    x0, x1, x2, x3 = x[..., 0:1], x[..., 1:2], x[..., 2:3], x[..., 3:4]

    # 这个其实是按照 RzRyRx v 即 x,y,z 依次旋转的结果但是 return 反了，是 z-y-x。
    # 然而由于bvh 是intrinsic ，即执行的时候是 v RzRyRx , 其实还是相当于extrinsic 的 RzRyRx v 执行 
    if order == 'zyx':
        return np.concatenate([
            np.arctan2(2.0 * (x0 * x3 + x1 * x2), 1.0 - 2.0 * (x2 * x2 + x3 * x3)),
            np.arcsin(np.clip(2.0 * (x0 * x2 - x3 * x1), -1.0, 1.0)),
            np.arctan2(2.0 * (x0 * x1 + x2 * x3), 1.0 - 2.0 * (x1 * x1 + x2 * x2)),
        ], axis=-1)
    elif order == 'xzy':
        return np.concatenate([
            np.arctan2(2.0 * (x1 * x0 - x2 * x3), -x1 * x1 + x2 * x2 - x3 * x3 + x0 * x0),
            np.arctan2(2.0 * (x2 * x0 - x1 * x3), x1 * x1 - x2 * x2 - x3 * x3 + x0 * x0),
            np.arcsin(np.clip(2.0 * (x1 * x2 + x3 * x0), -1.0, 1.0))
        ], axis=-1)
    else:
        raise NotImplementedError('Cannot convert to ordering %s' % order)

注意上面的四元素第一个是 w，但是 sci rotation 的四元数是 (x, y, z, w)。

深入了解旋转和内外旋转

通过欧拉角构造旋转矩阵的时候，要按照 |x y z| 这样看坐标系（无所谓左右手，因为和你的向量有关）：规定从轴指向往下看轴变换，逆时针为正，都是按照行列式 |x y z| 排列。比如x: 看 y->z， y: 看 x->z ， z：看 x->y ，而没有其他区分。
欧拉角：就完全是看轴和正方向，不看其他的
轴角：要区分个左右手看正方向，即右手坐标系逆时针为正，左手是顺时针为正。

内外旋的问题：下面 a,b等价，其实就是把旋转矩阵反着写一下

angles=[12,23,56]
a=sciR.from_euler('XYZ', angles, degrees=True)
b=sciR.from_euler('zyx', angles[::-1], degrees=True)

UE FRotator 大坑

UE 中可以自己改变一下 x,y,z 看一下其造成的变换是啥，具体就是如下，大体上来说，x,y 都是右手逆时针为正，但是到了y 就成了右手顺时针为正。于是造成了疑惑
在这里插入图片描述

首先要知道 FRotator 是按照 yzx 存放，然后 zyx 依次执行（即RxRyRz v），写一个调试代码

FVector Axis; float Angle; FRotator Rotator;
FQuat(FRotator(0, 90, 0)).ToAxisAndAngle(Axis, Angle);
UE_LOG(LogTemp, Warning, TEXT("axis angle z: %.3f %.3f %.3f %.3f"), Angle, Axis[0], Axis[1], Axis[2]);
FQuat(FRotator(90, 0, 0)).ToAxisAndAngle(Axis, Angle);
UE_LOG(LogTemp, Warning, TEXT("axis angle y: %.3f %.3f %.3f %.3f"), Angle, Axis[0], Axis[1], Axis[2]);
FQuat(FRotator(0, 0, 90)).ToAxisAndAngle(Axis, Angle);
UE_LOG(LogTemp, Warning, TEXT("axis angle x: %.3f %.3f %.3f %.3f"), Angle, Axis[0], Axis[1], Axis[2]);

FMatrix RotationMatrix; FString MatrixString;
RotationMatrix = FQuat(FRotator(0, 90, 0)).ToMatrix();
MatrixString = RotationMatrix.ToString();
UE_LOG(LogTemp, Warning, TEXT("Rotation Matrix z:\n%s"), *MatrixString);
RotationMatrix = FQuat(FRotator(90, 0, 0)).ToMatrix();
MatrixString = RotationMatrix.ToString();
UE_LOG(LogTemp, Warning, TEXT("Rotation Matrix y:\n%s"), *MatrixString);
RotationMatrix = FQuat(FRotator(0, 0, 90)).ToMatrix();
MatrixString = RotationMatrix.ToString();
UE_LOG(LogTemp, Warning, TEXT("Rotation Matrix x:\n%s"), *MatrixString);

其结果如下

LogTemp: Warning: axis angle z: 1.571 0.000 -0.000 1.000 【其实可以看做定义，就是定义好 FRotator 怎么转，这个 axis angle 反应出来】
LogTemp: Warning: axis angle y: 1.571 0.000 -1.000 0.000
LogTemp: Warning: axis angle x: 1.571 -1.000 -0.000 0.000
LogTemp: Warning: Rotation Matrix z: 【看成负的，估计是左手的原因】
[-0 1 0 0] 
[-1 -0 0 0] 
[0 -0 1 0] [0 0 0 1] 
LogTemp: Warning: Rotation Matrix y: 【看成负的】
[-0 0 1 0] 
[-0 1 0 0] 
[-1 -0 -0 0] [0 0 0 1] 
LogTemp: Warning: Rotation Matrix x: 【看成负的，注意还有个 z->y 与排列是反的】
[1 0 0 0] 
[0 -0 -1 0] 
[-0 1 -0 0] [0 0 0 1]

其实轴角和旋转矩阵都可以看出来，内部其实是顺时针为正，都是统一的，因此就是 FRotator 造成的推导不一致，属于一个定义问题。