这里回顾GAMES101 HW6，这次作业的内容是加速结构。

课程主页：

课程作业：

课程视频：

参考资料：

Render

微调之前的代码即可：

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
        Vector3f dir = Vector3f(x, y, -1); // Don't forget to normalize this direction!
        dir = normalize(dir);
        // change
        Ray ray(eye_pos, dir);
        framebuffer[m++] = scene.castRay(ray, 0);
        // Don't forget to normalize this direction!

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

getIntersection

也是沿用之前的代码，注意这里需要传递物体的信息：

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

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

    double det_inv = 1. / det;
    // 对应s
    Vector3f tvec = ray.origin - v0;
    u = dotProduct(tvec, pvec) * det_inv;
    if (u < 0 || u > 1)
        return inter;
    // 对应s2
    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.happened = true;
    // 交点坐标
    inter.coords = ray.origin + t_tmp * ray.direction;
    // 交点所在平面法线
    inter.normal = normal;
    // 距离
    inter.distance = t_tmp;
    // 交点对应的物体(三角形)
    inter.obj = this;
    // 材质
    inter.m = m;

    return inter;
}

IntersectP

这里首先回顾判断物体的BVH求交的过程，关键思想为：

  • 光线只有在进入所有成对的平板时才会进入盒子；
  • 只要光线离开任何一对平板，光线就会离开盒子；
  • 所以对于每一对平板，计算tmin和tmax（负数也ok）
  • 对于3D box，$t_{\mathrm{enter}} = \max{t_{\min}} $，$t_{\mathrm{exit}} = \min{t_{\max}}$
  • 如果$t_{\mathrm{enter}} < t_{\mathrm{exit}}$，我们知道光线在盒子里停留了一段时间（所以必须相交！）；
  • 然而，射线不是直线；
    • 应该检查$t$对应的实际情形是否正确；
  • 如果$t_{\mathrm{exit}}$怎么办？
    • 盒子在光线“后面”——没有交点！
  • 如果$t_{\mathrm{exit}}\ge 0$并且$t_{\mathrm{enter}}<0$怎么办？
    • 光线的起点在盒子内——有交点！
  • 综上所述，ray和AABB当且仅当：
    • $t_{\mathrm{enter}} < t_{\mathrm{exit}}$并且$t_{\mathrm{exit}}\ge 0$；

所以首先是计算每个方向的进入和离开时间，利用轴对齐情形的特殊计算公式：

注意该公式默认了法线的各个分量为正，对于负数情形，需要交换离开和进入的时间：

inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,
                                const std::array<int, 3>& dirIsNeg) const
{
    // 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
    Vector3f origin = ray.origin;
    Vector3f t_enter = (pMin - origin) * invDir;
    Vector3f t_exit = (pMax - origin) * invDir;
    // 根据方向调整enter和exit
    if (dirIsNeg[0]) {
        float d = t_enter.x;
        t_enter.x = t_exit.x;
        t_exit.x = d;
    }
    if (dirIsNeg[1]) {
        float d = t_enter.y;
        t_enter.y = t_exit.y;
        t_exit.y = d;
    }
    if (dirIsNeg[2]) {
        float d = t_enter.z;
        t_enter.z = t_exit.z;
        t_exit.z = d;
    }
    float t_min = std::max(t_enter.x, std::max(t_enter.y, t_enter.z));
    float t_max = std::min(t_exit.x, std::min(t_exit.x, t_exit.z));

    return (t_min < t_max && t_max >= 0);
}

getIntersection

实现如下算法即可：

Intersect(Ray ray, BVH node) {
    if (ray misses node.bbox) 
    	return;
    if (node is a leaf node)
		test intersection with all objs;
    	return closest intersection;
    hit1 = Intersect(ray, node.child1);
    hit2 = Intersect(ray, node.child2);
    return the closer of hit1, hit2;
}

代码如下：

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
    // TODO Traverse the BVH to find intersection
    Intersection isect;
    Vector3f invDir = ray.direction_inv;
    std::array<int, 3> dirIsNeg {invDir.x < 0, invDir.y < 0, invDir.z < 0};
    // 相交
    if (node->bounds.IntersectP(ray, invDir, dirIsNeg)) {
        if (node == NULL) {
            return isect;
        }
        if (node->left == NULL && node->right == NULL) {
            isect = node->object->getIntersection(ray);
            return isect;
        }
        Intersection isect_left = getIntersection(node->left, ray);
        Intersection isect_right = getIntersection(node->right, ray);
        if (isect_left.distance < isect_right.distance) {
            isect = isect_left;
        } else {
            isect = isect_right;
        }
    }

    return isect;
}

SVH

SVH可以先参考博客)，具体代码如下：

BVHBuildNode* BVHAccel::recursiveBuildSVH(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 = recursiveBuildSVH(std::vector{objects[0]});
        node->right = recursiveBuildSVH(std::vector{objects[1]});

        node->bounds = Union(node->left->bounds, node->right->bounds);
        return node;
    }
    else {
        // 原则
        // 找到最宽松的轴进行划分
        Bounds3 centroidBounds;
        Bounds3 maxBounds;
        for (int i = 0; i < objects.size(); ++i) {
            centroidBounds = Union(centroidBounds, objects[i]->getBounds().Centroid());
            maxBounds = Union(maxBounds, objects[i]->getBounds());
        }
        float total_size = maxBounds.SurfaceArea();

        int dim = centroidBounds.maxExtent();
        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;
        }

        // 策略
        float min_cos = INFINITY;
        // 计算面积
        int n = objects.size();
        // 从左往右前i个
        std::vector<float> left_size(n + 1, 0);
        // 从右往左前i个
        std::vector<float> right_size(n + 1, 0);
        std::vector<float> score(n + 1, 0);
        Bounds3 left_bound;
        Bounds3 right_bound;
        for (int i = 1; i <= n; i++) {
            left_bound = Union(left_bound, objects[i - 1]->getBounds());
            left_size[i] = left_bound.SurfaceArea() / total_size;
        }
        for (int i = 1; i <= n; i++) {
            right_bound = Union(right_bound, objects[n - i]->getBounds());
            right_size[i] = right_bound.SurfaceArea() / total_size;
        }
        // 计算总分
        for (int i = 0; i <= n; i++) {
            score[i] = i * left_size[i] + (n - i) * right_size[n - i];
        }
        // 找到分数最小的位置
        float s = INFINITY;
        int index = -1;
        for (int i = 0; i <= n; i++) {
            if (score[i] < s) {
                s = score[i];
                index = i;
            }
        }
        // 保证划分后都有元素
        if (index == n) {
            index--;
        } else if (index == 0) {
            index++;
        }

        auto beginning = objects.begin();
        auto middling = objects.begin() + index;
        auto ending = objects.end();

        auto leftshapes = std::vector<Object*>(beginning, middling);
        auto rightshapes = std::vector<Object*>(middling, ending);

        assert(objects.size() == (leftshapes.size() + rightshapes.size()));

        node->left = recursiveBuildSVH(leftshapes);
        node->right = recursiveBuildSVH(rightshapes);

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

    return node;
}

示例

几个注意点

  • 第一个注意点是源代码注释中dirIsNeg=[int(x>0),int(y>0),int(z>0)]定义错误；
  • 第二个注意点是如果在getIntersection中讨论hit1和hit2是否为空，则会出现问题；