平面几何

三角函数

计算正余弦函数 sin: 当 `|x|` 充分小时, 近似有 `sin x ~~ x`; 否则使用三倍角公式递推. cos: 当 `|x|` 充分小时, 近似有 `cos x ~~ 1-x^2/2`; 否则使用二倍角公式递推.

Cordic 算法 是一种计算 `arctan x` 的查表法. 事先将 `theta_n = arctan 2^-n`, `n = 0, 1, 2, cdots` 存于表中, 然后逐步旋转逼近. 此方法可以逼近任意一个不超过 `sum_(n ge 0) arctan 2^-n ~~ 99.88^@` 的角度.
从平面向量的旋转公式 `bm v_1 = [cos theta, -sin theta; sin theta, cos theta] bm v` 两边同除以 `cos theta` 得到伪旋转公式: `bm v' = [1, -tan theta; tan theta, 1] bm v`. 取 `bm v_n = (x_n, y_n)`, `tan theta_n = 2^(-n)`. 上式写成 `{ x_(n+1) = x_n - y_n // 2^n; y_(n+1) = y_n + x_n // 2^n :}` 当 `|y_n|` 足够小时就中止迭代.

多边形

判断多边形方向是否为顺时针 (面积法) [来自 github copilot] 用 Green 公式求边界的曲线积分 `int_(del S) -2y dx = 2 iint_S dx dy`, 若结果为负, 说明曲线为顺时针.

判断点是否在多边形内 (射线法, Ray casting) [来自 chatgpt] 从待判断的点向右发射一条射线, 若交点数为奇数, 则点在多边形内. 该算法适用于顺时针和逆时针的多边形, 甚至也适用于自相交的多边形 (重叠部分按异或计算).

判断点是否在多边形内 (回转数算法, Winding number) 把待判断的点与多边形的各顶点相连, 计算相邻两边的有向夹角之和. 若结果为 0, 说明点在多边形外, 否则点在多边形内.

当点落在多边形边界或与某个顶点重合时, 射线法和回转数法均失效, 需作特殊判断.

凸包 (convex hull) 是指包围所有给定点的最小凸多边形. 在 Andrew 算法中, 每个点入栈两次, 出栈不超过两次, 复杂度为 `O(n)`. 加上排序, 总复杂度 `O(n log(n))`.

线段

判断线段是否相交 (解方程法) 直接计算两直线的交点坐标, 把它表示为参数 t 的函数. 两线段相交当且仅当直线相交, 且 0 ≤ t ≤ 1.

判断线段是否相交 (四边形法, 又称跨立实验) 以两条线段为对角线, 张成一个四边形. 则两线段相交当且仅当该四边形是凸四边形. 具体说, 就是 p3, p4 位于直线 p1-p2 的两侧, 且 p1, p2 位于直线 p3-p4 的两侧.

1. 相比于解方程法, 四边形法的优势是避免了除法运算: 众所周知浮点数的除法运算容易引起数值误差.
2. 以上两种算法均没有考虑线段重合 (或部分重合) 的情形. 可以这样改进: 将四边形法判断叉乘的条件从 > 0 改为 >= 0, 但这样一来会将两条共线但不重合的线段判断为相交, 因此还要结合包围盒来判断: 若两线段的包围盒不相交, 则线段也不相交.

包围盒

包围盒相交判断 (AABB, Axis Aligned Bounding Box) 这里的包围盒特指各边与坐标轴都平行的矩形: { min: [x, y], max: [x, y] }. 两个包围盒相交, 当且仅当它们在各个维度的投影都相交.

包围盒与线段相交判断 (Slab 方法, 又称 Liang-Barsky 算法) 设包围盒各边与坐标轴平行, 用一对向量 `(b^(-), b^+)` 表示这个包围盒. 线段则用参数表示为 `P(t) = A + D t`, `t in [0, 1]`.
1. 初始令 `t_min = 0`, `t_max = 1`, 表示线段进入和离开包围盒的时刻.
2. 对每个维度 `i = x, y, z`: 若 `D_i = 0`, 检查 `A_i in [b_i^(-), b_i^+]` 是否成立, 若否, 则线段与包围盒不相交. 若 `D_i != 0`, 解方程 `P_i(t_i^(+-)) = b_i^(+-)` 得到 `t_i^(+-) = (b_i^(+-) - A_i) // D_i`,
   `t_1 = min(t_i^+-)`, `quad t_2 = max(t_i^+-)`. 更新 `t_min = max(t_min, t_1)`, `t_max = min(t_max, t_2)`.
3. 若最终 `t_min le t_max`, 则线段与包围盒相交.
4. 把初始值换成 `t_min = -oo`, `t_max = +oo`, 则可以进行包围盒与直线的相交判断.

有向包围盒相交判断 (OBB, Oriented Bounding Box) 设 `A, B` 是平面上任意位置的两个矩形, 规定矩形每一条边的法向量指向矩形内部, 于是直线的正侧代表内部, 负侧代表外部. 两矩形 `A`, `B` 不相交, 当且仅当下面两条之一成立:
1. 存在 `A` 的一边 `l_A`, 使得 `B` 的所有顶点位于 `l_A` 的负侧;
2. 存在 `B` 的一边 `l_B`, 使得 `A` 的所有顶点位于 `l_B` 的负侧.
若上面两条都不成立则说明 `A, B` 相交. 点与直线的位置关系可以通过计算有向距离得到.
该算法可以推广到平面上的两个凸多边形, 其理论基础是凸集分离定理: 平面上任意两个不相交紧凸集, 存在一条直线将它们分离开. 可以证明, 对于两个不相交凸多边形, 必存在其中一个多边形的一条边将它们分离开. 如果 1, 2 两条均不成立, 说明两个多边形的边均不是 `A, B` 的分割线, 因此 `A, B` 必相交.

生成有向包围盒 (主成分分析) 给定空间中的 `n` 个点, 想要生成它们的有向包围盒, 首先计算这些点的协方差矩阵 `A = ["cov"(x, x), "cov"(x, y), "cov"(x, z); "cov"(y, x), "cov"(y, y), "cov"(y, z); "cov"(z, x), "cov"(z, y), "cov"(z, z)]`, 其中 `"cov"(x, y) = E[(x - E(x))(y - E(y))]`. `A` 是一个半正定的实对称矩阵, 可用 Jacobi 迭代法求出 `A` 的特征向量组成的 3 阶方阵 `T`. `T` 是正交矩阵, 且 `det(T) = 1`, 因此它是一个旋转变换. 现在只需将这些点按 `T^-1` 变换到主方向上, 生成轴向包围盒 AABB, 再将 AABB 加以变换 `T`, 就得到有向包围盒 OBB.

OBB.setfromPoints = (points) => {
  const A = cov(points)
  const T = eig(A).eigvecs
  points = points.map(p => p.clone().applyMatrix3(T))
  this.aabb.setFromPoints(points)
  this.transform.setFromMatrix3(T)
}

包围球 (bounding sphere) [来自 WhereIsHeroFrom, github] 给定平面上的 n 个点, 求它们的最小覆盖圆盘. 它常常是其中 3 个点的外接圆, 偶尔也会以其中 2 点连线为直径. 当 n ≥ 2 时, 至少有 2 点在圆上, 设为 pi 和 pj, 以它们为直径画圆: 如果其它点都在圆内, 我们已经找到最小覆盖圆. 否则枚举第三点 pk, 作 pi, pj, pk 的外接圆, 选取其中最小者即是. 在随机打乱的情况下, 新加入一点落在圆外的概率仅为 3/n, 因此算法的均摊复杂度为 O(n).

矩阵变换

缩放平移边界 设画框是边长为 1 的正方形, 有一幅正方形的图画初始大小位置与画框重合. 现在对图画 (相对于几何中心) 缩放 `k` 倍 (`k gt 1`), 再平移 `p = (p_x, p_y)`. 若要求画框被图画填满 (不能露出空白), 则平移向量 `p` 的范围是多少?

缩放后的图画边长与画框边长之差是 `k - 1`, 即图画的上下左右各有 `(k-1)//2` 的空间容许移动. 注意, 如果先平移再缩放, 则 `p` 的范围要除以 `k`.

缩放中心变换 矩阵 `[s(k), p; 0, 1]` 表示先以原点为中心缩放 `k` 倍, 然后平移 `p = (p_x, p_y)^(sf T)` 的一个线性变换. 其中 `s(k) = [k; , k]`. 简记 `(s, p) = [s, p; 0, 1]`, 用矩阵运算可以得到 `(s_2, p_2) (s_1, p_1)` `= (s_2 s_1, s_2 p_1 + p_2)`. 特别地 `(1, p) (s, 0) = (s, p)`,
`(s, 0) (1, p) = (s, s p)`. 现在先做变换 `(s_1, p_1)`, 再以 `p_2` 为中心缩放 `k_2` 倍. 这第二个步骤相当于先平移 `-p_2`, 再以原点为中心缩放 `k_2` 倍, 然后平移 `p_2`. 最终结果为 `(1, p_2) (s_2, 0) (1, p_2)^-1 (s_1, p_1)` `= (s_2, (1-s_2) p_2) (s_1 p_1)` `= (s_2 s_1, s_2 p_1 + (1 - s_2) p_2)`. 最后这个平移系数可以看作线性插值 `"lerp"(p_1, p_2, 1 - s_2)`.

旋转中心变换 矩阵 `[r(theta), p; 0, 1]` 表示先绕原点旋转角度 `theta`, 然后平移 `p = (p_x, p_y)^(sf T)` 的一个刚体变换, 其中 `r(theta) = [cos theta, -sin theta; sin theta, cos theta]`. 注意, 也可以理解成先从原点做位移 `p`, 再以 `p` 点为中心旋转 `theta` 角度.
下面简记 `(r, p) = [r, p; 0, 1]`, 先后进行两次变换 `(r_1, p_1)`, `(r_2, p_2)`, 用矩阵运算可以得到结果为 `(r_2, p_2) (r_1, p_1) = (r_2 r_1, r_2 p_1 + p_2)`, 特别地 `(1, p) (r, 0) = (r, p)`,
`(r, 0) (1, p) = (r, r p)`. 现在要求第二次变换改为以第一次变换的结果为参考系, 换言之第二次旋转以 `p_1` 为中心. 这时第二次变换 `(r_2, p_2)` 应当调整为相似矩阵 `(r_1, p_1) (r_2, p_2) (r_1, p_1)^-1` `= (r_2, (1-r_2) p_1 + r_1 p_2)`. 验算: `(r_1, p_1) (r_2, p_2) (r_1, p_1)^-1 (r_1, p_1)` `= (r_1, p_1) (r_2, p_2)` `= (r_2 r_1, p_1 + r_1 p_2)`.

立体几何

碰撞检测

简单几何体间的距离 [来自 李雪峰@知乎] 约定: length 是向量长度 dot, cross 分别是点乘与叉乘.
1. 点-点 两点 p1, p2 的距离是 length(p1 - p2). 由于球体就是有厚度的点, 此公式也可以计算两球体的距离.
2. 点-平面 设 n 是平面的单位法向量, q 是平面上一点, 则点 p 到平面的有向距离是 dot(p - q, n), 取绝对值就得到距离. 同上所述, 此公式可以计算球体到平面的距离.
   推而广之, 求凸体到平面的距离, 只需计算每个顶点到平面的有向距离, 如果这些距离有正有负, 说明凸体与平面相交; 否则取其中绝对值最小者, 就等于它到平面的距离. 许多几何体都是凸体: 线段, 长方体等等.
3. 点-直线 设 v 是直线的单位切向量, q 是直线上一点, 则点 p 到直线的有向距离是 cross(p - q, v), 叉乘的结果是一个向量, 取向量的模长就得到距离.
4. 点-线段 点 p 与线段 a-b 的位置关系分成三种:
   1. p 到 a 最近, p-a-b 为钝角, dot(p-a, b-a) < 0;
   2. p 到 b 最近, p-b-a 为钝角, dot(p-b, a-b) < 0;
   3. p 到垂足最近, 计算点到直线距离即可.
   因为胶囊体就是有厚度的线段, 此方法也可以判断点到胶囊体的距离.
5. 点-长方体 对点和长方体施加同一个旋转, 使长方体各边与坐标轴平行. 然后在 `x, y, z` 三个维度上分别求点到线段的距离 `d_x, d_y, d_z`, 例如, 设点的坐标为 `(x, y, z)`, 长方体在 `x` 轴上的最小、最大值分别为 `b_x^(-), b_x^+` 时, `d_x = { x - b_x^+, if x gt b_x^+; b_x^(-) - x, if x lt b_x^(-); 0, otherwise :}` 最终点到长方体的距离为 `sqrt(d_x^2+d_y^2+d_z^2)`.
6. 线段-长方体 假定长方体的边与坐标轴平行, 以一对向量 `(b^(-), b^+)` 表示长方体的范围.
   1. 如果线段长度为 0, 问题化为点到长方体的距离.
   2. 否则, 先执行线段与长方体的相交测试, 若相交则距离为 0.
   3. 将线段参数化: `P(t) = A + (B-A) t = A + D t`, `t in [0, 1]`. 然后在 `x, y, z` 三个维度上分别解方程, 例如 `x` 维度的方程为 `P_x(t_x^(-)) = b_x^(-)`, `quad P_x(t_x^+) = b_x^+`. 收集所有落在 `[0, 1]` 内的 `t` 值并去重, 得到: `t_0 = 0 lt t_1 lt cdots lt t_n = 1`. 这些点将线段分为 `n` 段.
   4. 线段上任一点 `P(t)` 到长方体的距离平方等于 `d^2(t) = d_x^2 + d_y^2 + d_z^2`, 其中 `d_x = { P_x(t) - b_x^+, if P_x(t) gt b_x^+; b_x^(-) - P_x(t), if P_x(t) lt b_x^(-); 0, otherwise :}` `d_x` 虽是分段函数, 但在每个小段 `[t_k, t_(k+1)]` 上总是取固定的分支, 于是函数 `d^2` 分段可微. `"d"/dt d^2(t)` `= "d"/dt sum_i (P_i(t) - b_i^(+-))^2`
      `= 2 sum_i (P_i(t) - b_i^(+-)) P_i'(t)`
      `= 2 sum_i (A_i + D_i t - b_i^(+-)) D_i`. 上式求和只对 `d_i != 0` 的分量进行. 至于 `b_i^(+-)` 的取值, 可以任取区间上一点 (例如中点) 进行判断.
   5. 令 `"d"/dt d^2(t) = 0`, 解得极小值点 `t^ast = (sum_i (b_i^(+-) - A_i) D_i)/(sum_i D_i^2)`. 若 `t^ast in [t_k, t_(k+1)]`, 则它对应着这一小段到长方体的距离. 否则舍弃这个点, 使用区间端点到长方体的距离. 算出每个小段到长方体的距离后, 取最小值即可.
7. 点-三角形 首先将点 p 投影到三角形所在平面, 记它的投影为 q. 设三角形的顶点是 v1, v2, v3, 则 q 可以写为它们的线性组合 q = x1 v1 + x2 v2 + x3 v3, 其中 x1 + x2 + x3 = 1. (x1, x2, x3) 就称为它的重心坐标. 重心坐标的一个重要性质是, 它的 3 个分量之比等于该顶点相对的小三角形的有向面积之比. 下面的公式中 area 代表三角形的有向面积: x1 / area(q, v2, v3) = x2 / area(q, v3, v1) = x3 / area(q, v1, v2) 3 个小三角形有向面积之和等于三角形的面积 A, 所以 A * x1 = area(q, v2, v3) = length(v2 - v3) * d23 / 2, d23 是 q 点到边 v2-v3 的有向距离. 同理可以求出 q 点到其他两边的有向距离. 综上, 点 q 到三角形的距离等于它到 3 个顶点、3 条边的距离的最小值.
   优化: 如果点 q 落在三角形外部, 重心坐标至少有一个分量为负, 说明 q 远离这个分量对应的顶点, 我们可以排除掉这个顶点及其关联的边.

随着场景复杂度的提高, 可以为场景生成 BVH (层次包围体树, Bounding Volume Hierarchy), 如四叉树, 八叉树等. 我们将待检测对象按空间位置区分管理, 每次只需检测同组对象的碰撞即可.