三角函数

计算正余弦函数 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|` 足够小时就中止迭代.

结果约为: 63.43494882292201 度

平面几何

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

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

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

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

判断线段是否相交 (解方程法) 直接计算两直线的交点坐标, 把它表示为参数 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] }. 两个包围盒相交, 当且仅当它们在各个维度的投影都相交.

    有向包围盒相交判断 (OBB, Oriented Bounding Box) 设 `A, B` 是平面上任意位置的两个矩形, 规定矩形每一条边的法向量指向矩形内部, 于是直线的正侧代表内部, 负侧代表外部. 两矩形 `A`, `B` 不相交, 当且仅当下面两条之一成立:
  1. 存在 `A` 的一边 `l_A`, 使得 `B` 的所有顶点位于 `l_A` 的负侧;
  2. 存在 `B` 的一边 `l_B`, 使得 `A` 的所有顶点位于 `l_B` 的负侧.
    3. 若上面两条都不成立则说明 `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)
}

凸包 (convex hull) 是指包围所有给定点的最小凸多边形. 在 Andrew 算法中, 每个点入栈两次, 出栈不超过两次, 复杂度为 `O(n)`. 加上排序, 总复杂度 `O(n log(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. 点-长方体 将长方体的棱无限延长, 空间被分为 27 个区域, 其中 8 个区域靠近顶点, 12 个区域靠近棱, 6 个区域靠近面, 还有 1 个区域在内部. 为判断 p 点落在哪个区域中, 取相互正交的 3 条棱, 用点到线段的算法分别判断 p 与三条棱的位置关系即可.
    1. 如果 p 靠近顶点, 问题化为求两点的距离;
    2. 如果 p 靠近棱, 问题化为求点到直线的距离;
    3. 如果 p 靠近面, 问题化为求点到平面的距离;
    4. 如果 p 在长方体内部, 需要计算点到每个面的距离, 再取最小.
  6. 点-三角形 首先将点 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), 如四叉树, 八叉树等. 我们将待检测对象按空间位置区分管理, 每次只需检测同组对象的碰撞即可.