Luna Flow
简体中文

core 教程

本教程借助 Luna-Flow/quaternion，带你从空项目出发，完成向量旋转、姿态的复合与插值，以及与欧拉角之间的转换。最后介绍基于 luna-generic trait 的泛型代码以及容易犯的错误。这里尽量少讲数学；推导见 core 设计。

快速入门

把模块加入项目：

moon add Luna-Flow/quaternion@0.2.0

在 moon.pkg 中导入本包。示例还用到 luna-generic（用于泛型代码）与 moonbitlang/core/math（用于 π\pi）：

import {
  "Luna-Flow/quaternion",
  "Luna-Flow/luna-generic" @lf_alg,
  "moonbitlang/core/math",
}

最小的实用程序：把 x 轴绕 z 轴旋转四分之一圈：

fn main {
  let quarter_turn = @quaternion.from_axis_angle((0.0, 0.0, 1.0), @math.PI / 2.0)
  println(quarter_turn)
  let (x, y, z) = quarter_turn.rotate((1.0, 0.0, 0.0))
  println("(\{x}, \{y}, \{z})")
}

它输出：

0.7071067811865476 + 0i + 0j + 0.7071067811865475k
(2.220446049250313e-16, 1, 0)

x 轴落在 y 轴上，第一个坐标有约 2×10−162 \times 10^{-16} 的舍入误差。

日常任务

构造四元数并相乘

Quaternion::from_vec((w, x, y, z)) 构造 w+xi+yj+zkw + xi + yj + zk。常用运算符都可使用，* 是 Hamilton 积。基元的乘法满足 ij=kij = k 而 ji=−kji = -k，因此因子顺序很重要：

test "units" {
  let i = @quaternion.Quaternion::from_vec((0, 1, 0, 0))
  let j = @quaternion.Quaternion::from_vec((0, 0, 1, 0))
  let k = @quaternion.Quaternion::from_vec((0, 0, 0, 1))
  inspect(i * j, content="0 + 0i + 0j + 1k")
  inspect(j * i, content="0 + 0i + 0j + -1k")
  inspect(i * i, content="-1 + 0i + 0j + 0k")
  inspect(i * j * k, content="-1 + 0i + 0j + 0k")
}

整数分量是精确的，这使 Quaternion[Int] 很适合验证恒等式。旋转请使用 Double 分量。

旋转向量并复合旋转

from_axis_angle(axis, angle) 给出旋转的单位四元数，rotate 把它作用于三维向量。若要先施加 first 再施加 second，就像矩阵作用于列向量那样，将它们相乘为 second * first：

test "compose rotations" {
  let first = @quaternion.from_axis_angle((0.0, 0.0, 1.0), @math.PI / 2.0) // z, 90 degrees
  let second = @quaternion.from_axis_angle((1.0, 0.0, 0.0), @math.PI / 2.0) // x, 90 degrees
  let both = second * first
  let (x, y, z) = both.rotate((1.0, 0.0, 0.0))
  // x -> y under `first`, then y -> z under `second`
  assert_true(x.abs() < 1.0e-15 && y.abs() < 1.0e-15 && (z - 1.0).abs() < 1.0e-15)
  let (x2, y2, z2) = second.rotate(first.rotate((1.0, 0.0, 0.0)))
  assert_true((x - x2).abs() < 1.0e-15 && (y - y2).abs() < 1.0e-15 && (z - z2).abs() < 1.0e-15)
}

rotate 假定输入为单位四元数。from_axis_angle 构造的四元数是单位的；自行拼装的四元数请先调用 normalize()。

撤销旋转，以及求两个姿态之间的旋转

单位四元数的逆就是它的共轭。要求把姿态 from 转到姿态 to 的旋转 delta（即 delta * from == to），在右侧做除法：to / from 即 to⋅from−1to \cdot from^{-1}。left_div 求解另一个方程 from * x == to，得到在 from 坐标系中表示的旋转：

test "relative rotation" {
  let from = @quaternion.from_axis_angle((0.0, 1.0, 0.0), 0.3)
  let to = @quaternion.from_axis_angle((1.0, 1.0, 0.0), 1.2)
  let delta = to / from
  assert_true((delta * from - to).square_len() < 1.0e-30)
  let in_from_frame = to.left_div(from)
  assert_true((from * in_from_frame - to).square_len() < 1.0e-30)
  // a unit quaternion's inverse is its conjugate
  assert_true((from.inv() - from.conjugate()).square_len() < 1.0e-30)
}

比较两个旋转

q 与 -q 表示同一个旋转，因此 == 不适合用来判断“姿态相同”。对单位四元数，∣q⋅r∣=cos⁡(θ/2)|q \cdot r| = \cos(\theta/2)，其中 θ\theta 为两个姿态之间的夹角：

///|
fn rotation_angle_between(
  q : @quaternion.Quaternion[Double],
  r : @quaternion.Quaternion[Double],
) -> Double {
  let c = q.normalize().dot(r.normalize()).abs()
  2.0 * @math.acos(if c > 1.0 { 1.0 } else { c })
}

test "compare rotations" {
  let q = @quaternion.from_axis_angle((0.0, 0.0, 1.0), 0.5)
  assert_true(q != -q)
  inspect(rotation_angle_between(q, -q), content="0")
  let r = @quaternion.from_axis_angle((0.0, 0.0, 1.0), 0.75)
  assert_true((rotation_angle_between(q, r) - 0.25).abs() < 1.0e-12)
}

在姿态之间插值

slerp(q1, q2, t) 以恒定角速度沿最短的大圆弧运动：

test "slerp steps" {
  let start = @quaternion.id()
  let end = @quaternion.from_axis_angle((0.0, 0.0, 1.0), @math.PI / 2.0)
  for step in 0..<=4 {
    let t = step.to_double() / 4.0
    let (x, y, _) = @quaternion.slerp(start, end, t).rotate((1.0, 0.0, 0.0))
    let degrees = @math.atan2(y, x) * 180.0 / @math.PI
    assert_true((degrees - 22.5 * step.to_double()).abs() < 1.0e-9)
  }
}

x 轴每步转过 22.5∘22.5^\circ，正如恒定角速度所保证的那样。

与欧拉角相互转换

from_euler(roll, pitch, yaw) 默认使用顺序 "XYZ"，即 qX(roll) qY(pitch) qZ(yaw)q_X(\text{roll})\,q_Y(\text{pitch})\,q_Z(\text{yaw})：先绕 X、再绕新的 Y、最后绕新的 Z 的内旋。对应的逆为 to_euler(order="XYZ", external=false)。to_euler 的默认值是另一种约定——外旋 Z-Y-X，它以相反顺序返回同样的三个角：

test "euler round trip" {
  let q = @quaternion.from_euler(0.1, 0.2, 0.3)
  let (roll, pitch, yaw) = q.to_euler(order="XYZ", external=false)
  assert_true((roll - 0.1).abs() < 1.0e-12)
  assert_true((pitch - 0.2).abs() < 1.0e-12)
  assert_true((yaw - 0.3).abs() < 1.0e-12)
  let (about_z, _, about_x) = q.to_euler() // extrinsic ZYX
  assert_true((about_z - 0.3).abs() < 1.0e-12 && (about_x - 0.1).abs() < 1.0e-12)
}

to_euler 返回的第 k 个角总是属于 order 的第 k 个轴。使用其他顺序前请先阅读 API 说明：部分组合存在已知问题。

读取分量

Quaternion 是抽象类型，没有字段访问器。与基四元数做点积可以读出一个分量，Debug 则显示全部四个分量：

///|
fn components(q : @quaternion.Quaternion[Double]) -> (Double, Double, Double, Double) {
  let e = (w, x, y, z) => q.dot(@quaternion.Quaternion::from_vec((w, x, y, z)))
  (e(1.0, 0.0, 0.0, 0.0), e(0.0, 1.0, 0.0, 0.0), e(0.0, 0.0, 1.0, 0.0), e(0.0, 0.0, 0.0, 1.0))
}

test "components" {
  let q = @quaternion.Quaternion::from_vec((0.5, -1.5, 2.0, 4.0))
  debug_inspect(components(q), content="(0.5, -1.5, 2, 4)")
  debug_inspect(q, content="{ r: 0.5, vec: (-1.5, 2, 4) }")
}

对有限的分量，与基四元数的点积会精确返回该分量。

进阶

基于 luna-generic trait 的泛型代码

只要 T 实现了 luna-generic 的 Ring，Quaternion[T] 也实现它，因此环上的泛型算法接受四元数。下面的函数对任意 Ring 用反复平方计算 xnx^n：

///|
fn[R : @lf_alg.Ring] power(x : R, n : Int) -> R {
  if n == 0 {
    @lf_alg.One::one()
  } else {
    let half = power(x, n / 2)
    if n % 2 == 0 { half * half } else { x * half * half }
  }
}

test "generic power" {
  inspect(power(3, 4), content="81")
  let q = @quaternion.Quaternion::from_vec((1, 1, 1, 1))
  inspect(power(q, 3), content="-8 + 0i + 0j + 0k")
  assert_eq(power(q, 3), q.pow_by_int(3))
}

四元数不是 Field，因此要求 Field（并可能假定 ab=baab = ba）的算法不接受它们。这是有意为之；原因见设计页面。

Int 上的精确恒等式

精确的分量使范数恒等式 ∣pq∣2=∣p∣2∣q∣2|pq|^2 = |p|^2|q|^2 可以在没有舍入的情况下验证。在 Int 上它就是欧拉四平方恒等式：两个四平方和之积仍是四平方和，而 Hamilton 积告诉你是哪一个：

test "four squares" {
  let p = @quaternion.Quaternion::from_vec((1, 2, 3, 4)) // 1+4+9+16 = 30
  let q = @quaternion.Quaternion::from_vec((2, 0, 1, 5)) // 4+0+1+25 = 30
  let pq = p * q
  inspect(pq, content="-21 + 15i + -3j + 15k")
  assert_eq(pq.square_len(), p.square_len() * q.square_len()) // 900
}

需要平方根或三角函数的函数（magnitude、normalize、inv 中的除法、slerp、欧拉角转换）在 Int 上会截断；请把 Int 留给环运算。

长的旋转链

Double 中的每次 Hamilton 积都有舍入，因此许多单位四元数之积的范数会偏离 1，nn 次乘积后大约偏离 n⋅10−16n \cdot 10^{-16}。此时 rotate 会轻微缩放向量。请不时重新归一化：

test "renormalize" {
  let step = @quaternion.from_axis_angle((1.0, 2.0, 3.0), 0.001)
  let mut q = @quaternion.id()
  for _ in 0..<100000 {
    q = step * q
  }
  let drift = (q.square_len() - 1.0).abs()
  assert_true(drift < 1.0e-10)
  q = q.normalize()
  assert_true((q.square_len() - 1.0).abs() < 1.0e-15)
}

性能

所有运算都作用于四个标量并在常数时间内完成，只有 pow_by_int 需要 O(log⁡n)O(\log n) 次乘法。q.rotate(v) 的代价是两次叉积，比用完整的 Hamilton 积计算 q * v * q.inv() 更便宜。每个函数都返回新值，不修改任何东西。

常见陷阱

  • / 在 0.2.0 中改变了含义。 q / r 现在是 q r−1q\,r^{-1}（右除）。为 0.1.x 编写、依赖 / 计算 r−1qr^{-1} q 的代码必须改用 q.left_div(r)。
  • 因子顺序。 a * b 先按 b 旋转。交换因子会得到不同的旋转，除非两者绕同一轴转动。
  • rotate 中的非单位四元数。 rotate 不做归一化。对非单位的 qq，其公式并不是 qvq−1q v q^{-1}，结果通常连长度都与输入不同。请先归一化。
  • 对旋转使用 ==。 q 与 -q 是同一个旋转，但值不相等。浮点结果本来也很少完全相等；请按容差比较。
  • 除以零。 inv、/ 与 left_div 都会除以 ∣r∣2|r|^2。零四元数在 Double 上得到 NaN 分量，在 Int 上触发运行时陷阱。normalize 与 pow_by_T 会原样返回零四元数。
  • Int 分量与超越函数。 magnitude、normalize、inv 等在 Int 上会截断；除单位元外，任何四元数的 inv 都为零。
  • 欧拉角约定。 from_euler 默认为 "XYZ"，to_euler 默认为外旋 "ZYX"。请显式写出两者。避免使用 "ZYX" 或 "YZX" 的 from_euler、外旋 "XZY" 或内旋 "YZX" 的 to_euler，以及接近万向节锁（∣pitch∣≈90∘|\text{pitch}| \approx 90^\circ）的结果，这些都有 API 中列出的已知问题。接近万向节锁时，to_euler 还会向标准输出打印一行警告。
  • 分数次幂。 对标量部为负的四元数，pow_by_T 的结果是错误的；当它表示旋转时，请先对其取负。
  • 自定义分量类型。 DoubleConvert 不允许在本包之外实现，因此只有 Int 与 Double 能用于需要它的函数。

后续阅读

  • core API 列出了每个函数的精确签名与边界情况。
  • core 设计 推导了 Hamilton 积、旋转公式、slerp 与欧拉角提取。
  • luna-generic 定义了这里用到的 Ring、Inverse 与 Conjugate trait。