MATLAB,Mathematica研究三维空间中的旋转

放弃更新，都看这篇文档吧 https://krasjet.github.io/quaternion/quaternion.pdf 👍
https://www.zhihu.com/question/23005815/answer/483589712 quaternion.pdf文档作者本人krasjet的知乎回答
放弃更新，可以去看这篇文章 https://edu.uwa4d.com/lesson-detail/188/1051/0
你看这个repo也行 https://github.com/Krasjet/quaternion

三维空间旋转还是挺实用的，可能在游戏制作，3D模型，航海中用到

3维空间中的【绕过原点的轴旋转】的自由度是3，这是因为确定一个转轴的取向需要2个自由度(\(\cos^2\alpha+\cos^2\beta+\cos^2\gamma=1\))，确定绕这个轴转动了多少角度需要1个自由度。我见过不同的描述方法，但是他们都可以用\(3\times3\)的旋转矩阵\(R\)描述。\(\vec{x}^{\prime}=R\vec{x}\)

一般来说你见到的R，行列式都是1。如果是-1的话说明还来了次镜像

RollPitchYawMatrix
EulerMatrix
四元数形式
单次旋转
它们间的变换关系
MATLAB编程练习(✅验证了没问题)
把上面的代码用Mathematica重写一遍(✅验证了没问题)
参考

RollPitchYawMatrix

在全局框架中的轴向旋转

推导矩阵\(R\)可以用mathematica来辅助进行

RollPitchYawMatrix[{\[Alpha], \[Beta], \[Gamma]}, {1, 2, 
    3}] /. {Cos[x_] :> Subscript[c, x], 
   Sin[x_] :> Subscript[s, x]} // MatrixForm

EulerMatrix

在内禀框架(理解为框架长在旋转物体上)的轴向旋转

推导矩阵\(R\)可以用mathematica来辅助进行

(m = EulerMatrix[{\[Alpha], \[Beta], \[Gamma]}, {1, 2, 3}]) /. {Cos[
     x_] :> Subscript[c, x], Sin[x_] :> Subscript[s, x]} // MatrixForm

四元数形式

拿{a,b,c,d}来描述。把旋转对应到范数为1的四元数\(a+b\mathbf{i}+c\mathbf{j}+d\mathbf{k}\)。

满足\(a^{2}+b^{2}+c^{2}+d^{2}=1\),所以实际是3个自由度

\[\vec{x}^{\prime}=\left(\begin{array}{ccc} a^{2}+b^{2}-c^{2}-d^{2} & 2(b c-a d) & 2(b d+a c) \\ 2(b c+a d) & a^{2}+c^{2}-b^{2}-d^{2} & 2(c d-a b) \\ 2(b d-a c) & 2(c d+a b) & a^{2}+d^{2}-b^{2}-c^{2} \end{array}\right) \vec{x} \]

上式写成向量形式就是

\[\vec{x}^{\prime}=\vec{x}+2 a(\vec{\omega} \times \vec{x})+2(\vec{\omega} \times(\vec{\omega} \times \vec{x})) \]

拿四元数来写就是这个公式：

\[\mathbf{p}^{\prime}=\mathbf{q} \mathbf{p q}^{-1} \]

其中，
\( \mathbf{p}=\left(p_{x}, p_{y}, p_{z}\right)=p_{x} \mathbf{i}+p_{y} \mathbf{j}+p_{z} \mathbf{k} \)

\(\mathbf{q}=e^{\frac{\theta}{2}\left(u_{x} \mathbf{i} + u_{y} \mathbf{j}+u_{z} \mathbf{k}\right)}=\cos \frac{\theta}{2}+\left(u_{x} \mathbf{i}+u_{y} \mathbf{j}+u_{z} \mathbf{k}\right) \sin \frac{\theta}{2}\)

\(\mathbf{p}^{\prime}=\left(p_{x}^{\prime}, p_{y}^{\prime}, p_{z}^{\prime}\right)\)

单次旋转

给出空间中一个点，将该点绕单位向量(k1 k2 k3)旋转一定角度theta，从而得到新点

变换矩阵\(\boldsymbol{R}=\boldsymbol{I}+\sin (\theta) \boldsymbol{K}+(1-\cos (\theta)) \boldsymbol{K}^{2}\)

其中\(\boldsymbol{K}=\left[\begin{array}{ccc} 0 & -k_{3} & k_{2} \\ k_{3} & 0 & -k_{1} \\ -k_{2} & k_{1} & 0 \end{array}\right]\)

它们间的变换关系

在内禀框架(理解为框架长在旋转物体上)的轴向旋转方案：三个转动角度
在全局框架中的轴向旋转的方案：三个转动角度
四元数的方案：存4个实数
单次旋转的方案：转轴方向单位向量(k1 k2 k3)，转动角度theta
旋转矩阵的方案：存9个实数，存在冗余

一般你遇到的都是这5种之间转换来转换去

略（会出现死锁吗？）
也可以看这个 https://zhuanlan.zhihu.com/p/45404840

MATLAB编程练习(✅验证了没问题)

计算点\(p\)绕直线\(\frac{x-x_0}{k_x}=\frac{y-y_0}{k_y}=\frac{z-z_0}{k_z}\) （\(\color{red}{\vec{k}\text{是法向单位矢量}}\)）按右手螺旋规则旋转\(\theta\)角得到的新点

main.m

clc;
clear all;
close all;


p=[2,0,0]; %旋转前的点
r0=[1,0,0]; %直线上的一点r0
x0=r0(1);
y0=r0(2);
z0=r0(3);
param=[0,1,0,pi/3];%分别是kx,ky,kz,theta
p_relative=p-r0;  %旋转前的点相对r0的矢量
p_rotated=Rotation(p_relative,param)+r0 %p_rotated就是结果

Rotation.m

function p1=Rotation(p,param)
[row,col]=size(param);
if row==1&col==4 
    ;
else
    warndlg('请检查Rotation函数的第二个参数param的size是否是1*4');
end
if norm(param(1,end-1))-1.0<1e-8
    ;
else
    warndlg('请检查Rotation函数的第二个参数param的kx,ky,kz构成的向量是否模长1')
end   
    
[row,col]=size(p);
if row==1&col==3
    ;
else 
    warndlg('请检查Rotation函数的第一个参数p的size是否是1*3');
end

kx=param(1);
ky=param(2);
kz=param(3);
theta=param(4);

K=[0 -kz ky;
   kz 0 -kx;
   -ky kx 0];
R = eye(3) + sin(theta)*K + (1- cos(theta))*K*K;    %罗德里格斯公式
p1= R*p';                           %利用罗德里格斯公式计算的旋转后的点
p1=p1';
end

把上面的代码用Mathematica重写一遍(✅验证了没问题)

ClearAll["Global`*"];
On[Assert];
(*Define the point to be rotated,the line point,and rotation parameters*)
p = {2, 0, 0};(*Point before rotation*)
r0 = {1, 0, 0}(*Point on the line*)
param = {0, 1, 0, Pi/3};(*Rotation parameters:{kx,ky,kz,theta}*)

(*Define the Rotation function using Rodrigues' rotation formula*)
Rotation[p_, param_] := 
 Module[{kx, ky, kz, theta, K, R, p1}, {kx, ky, kz, theta} = param;
  Assert[kx^2 + ky^2 + kz^2 == 1];
  K = {{0, -kz, ky}, {kz, 0, -kx}, {-ky, kx, 0}};
  R = IdentityMatrix[3] + 
    Sin[theta]*K + (1 - Cos[theta])*MatrixPower[K, 2];
  p1 = R . p;
  p1]

(*Calculate the rotated point*)
pRelative = p - r0;
pRotated = Rotation[pRelative, param] + r0;

(*Display the result*)
pRotated