【Frustum Culling】视锥体剪裁数学原理和代码实现
前言
剪裁是渲染中常用的手段，避免将渲染资源浪费在无意义的片段中，在渲染管线的齐次除法，渲染管线就会帮我们做一次剪裁，防止在视锥体外的顶点跑到像素着色器被渲染。
但这终究会进入渲染管线，会进入顶点着色器乃至曲面细分和集合着色器，因此我们最好在进入渲染管线前就先做一次粗糙一些的剪裁，在渲染管线阶段之前的CPU阶段就将物体拒之门外。
最常见的就是给物体一个包围盒，然后我们用视锥体和包围盒进行碰撞检测，假如碰撞失效就不进行渲染。
包围盒
顾名思义，就是包裹当前物体或一组物体的盒子，最常用的就是轴对齐包围盒（AABB），其每个边都与坐标轴平行，没有旋转。
包围盒的表示
包围盒的表达方式不止一种。
其一：记录中点和边长，既：
struct BoundingBox
{
XMFLOAT3 Center; // Center of the box.
XMFLOAT3 Extents; // Distance from the center to each side.
}
其二，记录最大点、最小点：
struct BoundingBox
{
XMFLOAT3 MaxPoint;
XMFLOAT3 MinPoint;
}
两者很容易互相转换：
//1->2
max = center + extents/2;
min = center - extents/2;
//2->1
center = (max+min)/2
extents = max - min
Unity的包围盒类和DirectX的BoundingBox以及我们都是用第一种存储方式。
包围盒的计算
我们可以给物体手动加上包围盒，但这肯定不是最好的选择，毕竟人有失误，做的包围盒肯定有些许误差，所以我们可以让程序自己计算。
下面是Unity的程序实例：
using UnityEngine;
[RequireComponent(typeof(MeshFilter))]
public class ObjData : MonoBehaviour
{
public bool isActive = false;//是否被激活
private bool calc = false;//是否已经计算包围盒
public Bounds bound
{
get
{
if(!calc)//如果没有计算包围盒
{
Mesh mesh = GetComponent<MeshFilter>().sharedMesh;
Vector3 maxP = Vector3.one * float.MinValue;
Vector3 minP = Vector3.one * float.MaxValue;
foreach (Vector3 v in mesh.vertices)//遍历每个顶点
{
Vector4 world = transform.localToWorldMatrix * new Vector4(v.x, v.y, v.z, 1);//将顶点转换到世界空间
//用xyz分量更新最大和最小值
for (int i = 0; i < 3; ++i)
{
maxP[i] = Mathf.Max(maxP[i], world[i]);
minP[i] = Mathf.Min(minP[i], world[i]);
}
}
//构造函数，参数分别是中心center和xyz三边长度
_bound = new Bounds((maxP + minP) / 2, maxP - minP);
calc = true;
}
return _bound;
}
}
public Bounds _bound;
private void OnDrawGizmos()//可以让Unity给我们画出来
{
Gizmos.color = Color.cyan * (isActive ? 1 : 0.5f);//如果激活，则显示天蓝色，否则只有一半亮度
Gizmos.DrawWireCube(bound.center, bound.size + Vector3.one * 0.1f);
}
}
视锥体
顾名思义，就是我们的摄像头，又叫平截锥体、截头椎体、Frustum，我们先把它来表示出来。
我是采用DirectX标准库BoundingFrustum的存储方法：
struct BoundingFrustum
{
XMFLOAT3 Origin; // 位移
XMFLOAT4 Orientation; // 旋转
float RightSlope; // Positive X (X/Z)
float LeftSlope; // Negative X
float TopSlope; // Positive Y (Y/Z)
float BottomSlope; // Negative Y
float Near, Far; // Z of the near plane and far plane.
}
看起来可能有些蒙，视点和旋转和容易理解，我们看下面几个参数。
众所周知，视锥体有6个平面:
近平面很小，这不是椎体
我们想像初始状态，摄像头对着Z轴正半轴（DirectX和Unity的左手坐标系，Opengl是负半轴）。对于远、近平面，其初始状态是垂直于Z轴的，我们只要存储到Z轴的距离即可。而四个平面我们存储的是斜率。如下图：
我们高中做题时，斜率是y除以x（如上图左上角），而三维也差不多，首先看视锥体的右平面，其斜率就是X/Z，上平面斜率时Y/Z，对称平面斜率互为相反数。回过头来说一下视点和旋转是DirectX标准库中的定义，但Unity有Camera对象，我们可以直接通过Camera得到View矩阵，所以我这样组织视锥体对象：
class BoundingFrustum
{
public float RightSlope; // X/Z
public float LeftSlope; //-X/z
public float TopSlope; // Y/Z
public float BottomSlope;//-Y/Z
public float Near, Far;
public Matrix4x4 viewMat;
}
视锥体的构造方法
理想状态下就像上图画的那样，但再加上位移、旋转，思考复杂度瞬间上去了，有没有办法能一直保持为理想状态下呢？当然有！
我们变换空间操作，从模型空间->世界空间->观察者空间->投影空间->齐次除法->标准化设备空间(ndc)，其中观察者空间又称摄像机空间，顾名思义，就是让我们想象中，摄像机是不动的，摄像机的旋转和移动都看做是其他物体的反向移动，因此在摄像机空间下，我们的视锥体就一直是理想初始状态。
我们需要通过什么方法获取初始状态的视锥体？通过上图可知，我们主要需要6个点，分别是：远、近平面与Z轴相交的点、远平面矩形与XOZ平面的两个交点、远平面矩形与YOZ平面的两个交点；通过这六个点，斜率就能算出来。
我们不知道这些点在观察者空间下的坐标，但我们知道NDC空间下的坐标，复习一下，投影变换是为了将视锥体从截头椎体变为一个长方体，做到近大远小的效果，然后再转换到NDC空间来适应设备（不过实际上被分为投影矩阵和齐次除法），NDC空间就是一个[-1,1]x[-1,1]x[0,1]的长方体，如果变换后顶点不在其中，那就说明之前顶点也不在视锥体中，也就表明这个顶点要被剪裁掉，这都是渲染管线帮我们做的。
NDC坐标肯定不会变，那我们只要反向操作，就能得到视锥体的view坐标。
不过这里要提一下，DirectX的NDC空间范围和如上所说，OpenGL的是[-1,1]x[-1,1]x[-1,1]，这也这也就导致变换所需的投影矩阵不同，Unity表面上是DirectX做的，坐标也是左手坐标系，但如果你用camera的API：Camera.projectionMatrix来得到投影矩阵，抱歉，他给你OpenGL的。
为了和龙书保持一致，我自己写了一个构造DirectX投影矩阵的方法：
public static Matrix4x4 GetProjectionMatrixDX(this Camera camera)
{
float r = camera.aspect;
float a = camera.fieldOfView * Mathf.Deg2Rad;
float f = camera.farClipPlane;
float n = camera.nearClipPlane;
Matrix4x4 Out = new Matrix4x4();
Out.m00 = 1 / (r * Mathf.Tan(a / 2));
Out.m01 = 0;
Out.m02 = 0;
Out.m03 = 0;
Out.m10 = 0;
Out.m11 = 1 / Mathf.Tan(a / 2);
Out.m12 = 0;
Out.m13 = 0;
Out.m20 = 0;
Out.m21 = 0;
Out.m22 = f / (f - n);
Out.m23 = (f * n) / (n - f);
Out.m30 = 0;
Out.m31 = 0;
Out.m32 = 1;
Out.m33 = 0;
return Out;
}
以及view矩阵在z轴基向量也与DirectX不一致，不过差别不大，只有z轴基向量和位移量相反，构造之：
public static Matrix4x4 GetViewMatrixDX(this Camera camera)
{
Matrix4x4 Out = camera.worldToCameraMatrix;
Out.m20 = -Out.m20;
Out.m21 = -Out.m21;
Out.m22 = -Out.m22;
Out.m23 = -Out.m23;
return Out;
}
坐标推导
这里只要看就行了，我用Python的Sympy和Jupyter推导一下坐标变化。
假如我们在View空间有一点p=[x,y,z]，我们经过投影变换：
经过齐次除法到NDC：
这个看起来特别奇怪的点就是我们要自行设定的NDC点；我们现在要来一个逆过程，但是很明显，我们很难找出z来一个“齐次乘法”，不过我看DirectX标准库的实现很有意思，它先逆了投影变换的过程(既乘上投影矩阵的逆矩阵)，我们用代码推导一下：
实际上是sympy库没用整理好，结果是[x/z, y/z, 1, 1/z]，这一下就豁然开朗了。
视锥体的构造方法（继续）
回过来看代码：
public static BoundingFrustum CreateFromCamera(Camera camera)
{
BoundingFrustum Out = new BoundingFrustum();
//首先构建NDC空间下的视锥体，是一个[-1,1]x[-1,1]x[0,1]的盒子
Vector4[] HomogenousPoints = new Vector4[6];
HomogenousPoints[0] = new Vector4(1, 0, 1, 1);//右（在远平面）
HomogenousPoints[1] = new Vector4(-1, 0, 1, 1);//左
HomogenousPoints[2] = new Vector4(0, 1, 1, 1);//上
HomogenousPoints[3] = new Vector4(0, -1, 1, 1);//下
HomogenousPoints[4] = new Vector4(0, 0, 0, 1);//近平面
HomogenousPoints[5] = new Vector4(0, 0, 1, 1);//远平面
Matrix4x4 matInverse = camera.GetProjectionMatrixDX().inverse;//投影矩阵的逆矩阵
//将ndc空间的各点计算到view空间下
Vector4[] Points = new Vector4[6];
for(int i = 0; i < 6; ++i)
{
Points[i] = matInverse * HomogenousPoints[i];
}
//由于view->proj->齐次除法->ndc间还有齐次除法，要把齐次除法过程逆转
//ndc * projInv的结果是[x/z, y/z, 1, 1/z]，前两个刚好是斜率
Out.RightSlope = Points[0].x;
Out.LeftSlope = Points[1].x;
Out.TopSlope = Points[2].y;
Out.BottomSlope = Points[3].y;
Out.Near = (Points[4] / Points[4].w).z;
Out.Far = (Points[5] / Points[5].w).z;
Out.viewMat = camera.GetViewMatrixDX();
return Out;
}
前面的理论都搞清除，这部分代码应该就没什么问题了。
碰撞检测
正常人想到的方法，将包围盒每个顶点往视锥体里过一遍，假如有一个顶点在视锥体内部，就说明两者相交。
考虑到特殊情况，包围盒把视锥体包住了，那就检测不到了，还要特殊逻辑处理？
不过我们不用这种方法，而采取另一种方式。
平面表示
我们不要把视锥体看做8个顶点构成的形状，而是6个无限大的平面，现在我们来试着表示一下平面。
我们用数学上称之为一般式的表达方式：A*x+B*y+C*z+D=0。
这种方式有诸多好处：
- 首先我们只要记录v=[A, B, C, D]这四维向量即可表达一个平面(又记做v=(n,d))。
- 其次，[A, B, C]就是这个平面的法向量N（不保证标准化）；
- 其三，D的绝对值是平面到原点的最短距离。
- 其四，当空间中有一点p=[x, y, z]，我们只要将其拓展为齐次向量p=[x,y,z,1]，然后将其与v=(n,d)点乘（保证v=(n,d)中法向量n经过标准化），既v·p = d'，|d'|为点到平面的距离，并且如果d'>0，则点在平面的正半空间，否则在负半空间
回到视锥体，我们将其看做6个平面包裹住的空间，其法向量朝着内部，通过前面的参数，我们可以轻松构造6个平面：
Vector4[] Planes = new Vector4[6];
Planes[0] = new Vector4(0, 0, 1, -Near);
Planes[1] = new Vector4(0, 0, -1, Far);
Planes[2] = new Vector4(-1, 0, RightSlope, 0);
Planes[3] = new Vector4(1, 0, -LeftSlope, 0);
Planes[4] = new Vector4(0, -1, TopSlope, 0);
Planes[5] = new Vector4(0, 1, -BottomSlope, 0);
不信可以看一看，这些数值刚好符合标准式，不过要注意后四个平面还没有进行标准化。
平面的变换
我们常常用矩阵对点和向量进行变换，但对平面变换几乎没有过吧？
对平面的变换在这里指直接对一般式进行平移和旋转，不过很遗憾，直接用矩阵恐怕是不行。
如之前所说，一般式由v=(n,d)组成，对向量的操作我们很娴熟了，直接n=(n, 0)再n=Mn即可。
对于d，本身就是面距离原点最近的距离，从原点到此点的方向向量就是法向量n，我们得到这个最近点:d=normal * (-n.z)(注意，n是没变换之前的)，拓展至齐次向量d=(d,1)，将其旋转，旋转后这个点依然在平面上，
那么算这个变换后的新点到原点的距离？不不不，原来是最近的点，变换后就不一定了
我们原来的d点变为了d'，很显然，最近点是d''点，那么问题就转为了已知平面法线n'=(nx, ny, nz)和平面上一点p=(a,b,c)，求平面的一般式方程。推导：设平面上任意一点p0=(x,y,z)，向量p0-p与平面平行，则向量p0-p与法线n'垂直，既dot(n', (p0-p))=0.
得到nx*(x-a) + ny*(y-b) + nz*(z-c)=0，整理得nx*x+ny*y+nz*z-(nx*a+ny*b+nz*c)=0既一般式。
可以看到，这个新的距离d''正是-dot(d', n')
由此写出函数：
Vector4 TransformPlane(Vector4 Plane, Matrix4x4 M)
{
Vector4 normal = new Vector4(Plane.x, Plane.y, Plane.z, 0);
normal = normal.normalized;
Vector4 dVector = normal * -Plane.w;
dVector.w = 1;
Vector4 newN = M * normal;
Vector4 newD = M * dVector;
float d = -Vector3.Dot(newD, newN);
Vector4 Out = new Vector4(newN.x, newN.y, newN.z, d);
return Out;
}
注意一点，如果变换还有缩放，那个法向量变换所需的矩阵要经过逆转置操作，既n = transpose(inverse(M))*n
n=((M-1)T)n，不过view矩阵没有缩放操作，所以这里就不实现了。
AABB包围盒碰撞
之前说我们为了减小复杂度才把ndc转为view的，现在又要旋转和移动，这不矛盾了吗？
之前那么说是为了想像view空间的样子，但如果我们真的在view空间检测了，那么原来还和轴对齐的AABB包围盒肯定就不和轴对齐了，那么算法复杂度就会提升很多，因此为了迁就包围盒，还是至少将视锥体平面转换到世界空间下，就用view矩阵的逆矩阵就可以。
检测算法：每个平面都有正半空间和负半空间，假如存在一个平面，包围盒完全在平面的负半空间中，就说明平面在视锥体之外。如果不存在这样的平面，那就说明两者相交。
复述一遍，包围住视锥体的所有平面法向量都向内。
看上图，绿色立方体在视锥体之外，因为绿色立方体在上平面的负半空间，满足至少存在一个平面的条件，因此立方体在视锥体之外。
那么如何判断立方体在负半空间呢？看下图：
- 图a，立方体在平面正半空间，视锥体不一定与包围盒相交（情况类似上图的立方体和视锥体右平面）。
- 图b，立方体完全在平面负半空间，满足条件，视锥体和包围盒一定不相交
- 图c，立方体跨越正半空间和负半空间，两者一定相交？不不不，想像视锥体在很高很高，而包围盒很低，你从上面看，以为两者相交，但其实并不相交，如下图：
关键点在于找到上图所示的点P和点Q，向量PQ是和平面法向量n方向最接近的对角线向量。什么是最接近？假如向量n在x轴是从小到大，那么PQ向量也一样是从小到大，同理y、z轴也一样。
所以得到求法：
// 分别找x, y, z轴
for(int j = 0; j < 3; ++j)
{
//如果法向量当前轴是增大方向
if( planeNormal[j] >= 0.0f )
{//则向量PQ也是增大
P[j] = box.minPt[j];
Q[j] = box.maxPt[j];
}
else
{//否则PQ在当前轴向减小
P[j] = box.maxPt[j];
Q[j] = box.minPt[j];
}
这样就得到P和Q点了，那么怎么知道点在正半空间还是负半空间呢，这就用到之前的原理，将平面四维向量点坐标，结果的正负就关系到正负空间。
综上所述，写出如下代码：
public bool Intersects(Bounds bound)
{
//构建视锥体平面，用一般式（Ax+By+Cz+D=0）的各项系数表示，其中[A, B, C]可看做未标准化的法向量
Vector4[] Planes = new Vector4[6];
Planes[0] = new Vector4(0, 0, 1, -Near);
Planes[1] = new Vector4(0, 0, -1, Far);
Planes[2] = new Vector4(-1, 0, RightSlope, 0);
Planes[3] = new Vector4(1, 0, -LeftSlope, 0);
Planes[4] = new Vector4(0, -1, TopSlope, 0);
Planes[5] = new Vector4(0, 1, -BottomSlope, 0);
//包围盒在world空间下的最大、最小点
Vector3 boxMin_World = bound.center - bound.extents;
Vector3 boxMax_World = (bound.center + bound.extents);
for (int i = 0; i < 6; ++i)
{
//将平面转换到world空间
Vector4 Plane = TransformPlane(Planes[i], viewMat.inverse, i);
//计算相对于当前平面时，包围盒的PQ点
Vector4 P = new Vector4(0, 0, 0, 1);
Vector4 Q = new Vector4(0, 0, 0, 1);
for (int j = 0; j < 3; ++j)
{
if (Plane[j] >= 0)
{
P[j] = boxMin_World[j];
Q[j] = boxMax_World[j];
}
else
{
P[j] = boxMax_World[j];
Q[j] = boxMin_World[j];
}
}
//P、Q距离平面的距离
float P_Dist = Vector4.Dot(Plane, P);
float Q_Dist = Vector4.Dot(Plane, Q);
//如果P和Q都在平面负半空间，则包围盒在视锥体外面
if (P_Dist <= 0 && Q_Dist <= 0)
{
return false;
}
}
//全都判断完，不存在平面将包围盒置于负半空间，说明相交
return true;
}
结果展示
于此，视锥体碰撞完成。
https://imgchr.com/i/Gl63sx
中间出错很多次，我踩了很多坑，所以图上能看到不少调试的痕迹。
黄色的是视锥体；红线是射线，原点是当前平面离原点最近的点，方向是法线方向；白线是最近点和原点的连线。
如果感觉：不对啊，我看到红线都悬空了，还在视锥体外！这是正确的结果，平面离原点的最近点不一定在视锥体内。
给每一个物体挂一个包围盒，每帧判断，如果碰撞检测失败，就将渲染队列设为0，就将其剪裁了。
不过对所有物体都碰撞检测还是有些蠢了，实际上还要加上场景管理，比如我们加点特效：
https://imgchr.com/i/Gl68L6
嗯，就是知名的四叉树场景管理，开局自动根据物体划分四叉树，这样很容易就能剔除不需要渲染的物体（不过我东西少很显然不划算）。
现在还有些bug，就是如上图，激活很及时，但不能及时让激活的四叉树节点关闭，等修完BUG再把代码放上来。
