计算机图形学笔记(四)：光线追踪与全局光照#

1
bool intersect_sphere(const Ray& ray, const Vector3f& center, float radius,
2
                     float& t_near, float& t_far) {
3
    Vector3f oc = ray.origin - center;
4
    float a = ray.direction.dot(ray.direction);
5
    float b = 2.0f * oc.dot(ray.direction);
6
    float c = oc.dot(oc) - radius * radius;
7

8
    float discriminant = b * b - 4 * a * c;
9
    if (discriminant < 0) return false;
10

11
    float sqrt_discriminant = std::sqrt(discriminant);
12
    t_near = (-b - sqrt_discriminant) / (2.0f * a);
13
    t_far = (-b + sqrt_discriminant) / (2.0f * a);
14

15
    if (t_near > t_far) std::swap(t_near, t_far);
16

17
    return t_near >= ray.t_min && t_near <= ray.t_max;
18
}

1
bool intersect_triangle(const Ray& ray, const Vector3f& v0, const Vector3f& v1, const Vector3f& v2,
2
                       float& t, float& u, float& v) {
3
    Vector3f edge1 = v1 - v0;
4
    Vector3f edge2 = v2 - v0;
5
    Vector3f h = ray.direction.cross(edge2);
6
    float a = edge1.dot(h);
7

8
    if (a > -EPSILON && a < EPSILON) return false;  // 光线平行于三角形
9

10
    float f = 1.0f / a;
11
    Vector3f s = ray.origin - v0;
12
    u = f * s.dot(h);
13

14
    if (u < 0.0f || u > 1.0f) return false;
15

16
    Vector3f q = s.cross(edge1);
17
    v = f * ray.direction.dot(q);
18

19
    if (v < 0.0f || u + v > 1.0f) return false;
20

21
    t = f * edge2.dot(q);
22

23
    return t > EPSILON && t >= ray.t_min && t <= ray.t_max;
24
}

1
bool intersect_plane(const Ray& ray, const Vector3f& point, const Vector3f& normal,
2
                    float& t) {
3
    float denom = normal.dot(ray.direction);
4
    if (std::abs(denom) < EPSILON) return false;  // 光线平行于平面
5

6
    t = normal.dot(point - ray.origin) / denom;
7
    return t >= ray.t_min && t <= ray.t_max;
8
}

1
struct OBB {
2
    Vector3f center;
3
    Vector3f axes[3];      // 三个正交轴
4
    Vector3f half_extents; // 沿各轴的半长度
5

6
    AABB to_aabb() const {
7
        Vector3f extent = Vector3f::Zero();
8
        for (int i = 0; i < 3; ++i) {
9
            extent += axes[i].cwiseAbs() * half_extents[i];
10
        }
11
        return AABB(center - extent, center + extent);
12
    }
13
};

1
class TriangleMesh {
2
private:
3
    std::vector<Vector3f> vertices;
4
    std::vector<Vector3i> triangles;
5
    AABB bounding_box;
6

7
public:
8
    bool intersect(const Ray& ray, Intersection& hit) {
9
        if (!intersect_aabb(ray, bounding_box)) return false;
10

11
        bool hit_found = false;
12
        float closest_t = ray.t_max;
13

14
        for (const auto& triangle : triangles) {
15
            Vector3f v0 = vertices[triangle[0]];
16
            Vector3f v1 = vertices[triangle[1]];
17
            Vector3f v2 = vertices[triangle[2]];
18

19
            float t, u, v;
20
            if (intersect_triangle(ray, v0, v1, v2, t, u, v) && t < closest_t) {
21
                closest_t = t;
22
                hit.t = t;
23
                hit.position = ray.at(t);
24
                hit.normal = (v1 - v0).cross(v2 - v0).normalized();
25
                hit.uv = Vector2f(u, v);
26
                hit_found = true;
27
            }
28
        }
29

30
        return hit_found;
31
    }
32
};

1
bool intersect_implicit_surface(const Ray& ray,
2
                               std::function<float(Vector3f)> f,
3
                               std::function<Vector3f(Vector3f)> gradient,
4
                               float& t) {
5
    float t_current = ray.t_min;
6
    const int max_iterations = 100;
7
    const float tolerance = 1e-6f;
8

9
    for (int i = 0; i < max_iterations; ++i) {
10
        Vector3f point = ray.at(t_current);
11
        float value = f(point);
12

13
        if (std::abs(value) < tolerance) {
14
            t = t_current;
15
            return true;
16
        }
17

18
        Vector3f grad = gradient(point);
19
        float denominator = grad.dot(ray.direction);
20

21
        if (std::abs(denominator) < tolerance) break;
22

23
        t_current -= value / denominator;
24

25
        if (t_current < ray.t_min || t_current > ray.t_max) break;
26
    }
27

28
    return false;
29
}

1
Intersection BVHAccel::getIntersection(BVHBuildNode* node, const Ray& ray) const {
2
    Intersection isect;
3
    if (!node) return isect;
4

5
    // 检查光线是否与节点包围盒相交
6
    Vector3f inv_dir = Vector3f(1.0f / ray.direction.x(),
7
                               1.0f / ray.direction.y(),
8
                               1.0f / ray.direction.z());
9
    std::array<int, 3> dir_is_neg = {
10
        ray.direction.x() > 0 ? 0 : 1,
11
        ray.direction.y() > 0 ? 0 : 1,
12
        ray.direction.z() > 0 ? 0 : 1
13
    };
14

15
    if (!node->bounds.intersect_p(ray, inv_dir, dir_is_neg)) {
16
        return isect;
17
    }
18

19
    // 叶子节点：与对象求交
20
    if (node->left == nullptr && node->right == nullptr) {
21
        return node->object->getIntersection(ray);
22
    }
23

24
    // 内部节点：递归遍历子节点
25
    Intersection hit1 = getIntersection(node->left, ray);
26
    Intersection hit2 = getIntersection(node->right, ray);
27

28
    return hit1.distance < hit2.distance ? hit1 : hit2;
29
}

1
Cost = C_trav + P_left × N_left × C_isect + P_right × N_right × C_isect

1
P_left = SA_left / SA_parent
2
P_right = SA_right / SA_parent

1
struct SAHBucket {
2
    int count = 0;
3
    AABB bounds;
4
};
5

6
float evaluate_sah(const std::vector<Object*>& objects, int split_axis, float split_pos) {
7
    AABB left_bounds, right_bounds;
8
    int left_count = 0, right_count = 0;
9

10
    for (const auto& obj : objects) {
11
        Vector3f centroid = obj->getBounds().centroid();
12
        if (centroid[split_axis] < split_pos) {
13
            left_bounds.expand(obj->getBounds());
14
            left_count++;
15
        } else {
16
            right_bounds.expand(obj->getBounds());
17
            right_count++;
18
        }
19
    }
20

21
    float cost = 1.0f + // 遍历成本
22
                (left_count * left_bounds.surface_area() +
23
                 right_count * right_bounds.surface_area()) /
24
                parent_bounds.surface_area();
25

26
    return cost;
27
}
28

29
int find_best_split_sah(const std::vector<Object*>& objects) {
30
    const int n_buckets = 12;
31
    float min_cost = INFINITY;
32
    int best_split = -1;
33

34
    for (int dim = 0; dim < 3; ++dim) {
35
        // 创建桶
36
        std::vector<SAHBucket> buckets(n_buckets);
37

38
        // 将对象分配到桶中
39
        for (const auto& obj : objects) {
40
            int bucket = n_buckets * centroid_bounds.offset(obj->getBounds().centroid())[dim];
41
            bucket = std::min(bucket, n_buckets - 1);
42
            buckets[bucket].count++;
43
            buckets[bucket].bounds.expand(obj->getBounds());
44
        }
45

46
        // 计算每个分割位置的成本
47
        for (int i = 0; i < n_buckets - 1; ++i) {
48
            AABB b0, b1;
49
            int count0 = 0, count1 = 0;
50

51
            for (int j = 0; j <= i; ++j) {
52
                b0.expand(buckets[j].bounds);
53
                count0 += buckets[j].count;
54
            }
55

56
            for (int j = i + 1; j < n_buckets; ++j) {
57
                b1.expand(buckets[j].bounds);
58
                count1 += buckets[j].count;
59
            }
60

61
            float cost = 1.0f + (count0 * b0.surface_area() + count1 * b1.surface_area()) /
62
                         bounds.surface_area();
63

64
            if (cost < min_cost) {
65
                min_cost = cost;
66
                best_split = i;
67
            }
68
        }
69
    }
70

71
    return best_split;
72
}

1
class Octree {
2
private:
3
    struct OctreeNode {
4
        AABB bounds;
5
        std::vector<Object*> objects;
6
        std::array<std::unique_ptr<OctreeNode>, 8> children;
7
        bool is_leaf;
8

9
        OctreeNode(const AABB& b) : bounds(b), is_leaf(true) {}
10
    };
11

12
    std::unique_ptr<OctreeNode> root;
13
    int max_depth;
14
    int max_objects_per_node;
15

16
public:
17
    void build(const std::vector<Object*>& objects, const AABB& scene_bounds) {
18
        root = std::make_unique<OctreeNode>(scene_bounds);
19
        build_recursive(root.get(), objects, 0);
20
    }
21

22
private:
23
    void build_recursive(OctreeNode* node, const std::vector<Object*>& objects, int depth) {
24
        if (depth >= max_depth || objects.size() <= max_objects_per_node) {
25
            node->objects = objects;
26
            return;
27
        }
28

29
        // 创建8个子节点
30
        Vector3f center = node->bounds.center();
31
        Vector3f min_pt = node->bounds.min_point;
32
        Vector3f max_pt = node->bounds.max_point;
33

34
        // 8个子立方体的包围盒
35
        std::array<AABB, 8> child_bounds = {
36
            AABB(Vector3f(min_pt.x(), min_pt.y(), min_pt.z()),
37
                 Vector3f(center.x(), center.y(), center.z())),  // 000
38
            AABB(Vector3f(center.x(), min_pt.y(), min_pt.z()),
39
                 Vector3f(max_pt.x(), center.y(), center.z())),  // 100
40
            AABB(Vector3f(min_pt.x(), center.y(), min_pt.z()),
41
                 Vector3f(center.x(), max_pt.y(), center.z())),  // 010
42
            AABB(Vector3f(center.x(), center.y(), min_pt.z()),
43
                 Vector3f(max_pt.x(), max_pt.y(), center.z())),  // 110
44
            AABB(Vector3f(min_pt.x(), min_pt.y(), center.z()),
45
                 Vector3f(center.x(), center.y(), max_pt.z())),  // 001
46
            AABB(Vector3f(center.x(), min_pt.y(), center.z()),
47
                 Vector3f(max_pt.x(), center.y(), max_pt.z())),  // 101
48
            AABB(Vector3f(min_pt.x(), center.y(), center.z()),
49
                 Vector3f(center.x(), max_pt.y(), max_pt.z())),  // 011
50
            AABB(Vector3f(center.x(), center.y(), center.z()),
51
                 Vector3f(max_pt.x(), max_pt.y(), max_pt.z()))   // 111
52
        };
53

54
        // 将对象分配到子节点
55
        std::array<std::vector<Object*>, 8> child_objects;
56
        for (const auto& obj : objects) {
57
            for (int i = 0; i < 8; ++i) {
58
                if (child_bounds[i].intersects(obj->getBounds())) {
59
                    child_objects[i].push_back(obj);
60
                }
61
            }
62
        }
63

64
        // 递归构建非空子节点
65
        node->is_leaf = false;
66
        for (int i = 0; i < 8; ++i) {
67
            if (!child_objects[i].empty()) {
68
                node->children[i] = std::make_unique<OctreeNode>(child_bounds[i]);
69
                build_recursive(node->children[i].get(), child_objects[i], depth + 1);
70
            }
71
        }
72
    }
73

74
public:
75
    Intersection intersect(const Ray& ray) const {
76
        return intersect_recursive(root.get(), ray);
77
    }
78

79
private:
80
    Intersection intersect_recursive(OctreeNode* node, const Ray& ray) const {
81
        if (!node || !node->bounds.intersect(ray)) {
82
            return Intersection();
83
        }
84

85
        if (node->is_leaf) {
86
            Intersection closest_hit;
87
            float closest_t = INFINITY;
88

89
            for (const auto& obj : node->objects) {
90
                Intersection hit = obj->getIntersection(ray);
91
                if (hit.happened && hit.distance < closest_t) {
92
                    closest_hit = hit;
93
                    closest_t = hit.distance;
94
                }
95
            }
96
            return closest_hit;
97
        }
98

99
        // 遍历子节点
100
        Intersection closest_hit;
101
        float closest_t = INFINITY;
102

103
        for (const auto& child : node->children) {
104
            if (child) {
105
                Intersection hit = intersect_recursive(child.get(), ray);
106
                if (hit.happened && hit.distance < closest_t) {
107
                    closest_hit = hit;
108
                    closest_t = hit.distance;
109
                }
110
            }
111
        }
112

113
        return closest_hit;
114
    }
115
};