3D数学实例之一：解析法求最大可见平面

前言

3D数学在WebGL中有着非常重要的作用，然而系统性的对数学的学习又显得枯燥和乏味，相信很多同学都在这里被劝退。因此我准备开一个专题，通过工程中遇到的实际数学问题进行解剖，在此过程中培养3D数学思维，学习3D数学。

引出

当我们开发可编辑的3D软件时，通常会使用到坐标轴。一般情况下我们的坐标轴可以使用三维几何画出，如下图。

但有时候，为了更加美观和易用，坐标轴可能需要使用平面（图片）来完成，如下图。在游戏中也常常会有基于平面的交互。

下图为缩放物体的控件，要求 $x 、 y 、 z$ 坐标轴通过平面上的封闭图形画出，且三个轴分别垂直于所选物体包围盒的 $x 、 y 、 z$ 方向，以方便进行不同轴的缩放操作。

three.js采用右手坐标系，对于一个平面，其坐标轴如下图a所示。

图a

首先，我们需要将箭头平面的 $x$ 轴与物体的某一轴对齐，同时将轴移动到物体中心，如果我们要画物体的 $x$ 轴，就将其与 $x$ 轴对齐， $y 、 z$ 两轴亦然。以 $x$ 轴为例：

  mesh.getWorldPosition(_vector3);
  // 将x轴移动到选中物体中心
  axisX.position.copy(_vector3);
  axisX.updateMatrixWorld();
  axis.updateMatrixWorld();
  mesh.updateMatrixWorld();

  // 轴的x方向，转换成世界坐标的方向
  _xDir.set(1, 0, 0).applyMatrix4(axis.matrixWorld);
  _xDir.sub(_vector3);
  _xDir.normalize();

  // 模型的x方向，转换成世界坐标的方向
  _vector2.set(1, 0, 0).applyMatrix4(mesh.matrixWorld);
  _vector2.sub(_vector3).normalize();

  // 轴的x方向到模型的x方向，需要的四元数
  _quaternion.setFromUnitVectors(_xDir, _vector2);

  // 对齐
  axis.quaternion.multiply(_quaternion);

  axis.updateMatrixWorld();

当我们得到了物体 $x$ 轴和箭头平面的 $x$ 轴在世界坐标下的方向后，我们需要将后者与前者对齐，我们可以引入一个四元数，通过起始和目标方向设置，这个四元数即记录了这个变换，然后我们应用这个变换即可，以上完成了坐标对齐的操作。

但是当我们转动视角时，很可能存在下面的问题：某个或多个轴的可见区域变得很小，这是因为我们的箭头不会随着视角进行变化，想象你从纸的横截面方向看它，能看到的只是一根线。如果我们把纸沿着长度方向的轴旋转90°不是能看到了吗？所以我们需要根据视角方向实时调整控件的方向，下面重点讲一下如何去动态的调整这个方向。

问题分析和解决

如何表征可见区域？

我们知道当我们垂直看向平面时看到的区域最大，而平行看过去时区域最小，那么通过什么来衡量这个量呢？稍作分析，1.我们如何表示一个平面？--平面的法线加上平面上的一点即可表示。2.垂直和平行平面有什么区别？--垂直时视线与平面的法线平行，平行时视线与平面的法线垂直。那么我们只需要引入一个变量 $s=∣N⋅V∣s=\mid\mathbf{N}\cdot\mathbf{V}\mid$ ， $N\mathbf{N}$ 为平面的法线， $V\mathbf{V}$ 为视角方向，求这两个方向的点乘即可，我们需要做的就是调整法线 $N\mathbf{N}$ （视角方向为已知量），使得 $s$ 最大。

数学模型

现在我们的问题转化为：当我们将箭头平面与物体指定的轴对齐后，我们要根据视线方向，调整物体沿着自身 $x$ 轴旋转某个角度，使得 $s$ 得到最大值。为何是沿着自身 $x$ 轴旋转？见上图a， $x$ 轴为箭头平面的长度方向，我们最终是让箭头的长度方向对齐物体的某个轴。

上面我们已经得到了对齐操作的四元数和平移变换，为了降维，我们先把平面移动到原点，最后再将其移动到原来的位置，这样在计算过程中就不必考虑平移，且矩阵由4×4转为了3×3。

我们将对齐操作的旋转四元数转为一个矩阵 $M\mathbf{M}$ ，而在此变换之前，我们首先要让平面沿着自身的 $x$ 轴旋转一个角度 $θ\theta$ 。为什么上面的旋转需要在 $M\mathbf{M}$ 变换之前呢？因为 $M\mathbf{M}$ 变换之后平面的位置和旋转变了，其 $x$ 轴也跟着变了，为了方便计算，我们先做 $x$ 轴的旋转变换 $Mrx\mathbf{M_{rx}}$ ，且根据旋转矩阵的性质，可以得到：

$Mrx=[1000cosθ−sinθ0sinθcosθ]\mathbf{M_{rx}}=\left[\begin{array}{c}1 & 0 & 0 \\ 0 & cos\theta & -sin\theta\\ 0 & sin\theta & cos\theta \end{array}\right]$

初始时平面的法线为 $N=[001]\mathbf{N}=\left[\begin{array}{c} 0\\ 0 \\ 1 \end{array}\right]$ ，作用以上两个矩阵后得到的世界坐标系下的法线为:

$N′=MMrxN=[−T12sinθ+T13cosθ−T22sinθ+T23cosθ−T32sinθ+T33cosθ]\mathbf{N^{'}}=\mathbf{M}\mathbf{M_{rx}}\mathbf{N}=\left[\begin{array}{c} -T_{12}sin\theta+T_{13}cos\theta \\ -T_{22}sin\theta+T_{23}cos\theta\\ -T_{32}sin\theta+T_{33}cos\theta \end{array}\right]$

式中 $T_{ij}$ 为 $M\mathbf{M}$ 的 $i$ 行 $j$ 列元素。

因此，先不考虑绝对值时：

$s=N′Vs=\mathbf{N^{'}}\mathbf{V}$

$=(−T12sinθ+T13cosθ)Vx+(−T22sinθ+T23cosθ)Vy+(−T32sinθ+T33cosθ)Vz=(-T_{12}sin\theta+T_{13}cos\theta)V_x+(-T_{22}sin\theta+T_{23}cos\theta)V_y+(-T_{32}sin\theta+T_{33}cos\theta)V_z$

$=(−T12Vx−T22Vy−T32Vz)sinθ+(T13Vx+T23Vy+T33Vz)cosθ=(-T_{12}V_x-T_{22}V_y-T_{32}V_z)sin\theta + (T_{13}V_x+T_{23}V_y+T_{33}V_z)cos\theta$

令：

$a=-T_{12}V_x-T_{22}V_y-T_{32}V_z,b=T_{13}V_x+T_{23}V_y+T_{33}V_z$ ，则：

$s=asinθ+bcosθs=asin\theta+bcos\theta$ ，由高中的数学知识：

$s=a2+b2sin(θ+ϕ)s=\sqrt{a^2+b^2}sin(\theta+\phi)$

$cosϕ=a/a2+b2cos\phi=a/\sqrt{a^2+b^2}$

当我们得到了以上的算式，那么问题就迎刃而解了，即：当 $θ+ϕ=π/2\theta+\phi=\pi/2$ 时， $s$ 得到最大值，由于 $a2+b2>0\sqrt{a^2+b^2}>0$ ，前面的绝对值也不用考虑了。

需要注意的是，我们求得

$θ=π/2−ϕ\theta=\pi/2-\phi$

$ϕ∈[0,2π]\phi\in[0,2\pi]$

而在计算 $ϕ\phi$ 时， $a rccos$ 的值域为 $[0,π][0,\pi]$ ，显然值域无法覆盖所有的情况，我们再看看另一个条件：

$sinϕ=b/a2+b2sin\phi=b/\sqrt{a^2+b^2}$

那么我们就可以结合上述条件得到 $ϕ\phi$ 的值：

ϕ={arcsin(b/a2+b2),sinϕ<0,cosϕ>02π−arccos(a/a2+b2),sinϕ<0,cosϕ<0arccos(a/a2+b2),ohters\phi= \begin{cases} arcsin(b/\sqrt{a^2+b^2}),\quad sin\phi<0,cos\phi>0 \\[2ex] 2\pi-arccos(a/\sqrt{a^2+b^2}),\quad sin\phi<0,cos\phi<0\\[2ex] arccos(a/\sqrt{a^2+b^2}), ohters \end{cases}

在解出 $θ\theta$ 后我们将结果带回原 $Mrx\mathbf{M_{rx}}$ 矩阵，最后即可算出最后的变换。

代码实现

如下为完整代码：

const _vector = new Vector3();
const _vector2 = new Vector3();
const _vector3 = new Vector3();
const _xDir = new Vector3();
const _matrix = new Matrix4();
const _matrix2 = new Matrix4();

const _quaternion = new Quaternion();
const _quaternion2 = new Quaternion();

// 更新X轴
const updateAxisX = () => {
  mesh.getWorldPosition(_vector3);
  axis.position.copy(_vector3);
  axis.quaternion.identity();

  axis.updateMatrixWorld();
  mesh.updateMatrixWorld();

  // 轴的x方向，转换成世界坐标的方向
  _xDir.set(1, 0, 0).applyMatrix4(axis.matrixWorld);
  _xDir.sub(_vector3);
  _xDir.normalize();

  // 模型的x方向，转换成世界坐标的方向
  _vector2.set(1, 0, 0).applyMatrix4(mesh.matrixWorld);
  _vector2.sub(_vector3).normalize();

  // 轴的x方向到模型的x方向，需要的四元数
  _quaternion.setFromUnitVectors(_xDir, _vector2);

  // 对齐
  axis.quaternion.multiply(_quaternion);
  axis.updateMatrixWorld();
  axis.position.set(0, 0, 0);
  _xDir.set(1, 0, 0);

  // // 获取视线方向
  camera.getWorldDirection(_vector);
  _vector.multiplyScalar(-1);
  _vector.normalize();

  _matrix2.copy(axis.matrixWorld);
  const elments = _matrix2.elements;

  const b =
    elments[8] * _vector.x + elments[9] * _vector.y + elments[10] * _vector.z;

  const a =
    - elments[4] * _vector.x - elments[5] * _vector.y - elments[6] * _vector.z;


  if (Math.abs(a) + Math.abs(b) < 1e-8) {
    return;
  }

  const cosPhi = a / Math.sqrt(a * a + b * b);
  const sinPhi = b / Math.sqrt(a * a + b * b);
  let phi = Math.acos(cosPhi);

  // acos(x)值域为0~PI，asin(x)值域为-PI/2~PI/2,两者的值域不重叠
  // 此时x->[-PI/2,0]
  if (sinPhi < 0 && cosPhi > 0) {
    phi = Math.asin(sinPhi);
  } else if (sinPhi < 0 && cosPhi < 0) {
    //此时x->[PI, 3PI/2]
    phi = 2 * Math.PI - phi;
  }

  let theta = Math.PI / 2 - phi;
  _matrix.makeRotationAxis(_xDir, theta);
  _matrix2.multiply(_matrix);
  _quaternion2.setFromRotationMatrix(_matrix);
  axis.quaternion.multiply(_quaternion2);
  // 将轴移动到原来的位置
  axis.position.copy(_vector3);
  axis.updateMatrixWorld();
}