美文网首页
【Frustum Culling】视锥体剪裁数学原理和代码实现

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

作者: crossous | 来源:发表于2020-03-31 17:52 被阅读0次

    前言

      剪裁是渲染中常用的手段,避免将渲染资源浪费在无意义的片段中,在渲染管线的齐次除法,渲染管线就会帮我们做一次剪裁,防止在视锥体外的顶点跑到像素着色器被渲染。
      但这终究会进入渲染管线,会进入顶点着色器乃至曲面细分和集合着色器,因此我们最好在进入渲染管线前就先做一次粗糙一些的剪裁,在渲染管线阶段之前的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方向最接近的对角线向量。什么是最接近?假如向量nx轴是从小到大,那么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];
        }
    

      这样就得到PQ点了,那么怎么知道点在正半空间还是负半空间呢,这就用到之前的原理,将平面四维向量点坐标,结果的正负就关系到正负空间。
      综上所述,写出如下代码:

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

    相关文章

      网友评论

          本文标题:【Frustum Culling】视锥体剪裁数学原理和代码实现

          本文链接:https://www.haomeiwen.com/subject/emjwuhtx.html