GAMES101-现代计算机图形学学习笔记（作业06）

参考：https://blog.csdn.net/qq_36242312/article/details/115495482

思路

上节课的代码：

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;
//hero:这里自定义了一个eye_pos位置
    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 centerX = (float(i) + 0.5) / scene.width;
			float centerY = (float(j) + 0.5) / scene.height;

			float x = (2 * centerX - 1) * imageAspectRatio * scale;
			float y = (1 - 2 * centerY) * scale;
			
            // Ray 方向
			Vector3f dir = Vector3f(x, y, -1);
			dir = normalize(dir);

            // Ray，hero：eye到ray origin的连线
            Ray ray(eye_pos, dir);

			framebuffer[m++] = scene.castRay(ray, 0);
        }
        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);    
}

其中光线与三角形的相交的函数 Triangle::getIntersection() 采用了这次作业的框架进行改写：

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);
    double det = dotProduct(e1, pvec);
    if (fabs(det) < EPSILON)
        return inter;

    double det_inv = 1. / det;
    Vector3f tvec = ray.origin - v0;
    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;

    if (t_tmp < 0)
    {
        return inter;
    }
//interaction的distance= std::numeric_limits<double>::max();无限远
    inter.distance = t_tmp;
    inter.coords = ray(t_tmp);
    inter.happened = true;
    inter.m = m;
    inter.normal = normal;
    inter.obj = this;

    return inter;
}

 

 

 所以这里函数 IntersectP 的实现就可以写出下面的方式：

inline bool Bounds3::IntersectP(const Ray& ray, const Vector3f& invDir,
                                const std::array<int, 3>& dirIsNeg) const
{
	// hero：光线进入点，初始都定义为无限远，相当于无交点
	float tEnter = -std::numeric_limits<float>::infinity();
	// 光线离开点
	float tExit = std::numeric_limits<float>::infinity();
//hero：遍历x,y,z三个分量
	for (int i = 0; i < 3; i++)
	{
		float min = (pMin[i] - ray.origin[i]) * invDir[i];
		float max = (pMax[i] - ray.origin[i]) * invDir[i];
		// 坐标为负的话，需要进行交换
		if (!dirIsNeg[i])
		{
			std::swap(min, max);
		}
//hero：得到每个分量的交集
		tEnter = std::max(min, tEnter);
		tExit = std::min(max, tExit);
	}
    return tEnter < tExit && tExit >= 0;
}

② BVH 求交过程

BVH 本质上是把场景中的三角形用二叉树状的结构进行表示，中间结点包含这个结点所包含几何体的包围盒以及指向叶子结点的指针，叶子节点包含物体列表和包围盒，当判断光线与物体是否相交时，只需要递归的判断包围盒是否与光线相交即可：

代码描述如下：

Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const
{
    Intersection inter;

    // 光线方向
    float x = ray.direction.x;
    float y = ray.direction.y;
    float z = ray.direction.z;
	// 判断坐标是否为负，hero：false强转为int是0，true是1，这样方便操作
    std::array<int, 3> dirsIsNeg{ int(x > 0),int(y > 0),int(z > 0) };

	// 判断结点的包围盒与光线是否相交，hero：包围盒无交点终制
    if (node->bounds.IntersectP(ray, ray.direction_inv, dirsIsNeg) == false) return inter;
    //hero：到叶节点，查找交点
    if (node->left == nullptr && node->right == nullptr)
    {
        inter = node->object->getIntersection(ray);
        return inter;
    }
    
	// 递归判断子节点是否存在与光线相交的情况
	auto hit1 = getIntersection(node->left, ray);
	auto hit2 = getIntersection(node->right, ray);

    if (hit1.distance < hit2.distance)
        return hit1;
    return hit2;
}

③ SAH 的实现

SAH（Surface Area Heuristic）是一种启发式对 BVH 进行更好划分的方案，它是通过代价函数使得每一次划分都更加合理。由于 BVH 仅仅是按轴对图元进行划分，所有会有很多情况出现，不同的情况会导致计算量完全不同，如下图:

 

 

SAH 递归实现：

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

    // Compute bounds of all primitives in SVH 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;

        constexpr int nBuckets = 12;
        float minCost = std::numeric_limits<float>::max();
        int minCostSplitBucket = 0;

		// 记录按最小代价划分的索引
        int mid;
        // 记录按最小代价划分对应的维度
        int minDim = 0;
        
        int start = 0;
        int end = objects.size();

        for (int i = 0; i < end; ++i)
            centroidBounds = Union(centroidBounds, objects[i]->getBounds().Centroid());
           
        // 对三个维度分别计算最小代价
        for (int dim = 0; dim < 3; ++dim)
        {
        	// 如果图元数量小于4，直接取中间即可，减少运算时间
            if (objects.size() <= 4) {
                mid = (start + end) / 2;
                std::nth_element(&objects[start], &objects[mid], &objects[end - 1] + 1,
                    [dim](Object* a, Object* b) {
                        return a->getBounds().Centroid()[dim] < b->getBounds().Centroid()[dim];
                    });
            }
            else
            {
            	// 否则，按最小代价进行划分，首先初始化格子
                BucketInfo buckets[nBuckets];
                for (int i = 0; i < end; ++i) {
                    Vector3f offsetVec = centroidBounds.Offset(objects[i]->getBounds().Centroid());
                    int b = nBuckets * offsetVec[dim];
                    if (b == nBuckets) b = nBuckets - 1;
                    buckets[b].count++;
                    buckets[b].bounds = Union(buckets[b].bounds, objects[i]->getBounds().Centroid());
                }

				// 计算最小代价
                float cost[nBuckets - 1];
                for (int i = 0; i < nBuckets - 1; ++i) {
                    Bounds3 b0, b1;
                    int count0 = 0, count1 = 0;
                    for (int j = 0; j <= i; ++j) {
                        b0 = Union(b0, buckets[j].bounds);
                        count0 += buckets[j].count;
                    }
                    for (int j = i + 1; j < nBuckets; ++j) {
                        b1 = Union(b1, buckets[j].bounds);
                        count1 += buckets[j].count;
                    }
                    cost[i] = .125f + (count0 * b0.SurfaceArea() + count1 * b1.SurfaceArea()) / bounds.SurfaceArea();

					// 更新最小代价参数：最小代价，最小代价的格子索引，最小代价的维度
                    if (cost[i] < minCost) {
                        minCost = cost[i];
                        minCostSplitBucket = i; 
                        minDim = dim;
                    }
                }
            }
        }

		// 按最小代价进行划分
        auto pmid = std::partition(&objects[0], &objects[end - 1] + 1,
            [=](Object* pi) {
                int b = nBuckets * centroidBounds.Offset(pi->getBounds().Centroid())[minDim];
                if (b == nBuckets) b = nBuckets - 1;
                return b <= minCostSplitBucket;
            });
        mid = pmid - &objects[0];

        auto beginning = objects.begin();
        auto middling = objects.begin() + mid;
        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 = recursiveBuildSAH(leftshapes);
        node->right = recursiveBuildSAH(rightshapes);

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

    return node;
}

结果

BVH 时间：