视锥体平面公式推导

龙犁

理论

最近在使用Vulkan做一个GPU Frustum Cull。需要完成Frustum的上下左右远近六个平面和一个球形包围盒的相交测试。其中就需要快速的去构建Frustum的六个平面，在完成Forward Plus Render的操作当中。同样需要用到这个操作。在这里系统的总结一下，加强理解。
在这里我们是通过投影矩阵(Project Matrix)来推导视锥体六个平面方程。而且使用的矩阵不同，推导出的对应的视锥体平面也会在不同的坐标系当中。只使用投影矩阵推导出相应的视锥体平面是在视线空间(View Space)内，使用View-Project矩阵的话，推导出的平面在世界空间(World Space)中。使用World-View-Projection推导出的平面公式就是在对象空间当中(Object Space)。我们这里只使用投影矩阵。那就意味这这意味着摄像机位于世界坐标系的原点，而我们沿着正Z轴观察。
假设顶点为 V(x,y,z,1) 。 M_{(ij)} 是我们的投影矩阵是一个4X4矩阵。 V'=(x',y',z',w') 是顶点和投影矩阵相乘的结果。具体的公式如下所示: 其中的T符号代表其转置矩阵。
\mathbf{v}^{\prime}=\mathbf{v M} \Rightarrow\left(\begin{array}{ll}x^{\prime} , y^{\prime} ,z^{\prime} ,w^{\prime}\end{array}\right)=\left(\begin{array}{l}x \cdot m_{11}+y \cdot m_{21}+z \cdot m_{31}+w \cdot m_{41} \\x \cdot m_{12}+y \cdot m_{22}+z \cdot m_{32}+w \cdot m_{42} \\x \cdot m_{13}+y \cdot m_{23}+z \cdot m_{33}+w \cdot m_{43} \\x \cdot m_{14}+y \cdot m_{24}+z \cdot m_{34}+w \cdot m_{44}\end{array}\right)^{T}=\left(\begin{array}{c}\mathbf{v} \bullet \operatorname{col}_{1} \\\mathbf{v} \bullet \operatorname{col}_{2} \\\mathbf{v}  \bullet \operatorname{col}_{3} \\\mathbf{v} \bullet \operatorname{col}_{4}\end{array}\right)^{T}
其中 \bullet 表示点乘， \operatorname{col}_{j}=\left(\begin{array}{llll}m_{1 j},m_{2 j}, m_{3 j}, m_{4 j}\end{array}\right) 表示矩阵M的第j列所代表的向量。在变换之后顶点 v' 处于剪裁空间。在这个空间里，观察范围实际上是是一个轴对齐的盒子，其大小是各种不同图形API特定的。如果顶点 v' 在这个空间里面，那么未经变换的顶点 v 就在在 "未变换的 "坐标内。顶点 v 在这个坐标系内，那么 v' 的各个分量都要满足以下不等式:
\qquad\qquad\qquad\qquad \qquad\qquad\qquad\begin{array}{r}-w^{\prime}<x^{\prime}<w^{\prime} \\-w^{\prime}<y^{\prime}<w^{\prime} \\0<z^{\prime}<w^{\prime}\end{array}
得出的结论如下表所示：


接下来就开始我们的推导！
首先如果我们的满足我们的第一个式子 -w'<x' 接下来就用这个来进行推导。从上面最开始的式子可以推导下面式子
\qquad\qquad\qquad\qquad\qquad \qquad \qquad\qquad-(v \bullet col_{1}) < (v \bullet col_{4})
可以化简为下式：
\qquad\qquad\qquad\qquad\qquad \qquad \qquad\qquad 0 < (v \bullet col_{1})+(v \bullet col_{4})
最后得出下式：
\qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad 0 < (v \bullet (col_{1}+ col_{4}))
那么就可以推导下面这个式子：
\qquad \qquad x\left(m_{14}+m_{11}\right)+y\left(m_{24}+m_{21}\right)+z\left(m_{34}+m_{31}\right)+w\left(m_{44}+m_{41}\right)=0
这已经表示了在视线空间的视锥体左侧剪裁平面的平面方程。并且假如$w=1$的情况下。可得下式：
\qquad \qquad x\left(m_{14}+m_{11}\right)+y\left(m_{24}+m_{21}\right)+z\left(m_{34}+m_{31}\right)+\left(m_{44}+m_{41}\right)=0
这就相当于一个平面方程，如下所示：
\qquad \qquad\qquad  \qquad \qquad \qquad\qquad a x+b y+c z+d=0
其中的a,b,c参数如下
\qquad \qquad a=m_{14}+m_{11}, \quad b=m_{24}+m_{21}, \quad c=m_{34}+m_{31}, \quad d=m_{44}+m_{41}
到这里就可以提出一个结论，视锥体左剪裁平面可以直接从投影矩阵直接提取。需要注意的是，所得到的平面方程没有被归一化(即平面的法向量不是单位向量)，而且法向量是指向内侧空间的。这意味着顶点满足 0 <ax+ by+ cz +d 这个条件的话，则该顶点在左剪裁平面的内侧空间(也就是在视锥体当中)。
其他的各个平面方程就不一一推导。可以看下面这个表。有完整的六个平面的方程。


实际使用

视锥体和球型包围盒相交测试

在上面我们可以得出相应的理论基础。接下来我们就是要具体的运用了。在这里我们需要去视锥体是否和一个球形包围盒相交。
glm::vec4 normalizePlane(glm::vec4 p)
{
    return p / glm::length(glm::vec3(p));
}

glm::mat4 projection = params.projmat;
glm::mat4 projectionT = transpose(projection);
glm::vec4 frustumX = normalizePlane(projectionT[3] + projectionT[0]); // x + w < 0
glm::vec4 frustumY = normalizePlane(projectionT[3] + projectionT[1]); // y + w < 0

cullData.frustum[0] = frustumX.x;
cullData.frustum[1] = frustumX.z;
cullData.frustum[2] = frustumY.y;
cullData.frustum[3] = frustumY.z;
首先我们拿到投影矩阵的转置矩阵。之后projectionT[3] + projectionT[0]这里也就是在上面提到的 (col_{1}+ col_{4}) 参数。projectionT[3] + projectionT[0]对应的则是 (col_{2}+ col_{4}) 参数。分别用于左右和上下四个平面。之后将其对应的分量塞入一个数组中。后续传入Shader供计算使用。
vec3 center = sphereBounds.xyz;
// 将该center转到视线空间下
center = (cullData.view * vec4(center,1.f)).xyz;
float radius = sphereBounds.w;

// the left/top/right/bottom plane culling utilizes frustum symmetry to cull against two planes at the same time
visible = visible && center.z * cullData.frustum[1] - abs(center.x) * cullData.frustum[0] > -radius;
visible = visible && center.z * cullData.frustum[3] - abs(center.y) * cullData.frustum[2] > -radius;

if(cullData.distCull != 0)
{   // the near/far plane culling uses camera space Z directly
    // 是否满足在远近平面之间
    visible = visible && center.z + radius > cullData.znear && center.z - radius < cullData.zfar;
}
在这里我们首先将这个球形的圆心转换到了视线空间当中。因为在这里我们只使用的是投影矩阵来构建视锥体。首先来看看关于左、右平面的检测。center.z * cullData.frustum[1] - abs(center.x) * cullData.frustum[0] > -radius在这里有一个trick。在这里是假设了咱们的y分量在此刻是0。这样就可以减少相应的计算(这就是在Shader中的优化)。对于上、下平面的检测也是同样假设x分量是0。并且我们还需要注意这个abs()函数，这同样是一个trick。这使得可以通过这一行代码同时完成相对称平面的检测。在这里假如是左平面的检测，并且没有使用abs()函数的话。那么这里是不符合 (col_{1}+ col_{4}) 的要求。但是又因为在检测左平面的时候在这里的假设是所有的X分量都是小于0。通过这个abs()函数可以重新符合对应的公式。对于右平面来说X分量就是大于0的。那么abs()函数也就无影响。同样也符合右平面方程。在这里就还是减少了相应的计算量。对于上、下平面也是同样的操作。最后关于远近平面的检测应该就不需要多说了。
Forward+ Light Culling

static float2 screen_to_view_at_z1(constant AAPLFrameData & frameData, ushort2 screen)
{
    const float3 screenToViewSpace = frameData.screenToViewSpace;
    return float2(screen) * float2(screenToViewSpace.x, -screenToViewSpace.x) + float2(screenToViewSpace.y, -screenToViewSpace.z);
}

// Unproject depth from screen space to view space, where the culling is done.
// 将最大最小深度转化到视线空间当中
float minDepthView = unproject_depth(frameData, tile_data->minDepth);
float maxDepthView = unproject_depth(frameData, tile_data->maxDepth);
// 分别获取到左上和右下两个点在视线空间当中
float2 minTileViewAtZ1 = screen_to_view_at_z1(frameData, threadgroup_id * threadgroup_size);
float2 maxTileViewAtZ1 = screen_to_view_at_z1(frameData, (threadgroup_id + 1) * threadgroup_size);

// Calculate the normals of the tile bounding planes.
// 求六个平面组成的视锥体，由于左右上下平面都过原点，所以offset都是零。
// 但是远近平面是不相交于原点，所以有相应的offset
// 并且分别拿到tile的左上点以及右下点。我们在这里是去获取这个平面的法线。
// 我们要清楚我们在这个平面上任何一个点和这个法线点乘都是零。并且通过获取到这个左上和右下这两个点。计算出相应的法线。
AAPLPlane tile_planes[6] = {
    { normalize(float3(1.0, 0.0, -maxTileViewAtZ1.x)), 0.0f }, // right
    { normalize(float3(0.0, 1.0, -minTileViewAtZ1.y)), 0.0f }, // top
    { normalize(float3(-1.0, 0.0, minTileViewAtZ1.x)), 0.0f }, // left
    { normalize(float3(0.0, -1.0, maxTileViewAtZ1.y)), 0.0f }, // bottom
    { float3(0.0, 0.0, -1.0), -minDepthView },                 // near
    { float3(0.0, 0.0, 1.0), maxDepthView }                    // far
};
在Forward+的实现中，同样需要做一个关于球形光源和每个Tile的构建出的视锥体做一个相交测试。其他的就不多赘述。主要就是讲讲这个特殊的视锥体构建。在这里我们需要根据一个Tile去构建相应的视锥体。首先我们可以获取到相应的Tile编号。然后将其从屏幕空间当中转化到视线空间当中。在这里的远近平面是根据我们的最大和最小深度来决定的。同时我们这里还有一个特殊的点。在这里取得是这个Tile上左上和右下两个点。我们假设所有的点都在一个深度为1的平面。这样可以方便来构建我们的视锥体平面并且减少计算量。首先我们可以看到右平面。在这里我们构建的这个平面的法线。同样可以达到效果。只是方式不同。我们可以通过这样的计算得出。首先Tile 右下这个点，就是在这个右平面之上的。并且就可以通过 (maxTileViewatZ1.x,0,1) \bullet normal = 0 。这个式子推导出相应的法线。对于其他平面也是如此。
References

中央空调新风净水

佩服佩服！