games101 作業6

這是games101 現代計算機圖形學 作業06 by coolwx

感謝https://blog.csdn.net/miyu1994/article/details/107006010/ 這位大佬的文章，因為我一開始真的寫不出來SAH（啟發式搜索），看了這位大神的算法，我最終才理解了他的算法，其實還是比較簡單的，關於SAH，我應該會寫一篇筆記，來記錄我的感想，再次感謝以上這位大神的算法，因為我一開始真的想不出這個啟發式算法，所以代碼我會使用他的代碼，然后添上我的注釋，便於理解（哎，還是太菜，沒寫出來）。

作業要求：

在之前的編程練習中，我們實現了基礎的光線追蹤算法，具體而言是光線傳

輸、光線與三角形求交。我們采用了這樣的方法尋找光線與場景的交點：遍歷場景

中的所有物體，判斷光線是否與它相交。在場景中的物體數量不大時，該做法可以

取得良好的結果，但當物體數量增多、模型變得更加復雜，該做法將會變得非常低

效。因此，我們需要加速結構來加速求交過程。在本次練習中，我們重點關注物體

划分算法 Bounding Volume Hierarchy (BVH)。本練習要求你實現 Ray-Bounding

Volume 求交與 BVH 查找。

首先，你需要從上一次編程練習中引用以下函數：

• Render() in Renderer.cpp: 將你的光線生成過程粘貼到此處，並且按照新框

架更新相應調用的格式。

• Triangle::getIntersection in Triangle.hpp: 將你的光線-三角形相交函數

粘貼到此處，並且按照新框架更新相應相交信息的格式。

在本次編程練習中，你需要實現以下函數：

• IntersectP(const Ray& ray, const Vector3f& invDir,

const std::array<int, 3>& dirIsNeg) in the Bounds3.hpp: 這個函數的

作用是判斷包圍盒 BoundingBox 與光線是否相交，你需要按照課程介紹的算

法實現求交過程。

• getIntersection(BVHBuildNode* node, const Ray ray)in BVH.cpp: 建

立 BVH 之后，我們可以用它加速求交過程。該過程遞歸進行，你將在其中調

用你實現的 Bounds3::IntersectP.

 

 

其他，這次作業里面有兩個坑。

一個是castRay方法現在在Scene.cpp里面，然后castRay函數發生了一點變化，ray作為參數會傳進去，其實ray也就定義了direction和origin兩個參數。

細讀這次的所有代碼，也算是對整個光線追蹤有了更深刻的理解吧，可惜我本人還是太菜了。

第二個坑是我一開始沒畫出來我佛了，后來也是參考了這位大神的代碼，發現我對於Bounds包圍盒中的求解中，沒有考慮方向中有負數的情況，直接佛了。。

 

首先進入Triangle.h

修改getInsection函數，需要注意的是，這里需要理解以下參數

coords：光線經過的那一點的坐標，這里使用ray（double）（注意，我們在Ray.h中重載了（），用於計算t時刻的光線坐標）

happend：光線和三角形是否相交？

m：材質

normal：交點處三角形的法線

distance：光線經過的時間，也就是我們算出來的t_tmp。

object：事實上你進入Objects.h,就會發現基類Object（接口，里面一堆純虛函數），這里作為一個指針，指向本對象，也就是this（框架的編程還是很漂亮的），指代我們的三角形。

（注意我們作業5中我自己寫的求交，感覺low爆了，還是框架的寫的好看）

 

inline Intersection Triangle::getIntersection(Ray ray)
{
    Intersection inter;

    if (dotProduct(ray.direction, normal) > 0)
        return inter;
    double u, v, t_tmp = 0;
    Vector3f pvec = crossProduct(ray.direction, e2);//s1
    double det = dotProduct(e1, pvec);
    if (fabs(det) < EPSILON)
        return inter;//no solution

    double det_inv = 1. / det;
    Vector3f tvec = ray.origin - v0;//s
    u = dotProduct(tvec, pvec) * det_inv;
    if (u < 0 || u > 1)
        return inter;
    Vector3f qvec = crossProduct(tvec, e1);
    v = dotProduct(ray.direction, qvec) * det_inv;
    if (v < 0 || u + v > 1)
        return inter;
    t_tmp = dotProduct(e2, qvec) * det_inv;

    // TODO find ray triangle intersection
    if(t_tmp<0)
        return inter;
    inter.distance = t_tmp;
    inter.happened = true;
    inter.m = m;
    inter.obj = this;
    inter.normal = normal;
    inter.coords = ray(t_tmp);
    return inter;
}

　　

然后進入Renderer.cpp，修改Render函數，注意這里，我們上面兩行代碼計算x和y就是作業5的答案，然后下面對比作業5的代碼出現了變化，因為框架修改了

void Renderer::Render(const Scene& scene)
{
    std::vector<Vector3f> framebuffer(scene.width * scene.height);

    float scale = tan(deg2rad(scene.fov * 0.5));
    float imageAspectRatio = scene.width / (float)scene.height;
    Vector3f eye_pos(-1, 5, 10);
    int m = 0;
    for (uint32_t j = 0; j < scene.height; ++j) {
        for (uint32_t i = 0; i < scene.width; ++i) {
            // generate primary ray direction
            float x = (2 * (i + 0.5) / (float)scene.width - 1) *
                      imageAspectRatio * scale;
            float y = (1 - 2 * (j + 0.5) / (float)scene.height) * scale;
            // TODO: Find the x and y positions of the current pixel to get the
            // direction
            //  vector that passes through it.
            // Also, don't forget to multiply both of them with the variable
            // *scale*, and x (horizontal) variable with the *imageAspectRatio*
            Vector3f dir = Vector3f(x, y, -1); 
            dir = normalize(dir);
            Ray ray(eye_pos,dir);
            framebuffer[m++] = scene.castRay(ray,0);
            // Don't forget to normalize this direction!

        }
        UpdateProgress(j / (float)scene.height);
    }
    UpdateProgress(1.f);

    // save framebuffer to file
    FILE* fp = fopen("binary.ppm", "wb");
    (void)fprintf(fp, "P6\n%d %d\n255\n", scene.width, scene.height);
    for (auto i = 0; i < scene.height * scene.width; ++i) {
        static unsigned char color[3];
        color[0] = (unsigned char)(255 * clamp(0, 1, framebuffer[i].x));
        color[1] = (unsigned char)(255 * clamp(0, 1, framebuffer[i].y));
        color[2] = (unsigned char)(255 * clamp(0, 1, framebuffer[i].z));
        fwrite(color, 1, 3, fp);
    }
    fclose(fp);    
}

 

進入Bounds3.cpp ，修改IntersectP函數，此函數用於求包圍盒是否和光線相交，這里特別要注意負數的情況，一開始我想后面這個數組有什么用，然后一直沒用到，最后卡住了，佛了，其實這個就是判斷大小用的，算最大最小的時候注意負數的情況

inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,
                                const std::array<int, 3>& dirIsNeg) const
{
    float t_Min_x = (pMin.x - ray.origin.x)*invDir[0];
    float t_Min_y = (pMin.y - ray.origin.y)*invDir[1];
    float t_Min_z = (pMin.z - ray.origin.z)*invDir[2];
    float t_Max_x = (pMax.x - ray.origin.x)*invDir[0];
    float t_Max_y = (pMax.y - ray.origin.y)*invDir[1];
    float t_Max_z = (pMax.z - ray.origin.z)*invDir[2];
    if(!dirIsNeg[0])
    {
        float t = t_Min_x;
        t_Min_x = t_Max_x;
        t_Max_x = t;
    }
    if(!dirIsNeg[1])
    {
        float t = t_Min_y;
        t_Min_y = t_Max_y;
        t_Max_y = t;
    }
    if(!dirIsNeg[2])
    {
        float t = t_Min_z;
        t_Min_z = t_Max_z;
        t_Max_z = t;
    }

    float t_enter = std::max(t_Min_x,std::max(t_Min_y,t_Min_z));
    float t_exit =  std::min(t_Max_x,std::min(t_Max_y,t_Max_z));
    if(t_enter<t_exit&&t_exit>=0)
    return true;
    else
    return false;
    
    // invDir: ray direction(x,y,z), invDir=(1.0/x,1.0/y,1.0/z), use this because Multiply is faster that Division
    // dirIsNeg: ray direction(x,y,z), dirIsNeg=[int(x>0),int(y>0),int(z>0)], use this to simplify your logic
    // TODO test if ray bound intersects
    
}

　　最后進入BVH.cpp,修改我們的getInsection,注意，我們這里是為了求交，根據閆神所講的，肯定要到葉子節點才算求交，而中間如果和包圍盒都沒交點，那么和三角形一定無交點，簡單二叉樹遍歷就行。

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
    Intersection intersect;
    Vector3f invdir(1./ray.direction.x,1./ray.direction.y,1./ray.direction.z);
    std::array<int, 3> dirIsNeg;
    dirIsNeg[0] = ray.direction.x>0;
    dirIsNeg[1] = ray.direction.y>0;
    dirIsNeg[2] = ray.direction.z>0;
     if(!node->bounds.IntersectP(ray,invdir,dirIsNeg))
     {
         return intersect;
     }
    if(node->left == nullptr && node->right==nullptr)
    {
       return node->object->getIntersection(ray);
    }
    Intersection h1 = getIntersection(node->left,ray);
    Intersection h2 = getIntersection(node->right,ray);
    return h1.distance<h2.distance?h1:h2;
    return intersect;

    
    // TODO Traverse the BVH to find intersection

}

 

其實到這里，我們的任務已經完成了，效果如下：

 

 

 

但是其實還沒完，注意我們的作業提高要求：

SAH 查找：自學 SAH(Surface Area Heuristic) , 正 確實現 SAH 加速，並且提交結果圖片，並在 README.md 中說明 SVH 的實現 方法，並對比 BVH、SVH 的時間開銷。(可參考 http://15462.courses.cs .cmu.edu/fall2015/lecture/acceleration/slide_024，也可以查找其他資 料)

 

點進去ppt，發現是這個東西

 

 

告訴我們是一個啟發性算法，這里的思想是如果有兩個包圍盒，B包圍A，那么光線過B也過A的概率就是SA/SB，

注意，我們在這一次作業的時候，觀察前面的rebuild方法，是遞歸地將objects數組分成左右兩個部分，那么SA可以說代表左邊的objects形成的總包圍盒的面積，而SB可以說是代表右邊的划分總的包圍盒的面積，SN代表我們的最大的包圍盒的面積，我們的想法是啟發式搜索（哇本人看到啟發式搜索就頭疼），這里Ctrav表示當前划分一次的花費，這里取1/8（這個值應該是可以版的，具體看https://blog.csdn.net/qq_39300235/article/details/106999514，這篇文章說的很好）

我們的目標是分割這兩個數組（搜索），所花費的cost越少越好。（注意，我們一開始划分數組的時候，方法是就是取最大的跨度維度，然后從中間划分）

 

但是啟發式搜索目的是不要讓我們從中間划分，而是划分的時候讓cost（也就是上面的式子C越小越好），找出這種划分方法就行。

 

https://blog.csdn.net/qq_39300235/article/details/106999514這篇文章其實很清楚，但是里面一堆桶的說法很麻煩，記住，我們的目的就是找出某個划分方法，然后讓cost最小，這樣形成的BVH結構不會有很大的深度。

 

我截取這篇文章的核心思想就是在每一個維度（x，y，z）分別計算一種划分（左，右），目的是求左右之間的中間值

 

這道題坑的地方是其實他 在bound3s中寫了一個offset（），這個東西是就是用來求上面那篇文章的桶的標號的。

還有，這里面有一個bid = bid -1 ，其實是我們分了12個桶（數組最大數11），但是根據offset是可以算出12的（12棵樹中間的間隔其實有13個，那么就把12的地方換成11就行，也就是桶的大小最后一個桶會大一點）

 

感謝https://blog.csdn.net/miyu1994/article/details/107006010/大神的代碼，我這里直接用了他的代碼，不過標記了一點注釋，這樣比較容易理解。

 

BVHBuildNode* BVHAccel::recursiveBuild(std::vector<Object*> objects)
{
    BVHBuildNode* node = new BVHBuildNode();

    // Compute bounds of all primitives in BVH node
    Bounds3 bounds;
    for (int i = 0; i < objects.size(); ++i)
        bounds = Union(bounds, objects[i]->getBounds());
    if (objects.size() == 1) {
        // Create leaf _BVHBuildNode_
        node->bounds = objects[0]->getBounds();
        node->object = objects[0];
        node->left = nullptr;
        node->right = nullptr;
        return node;
    }
    else if (objects.size() == 2) {
        node->left = recursiveBuild(std::vector{objects[0]});
        node->right = recursiveBuild(std::vector{objects[1]});

        node->bounds = Union(node->left->bounds, node->right->bounds);
        return node;
    }
    else 
    {
        Bounds3 centroidBounds;
        //算出最大的包圍盒（通用的，因為兩個方法都要用到） 
        for (int i = 0; i < objects.size(); ++i)
            centroidBounds =
                Union(centroidBounds, objects[i]->getBounds().Centroid());

        std::vector<Object*> leftshapes;
        std::vector<Object*> rightshapes;

            switch (splitMethod)//這里注意了在BVH.h里面有個枚舉類，構造函數中的初始將決定是普通方法，還是SAH 
            {
            case SplitMethod::NAIVE:
            {
                int dim = centroidBounds.maxExtent();//算出最大的跨度對應的值，x為0，y為1，z為2 
                int index = objects.size() / 2;
                switch (dim)
                //排序，針對最大的跨度排序 
                {
                case 0:
                    std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                        return f1->getBounds().Centroid().x <
                               f2->getBounds().Centroid().x;
                    });
                    break;
                case 1:
                    std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                        return f1->getBounds().Centroid().y <
                               f2->getBounds().Centroid().y;
                    });
                    break;
                case 2:
                    std::sort(objects.begin(), objects.end(), [](auto f1, auto f2) {
                        return f1->getBounds().Centroid().z <
                               f2->getBounds().Centroid().z;
                    });
                    break;
                }

                auto beginning = objects.begin();
                auto middling = objects.begin() + index;
                auto ending = objects.end();
                 //遞歸算法，樞軸是中間元素。 
                leftshapes = std::vector<Object *>(beginning, middling);
                rightshapes = std::vector<Object *>(middling, ending);
            }
            break;
            case SplitMethod::SAH:
            {
                float nArea = centroidBounds.SurfaceArea();//算出最大的 

                int minCostCoor = 0;
                int mincostIndex = 0;
                float minCost = std::numeric_limits<float>::infinity();
                std::map<int, std::map<int, int>> indexMap;
                //indexmap用於記錄x，y，z（前一個int代表x，y，z，后一個map代表那個軸對應的map） 
                //遍歷x，y，z的划分 
                for(int i = 0; i < 3; i++)
                {
                    int bucketCount = 12;//桶的個數，這里定了12個桶，就是在某一個軸上面划分了12個區域 
                    std::vector<Bounds3> boundsBuckets;
                    std::vector<int> countBucket;
                    for(int j = 0; j < bucketCount; j++)
                    {
                        boundsBuckets.push_back(Bounds3());
                        countBucket.push_back(0); 
                    }

                    std::map<int, int> objMap;

                    for(int j = 0; j < objects.size(); j++)
                    {
                        int bid =  bucketCount * (centroidBounds.Offset(objects[j]->getBounds().Centroid()))[i];//算出對應x，y。z上的id值，這里【i】代表x，y，z 
                        if(bid > bucketCount - 1)//實質是可以划分13個區域的，將最后一個區域合並。 
                        {
                            bid = bucketCount - 1;
                        }
                        Bounds3 b = boundsBuckets[bid];
                        b = Union(b, objects[j]->getBounds().Centroid());
                        boundsBuckets[bid] = b;
                        countBucket[bid] = countBucket[bid] + 1;
                        objMap.insert(std::make_pair(j, bid)); 
                    }

                    indexMap.insert(std::make_pair(i, objMap));
                    //對於每一個划分，計算他所對應的花費，方法是對於桶中的每一個面積，計算他的花費，最后進行計算 
                    //找出這個划分。 
                    for(int j = 1; j < boundsBuckets.size(); j++)
                    {
                        Bounds3 A;
                        Bounds3 B;
                        int countA = 0;
                        int countB = 0;
                        for(int k = 0; k < j; k++)
                        {
                            A = Union(A, boundsBuckets[k]);
                            countA += countBucket[k];
                        }

                        for(int k = j; k < boundsBuckets.size(); k++)
                        {
                            B = Union(B, boundsBuckets[k]);
                            countB += countBucket[k];
                        }

                        float cost = 1 + (countA * A.SurfaceArea() + countB * B.SurfaceArea()) / nArea;//計算花費 
                        //找出這個花費。 
                        if(cost < minCost)
                        {
                            minCost = cost;
                            mincostIndex = j;
                            minCostCoor = i;
                        }
                    }
                }
                //加入左右數組，這里很重要，具體還是看那篇博客 
                for(int i = 0; i < objects.size(); i++)
                {
                    if(indexMap[minCostCoor][i] < mincostIndex)
                    {
                        leftshapes.push_back(objects[i]);
                    }
                    else
                    {
                        rightshapes.push_back(objects[i]);
                    }
                }
            }
            break;
            dafault:
            break;
        }

        assert(objects.size() == (leftshapes.size() + rightshapes.size()));
        //遞歸計算，同普通方法 
        node->left = recursiveBuild(leftshapes);
        node->right = recursiveBuild(rightshapes);

        node->bounds = Union(node->left->bounds, node->right->bounds);
    }

    return node;
}

 

 

 

最后終於完成了。

 

這一次還是很有挑戰性的，里面涉及啟發式搜索，BVH，bounds求交