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

Cheanfan
863 阅读5分钟

前言

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

引出

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

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

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

image.png

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

image.png

图a

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

  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();

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

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

image.png

问题分析和解决

如何表征可见区域？

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

数学模型

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

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

我们将对齐操作的旋转四元数转为一个矩阵M\mathbf{M}，而在此变换之前，我们首先要让平面沿着自身的xx轴旋转一个角度θ\theta。为什么上面的旋转需要在M\mathbf{M}变换之前呢？因为M\mathbf{M}变换之后平面的位置和旋转变了，其xx轴也跟着变了，为了方便计算，我们先做xx轴的旋转变换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]

式中TijT_{ij}为M\mathbf{M}的ii行jj列元素。

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

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=−T12Vx−T22Vy−T32Vz,b=T13Vx+T23Vy+T33Vza=-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时，ss得到最大值，由于a2+b2>0\sqrt{a^2+b^2}>0，前面的绝对值也不用考虑了。

需要注意的是，我们求得

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

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

而在计算ϕ\phi时，arccosarccos的值域为[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();
}

总结

以上通过3D数学的知识，从分析问题到数学建模，再到问题求解，是不是很优雅？

avatar
文章
阅读
粉丝