【Frustum Culling】视锥体剪裁数学原理和代码实现

【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().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再把代码放上来。

相关新闻

倩女幽魂PK值要怎么减少?
365bet手机投注网

倩女幽魂PK值要怎么减少?

🕒 07-14 👽 2927
胼胝体主要功能
365bet手机投注网

胼胝体主要功能

🕒 08-12 👽 3441
家常醋熘白菜
365bet手机投注网

家常醋熘白菜

🕒 09-03 👽 6429
上海显高装饰设计工程有限公司
365bet手机投注网

上海显高装饰设计工程有限公司

🕒 08-12 👽 8747
超越小嶋阳菜!谁是AKB48团体里的第一美胸?
Bet体育365提款要多久2022

超越小嶋阳菜!谁是AKB48团体里的第一美胸?

🕒 01-01 👽 8993
微信租车平台,微信汽车租赁app有哪些?
365体育投注一直进不去

微信租车平台,微信汽车租赁app有哪些?

🕒 01-18 👽 6881