CMU 15-462 A3 实现笔记

上次写了 A1 的实现笔记,A2 感觉基本都是有很多细节的大模拟所以不想写实现笔记了(代码在 GitHub 上),一开始写得感觉都还可以,最后写完 downsample 和 resample 时发现之前写的有一个 corner case 没有想到,导致 RE,搞得心态也有点崩。。。(不过也算是完成 A2 了吧,尽管不是完美完成。。。)

A3 是路径追踪的内容,算是挺喜欢的内容了。。。实现部分有参考 Ubpa 的实现:https://github.com/Ubpa/CMU_15_462

之前也跟过 Ray Tracing in One Weekend 三本,但感觉还是有不太理解的地方。想着之后可以把 PBR 的那本看一看好了。

笔记和代码有放在 GitHub 上:https://github.com/PepcyCh/cmu15462-notes

课程: Computer Graphics (CMU 15-462)

Task 1: Generating Camera Rays

要求实现摄像机射出光线的代码,在 camera.cpp 中:

1
2
3
4
5
6
7
8
9
10
11
Ray Camera::generate_ray(double x, double y) const {
// compute position of the input sensor sample coordinate on the
// canonical sensor plane one unit away from the pinhole.

Vector3D tar = Vector3D((x - .5) * screenW, (y - .5) * screenH, -screenDist);
Vector3D dir = c2w * tar;
dir.normalize();
Vector3D org = pos;

return Ray(org, dir);
}

- .5 是把原点从左下角调至中心。

一开始没有单位化,一切都很正常,直到我实现最后一项时出了问题。。。(当时想当然地以为方向向量是单位化的,却看到自己没有这么做。。。)

然后把各种 Sampler 实现了一波:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
Vector2D JitteredSampler::get_sample() const {
static std::vector<Vector2D> samples(ns_aa);
static int curr = ns_aa;
if (curr == ns_aa) {
int w = std::sqrt(ns_aa);
curr = 0;
for (int i = 0; i < w; i++) {
for (int j = 0; j < w; j++) {
double x = (i + double(std::rand()) / RAND_MAX) / w;
double y = (j + double(std::rand()) / RAND_MAX) / w;
samples[curr++] = Vector2D(x, y);
}
}
curr = 0;
}
return samples[curr++];
}

Vector2D MultiJitteredSampler::get_sample() const {
static std::vector<Vector2D> samples(ns_aa);
static int curr = ns_aa;
if (curr == ns_aa) {
int w = std::sqrt(ns_aa);
curr = 0;
int u = 0, v = 0;
for (int i = 0; i < w; i++) {
for (int j = 0; j < w; j++) {
double x = (u + double(std::rand()) / RAND_MAX) / w;
double y = (v + double(std::rand()) / RAND_MAX) / w;
samples[curr++] = Vector2D(x, y);
if (v == w - 1) {
++u;
v = 0;
} else {
++v;
}
}
}
std::random_shuffle(samples.begin(), samples.begin() + curr);
curr = 0;
for (int i = 0; i < w; i++) {
for (int j = 0; j < w; j++) {
double x = (i + samples[curr].x) / w;
double y = (j + samples[curr].y) / w;
samples[curr++] = Vector2D(x, y);
}
}
curr = 0;
}
return samples[curr++];
}

Vector2D NRooksSampler::get_sample() const {
static std::vector<Vector2D> samples(ns_aa);
static int curr = ns_aa;
if (curr == ns_aa) {
curr = 0;
std::vector<double> X(ns_aa), Y(ns_aa);
for (int i = 0; i < ns_aa; i++) {
X[i] = (i + double(std::rand()) / RAND_MAX) / ns_aa;
Y[i] = (i + double(std::rand()) / RAND_MAX) / ns_aa;
}
std::random_shuffle(Y.begin(), Y.end());
for (int i = 0; i < ns_aa; i++) {
samples[i] = Vector2D(X[i], Y[i]);
}
curr = 0;
}
return samples[curr++];
}

static double radicalInverse3(int i) {
double b = 1. / 3.;
double t = b;
double res = 0;
while (i) {
res += t * (i % 3);
t *= b;
i /= 3;
}
return res;
}
static double radicalInverse2(int i) {
double t = 1. / 2.;
double res = 0;
while (i) {
res += t * (i & 1);
t /= 2.;
i >>= 1;
}
return res;
}

Vector2D SobolSampler::get_sample() const {
static std::vector<Vector2D> samples(ns_aa);
static bool init = false;
static int curr = ns_aa;
if (!init) {
std::vector<unsigned> C(ns_aa), V(32);
for (int i = 0; i < ns_aa; i++) {
int w = i;
C[i] = 1;
while (w & 1) {
++C[i];
w >>= 1;
}
}
for (int i = 1; i <= 31; i++) V[i] = 1 << (32 - i);

std::vector<unsigned> X(ns_aa);
X[0] = 0;
for (int i = 1; i < ns_aa; i++) {
X[i] = X[i - 1] ^ V[C[i - 1]];
samples[i].x = X[i] / std::pow(2., 32.);
}

V[1] = 1 << 31;
for (int i = 2; i <= 31; i++) V[i] = V[i - 1] ^ (V[i - 1] >> 1);
X[0] = 0;
for (int i = 1; i < ns_aa; i++) {
X[i] = X[i - 1] ^ V[C[i - 1]];
samples[i].y = X[i] / std::pow(2., 32.);
}
init = true;
}
if (curr == ns_aa) curr = 0;
return samples[curr++];
}

Vector2D HaltonSampler::get_sample() const {
static std::vector<Vector2D> samples(ns_aa);
static bool init = false;
static int curr = ns_aa;
if (!init) {
for (int i = 0; i < ns_aa; i++) {
double x = radicalInverse2(i);
double y = radicalInverse3(i);
samples[i] = Vector2D(x, y);
}
init = true;
}
if (curr == ns_aa) curr = 0;
return samples[curr++];
}

Vector2D HammersleySampler::get_sample() const {
static std::vector<Vector2D> samples(ns_aa);
static bool init = false;
static int curr = ns_aa;
if (!init) {
for (int i = 0; i < ns_aa; i++) {
double x = double(i) / ns_aa;
double y = radicalInverse2(i);
samples[i] = Vector2D(x, y);
}
init = true;
}
if (curr == ns_aa) curr = 0;
return samples[curr++];
}

这个实现中,除 Sobol、Halton 和 Hammersley 都不能在多线程下正常工作。。。

之后我一直用的是 Sobol,在 pathtracer.cpp 中的构造函数部分做了修改:

1
2
3
4
5
6
7
// gridSampler = new UniformGridSampler2D();
// gridSampler = new JitteredSampler(ns_aa);
// gridSampler = new MultiJitteredSampler(ns_aa);
// gridSampler = new NRooksSampler(ns_aa);
gridSampler = new SobolSampler(ns_aa);
// gridSampler = new HaltonSampler(ns_aa);
// gridSampler = new HammersleySampler(ns_aa);

Task 2: Intersecting Triangles and Spheres

球体,在 static_sence/sphere.cpp:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
bool Sphere::test(const Ray& r, double& t1, double& t2) const {
// Implement ray - sphere intersection test.
// Return true if there are intersections and writing the
// smaller of the two intersection times in t1 and the larger in t2.

Vector3D oc = r.o - o;
double a = r.d.norm2();
double b = dot(oc, r.d);
double c = oc.norm2() - r2;
double disc = b * b - a * c;
if (disc > 0) {
double sqrtd = sqrt(disc);
t1 = (-b - sqrtd) / a;
t2 = (-b + sqrtd) / a;
return true;
}
return false;
}

bool Sphere::intersect(const Ray& r) const {
// Implement ray - sphere intersection.
// Note that you might want to use the the Sphere::test helper here.
double t1, t2;
if (test(r, t1, t2)) {
return !(t1 > r.max_t || t2 < r.min_t);
}

return false;
}

bool Sphere::intersect(const Ray& r, Intersection* isect) const {
// Implement ray - sphere intersection.
// Note again that you might want to use the the Sphere::test helper here.
// When an intersection takes place, the Intersection data should be updated
// correspondingly.

double t1, t2;
if (test(r, t1, t2)) {
if (t1 > r.max_t || t2 < r.min_t) return false;
isect->primitive = this;
isect->bsdf = get_bsdf();
if (t1 > r.min_t) {
r.max_t = isect->t = t1;
isect->n = (r.o + t1 * r.d - o);
isect->n.normalize();
} else {
r.max_t = isect->t = t2;
isect->n = (r.o + t2 * r.d - o);
isect->n.normalize();
}
return true;
}

return false;
}

三角形,在 static_sence/triangle.cpp:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
bool Triangle::intersect(const Ray& r) const {
Vector3D p1 = mesh->positions[v1];
Vector3D p2 = mesh->positions[v2];
Vector3D p3 = mesh->positions[v3];

Vector3D e1 = p2 - p1;
Vector3D e2 = p3 - p1;
Vector3D s = r.o - p1;

double det = dot(cross(e1, r.d), e2);
if (det != 0) {
double du = -dot(cross(s, e2), r.d);
double dv = dot(cross(e1, r.d), s);
double dt = -dot(cross(s, e2), e1);
double u = du / det;
double v = dv / det;
double t = dt / det;
if (u < 0 || v < 0 || 1 - u - v < 0) return false;
return !(t < r.min_t || t > r.max_t);
} else {
double t;
if (doesRayIntersectSegment(r, p1, p2, t)) return !(t < r.min_t || t > r.max_t);
if (doesRayIntersectSegment(r, p1, p3, t)) return !(t < r.min_t || t > r.max_t);
if (doesRayIntersectSegment(r, p2, p3, t)) return !(t < r.min_t || t > r.max_t);
}

return false;
}

bool Triangle::intersect(const Ray& r, Intersection* isect) const {
// implement ray-triangle intersection. When an intersection takes
// place, the Intersection data should be updated accordingly

Vector3D p1 = mesh->positions[v1];
Vector3D p2 = mesh->positions[v2];
Vector3D p3 = mesh->positions[v3];
Vector3D n1 = mesh->normals[v1];
Vector3D n2 = mesh->normals[v2];
Vector3D n3 = mesh->normals[v3];

Vector3D e1 = p2 - p1;
Vector3D e2 = p3 - p1;
Vector3D s = r.o - p1;

double det = dot(cross(e1, r.d), e2);
if (det != 0) {
double du = -dot(cross(s, e2), r.d);
double dv = dot(cross(e1, r.d), s);
double dt = -dot(cross(s, e2), e1);
double u = du / det;
double v = dv / det;
double t = dt / det;
if (u < 0 || v < 0 || 1 - u - v < 0) return false;
if (t < r.min_t || t > r.max_t) return false;
isect->n = u * n1 + v * n2 + (1 - u - v) * n3;
if (dot(isect->n, r.d) > 0) isect->n = -isect->n;
r.max_t = isect->t = t;
isect->primitive = this;
isect->bsdf = get_bsdf();
return true;
} else {
Matrix3x3 inv = Matrix3x3();
inv.column(0) = p1;
inv.column(1) = p2;
inv.column(2) = p3;
inv = inv.inv();
double t, u, v;
bool inter = false;
if (doesRayIntersectSegment(r, v1, v2, t)) {
if (t < r.min_t || t > r.max_t) return false;
Vector3D ret = inv * (r.o + t * r.d);
u = ret.x;
v = ret.y;
inter = true;
}
if (doesRayIntersectSegment(r, v1, v3, t)) {
if (t < r.min_t || t > r.max_t) return false;
Vector3D ret = inv * (r.o + t * r.d);
u = ret.x;
v = ret.y;
inter = true;
}
if (doesRayIntersectSegment(r, v2, v3, t)) {
if (t < r.min_t || t > r.max_t) return false;
Vector3D ret = inv * (r.o + t * r.d);
u = ret.x;
v = ret.y;
inter = true;
}
if (!inter) return false;
isect->n = u * n1 + v * n2 + (1 - u - v) * n3;
if (dot(isect->n, r.d) > 0) isect->n = -isect->n;
r.max_t = isect->t = t;
isect->primitive = this;
isect->bsdf = get_bsdf();
return true;
}
}

感觉当光线与三角面平行时的处理写得有点丑。

Task 3: Implementing a Bounding Volume Hierarchy (BVH)

先实现了包围盒的光线检测,在 bbox.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
bool BBox::intersect(const Ray &r, double &t0, double &t1) const {
// Implement ray - bounding box intersection test
// If the ray intersected the bounding box within the range given by
// t0, t1, update t0 and t1 with the new intersection times.

t0 = r.min_t;
t1 = r.max_t;
for (int d = 0; d < 3; d++) {
double tt0 = std::min((min[d] - r.o[d]) * r.inv_d[d], (max[d] - r.o[d]) * r.inv_d[d]);
double tt1 = std::max((min[d] - r.o[d]) * r.inv_d[d], (max[d] - r.o[d]) * r.inv_d[d]);
t0 = std::max(t0, tt0);
t1 = std::min(t1, tt1);
if (t0 > t1) return false;
}

return t0 < r.max_t;
}

写的时候去看了当时跟 Ray Tracing in One Weeken 的代码,发现当时包围盒检查有两处 min 和 max 写的是反的。。。难怪当时跑得那么慢。。。

建 BVH 树和查询都写成非递归的了,在 bvh.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
BVHAccel::BVHAccel(const std::vector<Primitive *> &_primitives,
size_t max_leaf_size) {
static const int B = 16;

this->primitives = _primitives;

// Construct a BVH from the given vector of primitives and maximum leaf
// size configuration. The starter code build a BVH aggregate with a
// single leaf node (which is also the root) that encloses all the
// primitives.

BBox bb;
for (int i = 0; i < primitives.size(); i++) {
bb.expand(primitives[i]->get_bbox());
}
root = new BVHNode(bb, 0, primitives.size());

std::queue<BVHNode *> q;
q.push(root);
while (!q.empty()) {
auto u = q.front();
q.pop();
if (u->range <= max_leaf_size) continue;

int bestD = -1, bestI = -1;
double SN = u->bb.surface_area(), bestC = DBL_MAX;

for (int d = 0; d < 3; d++) {
std::vector<BBox> boxes(B);
std::vector<std::vector<Primitive*>> prims(B);
double min = u->bb.min[d], max = u->bb.max[d];
double length = (max - min) / B;
if (length == 0) continue;

for (size_t i = 0; i < u->range; i++) {
BBox cb = primitives[u->start + i]->get_bbox();
double p = cb.centroid()[d];
int buc = clamp<int>((p - min) / length, 0, B - 1);
prims[buc].push_back(primitives[u->start + i]);
boxes[buc].expand(cb);
}

for (int i = 1; i < B; i++) {
BBox lb, rb;
int ln = 0, rn = 0;
for (int j = 0; j < i; j++) {
lb.expand(boxes[j]);
ln += prims[j].size();
}
for (int j = i; j < B; j++) {
rb.expand(boxes[j]);
rn += prims[j].size();
}
double SA = lb.surface_area(), SB = rb.surface_area();
double C = SA / SN * ln + SB / SN * rn;
if (C < bestC) {
bestD = d;
bestI = i;
bestC = C;
}
}
}

double min = u->bb.min[bestD], max = u->bb.max[bestD];
double length = (max - min) / B;
BBox lb, rb;
std::vector<Primitive *> lp, rp;
for (size_t i = 0; i < u->range; i++) {
BBox cb = primitives[u->start + i]->get_bbox();
double p = cb.centroid()[bestD];
int buc = clamp<int>((p - min) / length, 0, B - 1);
if (buc < bestI) {
lb.expand(cb);
lp.push_back(primitives[u->start + i]);
} else {
rb.expand(cb);
rp.push_back(primitives[u->start + i]);
}
}

if (lp.size() == 0 || lp.size() == u->range) {
lb = BBox(), rb = BBox();
int hn = u->range / 2;
for (int i = 0; i < hn; i++) lb.expand(primitives[u->start + i]->get_bbox());
for (int i = hn; i < u->range; i++) rb.expand(primitives[u->start + i]->get_bbox());
u->l = new BVHNode(lb, u->start, hn);
u->r = new BVHNode(rb, u->start + hn, u->range - hn);
q.push(u->l);
q.push(u->r);
} else {
int p = 0;
for (auto prim : lp) {
primitives[u->start + p] = prim;
++p;
}
int ln = p;
for (auto prim : rp) {
primitives[u->start + p] = prim;
++p;
}
u->l = new BVHNode(lb, u->start, ln);
u->r = new BVHNode(rb, u->start + ln, u->range - ln);
q.push(u->l);
q.push(u->r);
}
}
}


BVHAccel::~BVHAccel() {
// Implement a proper destructor for your BVH accelerator aggregate

std::queue<BVHNode *> q;
q.push(root);
while (!q.empty()) {
auto u = q.front();
q.pop();
if (u->l) q.push(u->l);
if (u->r) q.push(u->r);
delete u;
}
}

// ...

bool BVHAccel::intersect(const Ray &ray) const {
// Implement ray - bvh aggregate intersection test. A ray intersects
// with a BVH aggregate if and only if it intersects a primitive in
// the BVH that is not an aggregate.

bool hit = false;

std::stack<BVHNode *> s;
s.push(root);
while (!s.empty()) {
auto u = s.top();
s.pop();

double t0, t1;
if (u->bb.intersect(ray, t0, t1)) {
if (u->isLeaf()) {
for (size_t i = 0; i < u->range; i++) {
if (primitives[u->start + i]->intersect(ray)) {
hit = true;
break;
}
}
if (hit) break;
} else {
if (u->l) s.push(u->l);
if (u->r) s.push(u->r);
}
}
}

return hit;
}

bool BVHAccel::intersect(const Ray &ray, Intersection *isect) const {
// Implement ray - bvh aggregate intersection test. A ray intersects
// with a BVH aggregate if and only if it intersects a primitive in
// the BVH that is not an aggregate. When an intersection does happen.
// You should store the non-aggregate primitive in the intersection data
// and not the BVH aggregate itself.

bool hit = false;
isect->t = ray.max_t;

std::stack<BVHNode *> s;
s.push(root);
while (!s.empty()) {
auto u = s.top();
s.pop();

double t0, t1;
if (u->bb.intersect(ray, t0, t1)) {
if (t0 <= isect->t) {
if (u->isLeaf()) {
for (size_t i = 0; i < u->range; i++) {
if (primitives[u->start + i]->intersect(ray, isect)) {
hit = true;
ray.max_t = isect->t;
}
}
} else {
double lt0, lt1;
bool lhit = false;
if (u->l) lhit = u->l->bb.intersect(ray, lt0, lt1);
double rt0, rt1;
bool rhit = false;
if (u->r) rhit = u->r->bb.intersect(ray, rt0, rt1);

if (lhit && rhit) {
if (lt0 < rt0) {
s.push(u->r);
s.push(u->l);
} else {
s.push(u->l);
s.push(u->r);
}
} else if (lhit) {
s.push(u->l);
} else if (rhit) {
s.push(u->r);
}
}
}
}
}

return hit;
}

Task 4: Implementing Shadow Rays

pathtracer.cppPathTracer::trace_ray 方法中的相应部分增加光线的检测就好了:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
if (!isect.bsdf->is_delta()) {
Vector3D dir_to_light;
float dist_to_light;
float pr;

// ### Estimate direct lighting integral
for (SceneLight* light : scene->lights) {

// no need to take multiple samples from a point/directional source
int num_light_samples = light->is_delta_light() ? 1 : ns_area_light;

// integrate light over the hemisphere about the normal
for (int i = 0; i < num_light_samples; i++) {

// returns a vector 'dir_to_light' that is a direction from
// point hit_p to the point on the light source. It also returns
// the distance from point x to this point on the light source.
// (pr is the probability of randomly selecting the random
// sample point on the light source -- more on this in part 2)
const Spectrum& light_L = light->sample_L(hit_p, &dir_to_light, &dist_to_light, &pr);

// convert direction into coordinate space of the surface, where
// the surface normal is [0 0 1]
const Vector3D& w_in = w2o * dir_to_light;
if (w_in.z < 0) continue;

// note that computing dot(n,w_in) is simple
// in surface coordinates since the normal is (0,0,1)
double cos_theta = w_in.z;

// evaluate surface bsdf
const Spectrum& f = isect.bsdf->f(w_out, w_in);

// (Task 4) Construct a shadow ray and compute whether the intersected surface is
// in shadow. Only accumulate light if not in shadow.
Ray shadow(hit_p + hit_n * EPS_D, dir_to_light);
if (!bvh->intersect(shadow)) {
L_out += (cos_theta / (num_light_samples * pr)) * f * light_L;
}
}
}
}

其实这一段代码是做了重要性采样。

Task 5: Adding Path Tracing

实现 sampler.cpp 中的 cosine weighted sampler:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Vector3D CosineWeightedHemisphereSampler3D::get_sample(float *pdf) const {
// You may implement this, but don't have to.
double Xi1 = double(std::rand()) / RAND_MAX;
double Xi2 = double(std::rand()) / RAND_MAX;

double sintheta = sqrt(Xi1);
double costheta = sqrt(1 - Xi1);
double phi = 2.0 * PI * Xi2;

double xs = sintheta * cos(phi);
double ys = sintheta * sin(phi);
double zs = costheta;

*pdf = costheta / PI;

return Vector3D(xs, ys, zs);
}

然后在 raytracer.cpp 中的 PathTracer::trace_ray 方法中递归:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// (1) randomly select a new ray direction (it may be
// reflection or transmittence ray depending on
// surface type -- see BSDF::sample_f()
Vector3D w_in;
float pdf;
Spectrum fr = isect.bsdf->sample_f(w_out, &w_in, &pdf);

// (2) potentially terminate path (using Russian roulette)
double prob = 1.0;
if (fr.illum() < 0.5)
prob = 0.5;
if (double(rand()) / RAND_MAX > prob)
return L_out;

// (3) evaluate weighted reflectance contribution due
// to light from this direction
// Spectrum fr = isect.bsdf->f(w_out, w_in);
Ray ri(hit_p, o2w * w_in, int(r.depth + 1));
ri.min_t = EPS_D;
Spectrum Li = trace_ray(ri);
L_out += fr * Li * (std::abs(w_in.z) / (pdf * prob));

这里是修改 ri.min_t 而不是改变起点是为了处理要进入球体的内部的光线。

Task 6: Adding New Materials

全部在 bsdf.cpp 中。

diffuse 应该是前一个 task 个的?

1
2
3
4
5
6
7
8
9
10
11
12
// Diffuse BSDF //

Spectrum DiffuseBSDF::f(const Vector3D& wo, const Vector3D& wi) {
return albedo * (1.0 / PI);
}

Spectrum DiffuseBSDF::sample_f(const Vector3D& wo, Vector3D* wi, float* pdf) {
// Implement DiffuseBSDF
*wi = sampler.get_sample(pdf);
if (wo.z < 0) wi->z *= -1.;
return albedo * (1.0 / PI);
}

然后是折射与反射的函数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void BSDF::reflect(const Vector3D& wo, Vector3D* wi) {
// Implement reflection of wo about normal (0,0,1) and store result in wi.
wi->x = -wo.x;
wi->y = -wo.y;
wi->z = wo.z;
}

bool BSDF::refract(const Vector3D& wo, Vector3D* wi, float ior) {
// Use Snell's Law to refract wo surface and store result ray in wi.
// Return false if refraction does not occur due to total internal reflection
// and true otherwise. When dot(wo,n) is positive, then wo corresponds to a
// ray entering the surface through vacuum.

double k = wo.z >= 0 ? 1.0 / ior : ior;
double d = 1 - (1 - wo.z * wo.z) * k * k;
if (d < 0) return false;
wi->x = -wo.x * k;
wi->y = -wo.y * k;
wi->z = sqrt(d);
if (wo.z >= 0) wi->z *= -1.;

return true;
}

然后就是我有点懵的镜面与玻璃的 brdf 了。。。看了 Ubpa 的代码,发现结果确实和参考很相像,但我有一点不是很清楚为什么是这样的。。。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
// Mirror BSDF //

Spectrum MirrorBSDF::f(const Vector3D& wo, const Vector3D& wi) {
return 1.0 / fabs(wi.z) * reflectance;
}

Spectrum MirrorBSDF::sample_f(const Vector3D& wo, Vector3D* wi, float* pdf) {
// Implement MirrorBSDF
reflect(wo, wi);
*pdf = 1.0;
return 1.0 / fabs(wi->z) * reflectance;
}

// Glass BSDF //
double schlick(double cosine, double ior) {
double r0 = (1 - ior) / (1 + ior);
r0 = r0 * r0;
return r0 + (1 - r0) * pow(1 - cosine, 5);
}

Spectrum GlassBSDF::f(const Vector3D& wo, const Vector3D& wi) {
Vector3D temp;
bool canRefract = refract(wo, &temp, ior);
double fresnel = canRefract ? schlick(fabs(wo.z), ior) : 1.0;
double k = wo.z >= 0 ? 1.0 / ior : ior;
if (wo.z * wi.z >= 0) return fresnel / fabs(wi.z) * reflectance;
else return k * k * (1 - fresnel) / fabs(wi.z) * transmittance;
}

Spectrum GlassBSDF::sample_f(const Vector3D& wo, Vector3D* wi, float* pdf) {
// Compute Fresnel coefficient and either reflect or refract based on it.
Vector3D refle, refra;
bool canRefract = refract(wo, &refra, ior);
reflect(wo, &refle);
double fresnel = canRefract ? schlick(fabs(wo.z), ior) : 1.0;
double rnd = double(rand()) / RAND_MAX;

Spectrum retf;
if (rnd <= fresnel) {
*wi = refle;
retf = fresnel / fabs(wi->z) * reflectance;
*pdf = fresnel;
} else {
double k = wo.z >= 0 ? 1.0 / ior : ior;
*wi = refra;
retf = k * k * (1 - fresnel) / fabs(wi->z) * transmittance;
*pdf = 1.0 - fresnel;
}

return retf;
}

另外,实现的时候没有管 roughness 。。。应该是指把反射或折射光线在一定范围内扰动以形成类似金属表面样式的结果。

Task 7: Infinite Environment Lighting

和天空盒很像不过是球形的,在 static_sence/environment_light{.h, .cpp}

首先实现 AliasTable:

1
2
3
4
5
6
7
8
9
10
11
12
class AliasTable {
private:
struct Item {
int id0, id1;
double ratio;
};
std::vector<Item> items;

public:
void init(const std::vector<double>& vec);
int sample(double p) const;
};
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
void EnvironmentLight::AliasTable::init(const std::vector<double>& vec) {
int N = vec.size();
items.resize(N);
for (int i = 0; i < N; i++) {
items[i].id0 = i;
items[i].id1 = -1;
items[i].ratio = vec[i];
}

double mid = 1.0 / N;
int rich = -1, poor = -1;
for (int i = 0; i < N; i++) if (items[i].ratio < mid) {
poor = i;
break;
}
for (int i = 0; i < N; i++) if (items[i].ratio > mid) {
rich = i;
break;
}

int poor_max = poor;
while (rich != -1 && poor != -1) {
double diff = mid - items[poor].ratio;
items[poor].id1 = rich;
items[poor].ratio = mid;
items[rich].ratio -= diff;

int temp_poor = -1;
if (items[rich].ratio < mid && rich < poor) {
temp_poor = rich;
} else {
for (int i = poor_max; i < N; i++) if (items[i].ratio < mid) {
temp_poor = i;
break;
}
}
poor = temp_poor;
poor_max = std::max(poor_max, poor);

int temp_rich = -1;
for (int i = rich; i < N; i++) if (items[i].ratio > mid) {
temp_rich = i;
break;
}
rich = temp_rich;
}
}

int EnvironmentLight::AliasTable::sample(double p) const {
int id = p;
double left = p - id;
return left <= items[id].ratio ? items[id].id0 : items[id].id1;
}

EnvironmentLight::EnvironmentLight(const HDRImageBuffer* envMap)
: envMap(envMap) {
int w = envMap->w, h = envMap->h;

probs.resize(w * h);
double sum = 0;
for (int j = 0; j < h; j++) {
double theta = (j + 0.5) / h * PI;
double sintheta = sin(theta);
for (int i = 0; i < w; i++) {
int id = i + j * w;
probs[id] = envMap->data[id].illum() * sintheta;
sum += probs[id];
}
}
for (double& p : probs) p /= sum;
table.init(probs);
}

之后用 AliasTable 采样实现 sample_L:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
Spectrum EnvironmentLight::sample_L(const Vector3D& p, Vector3D* wi,
float* distToLight, float* pdf) const {
*distToLight = INF_F;

double rnd = double(rand()) / RAND_MAX * probs.size();
int id = table.sample(rnd);
*pdf = probs[id];

int x = id % envMap->w;
int y = id / envMap->w;
double theta = PI * (y + double(rand()) / RAND_MAX) / envMap->h;
double phi = 2 * PI * (x + double(rand() / RAND_MAX)) / envMap->w;

wi->x = sin(theta) * cos(phi);
wi->z = sin(theta) * sin(phi);
wi->y = cos(theta);

return sample_dir(*wi);
}

用双线性插值实现 sample_dir,因为完全没用上 Ray::o,所以直接把参数改成 Vector3D 了。。。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
Spectrum EnvironmentLight::sample_dir(const Vector3D& r) const {
int w = envMap->w, h = envMap->h;
double theta = acos(r.y);
double phi = atan2(r.z, r.x) + PI;
double tx = phi / 2 / PI * w;
double ty = theta / PI * h;

int x[2], y[2];
x[0] = round(tx) - 1;
x[1] = x[0] + 1;
x[0] = clamp(x[0], 0, w - 1);
x[1] = clamp(x[1], 0, w - 1);
double dx = tx - x[0] - 0.5;
y[0] = round(ty) - 1;
y[1] = y[0] + 1;
y[0] = clamp(y[0], 0, h - 1);
y[1] = clamp(y[1], 0, h - 1);
double dy = ty - y[0] - 0.5;

Spectrum mix(0, 0, 0);
for (int i = 0; i < 2; i++) {
for (int j = 0; j < 2; j++) {
int id = x[i] + y[j] * w;
mix += envMap->data[id] * (i * dx + (1 - i) * (1 - dx)) * (j * dy + (1 - j) * (1 - dy));
}
}

return mix;
}

其他(随便口胡)

感觉这份代码的 trace_ray 在相对坐标中做很舒服的样子。。。

写完之后更想学 PBR 了,但不知道会花多少时间,毕竟自己的英文阅读速度不是很快,以及不知道会让自己离想象中的「游戏工程师」是更近还是更远。有看到 CMU 15-466 的 Computer Game Programming,也很想学的样子,然后再这么一想,就发现自己还有很多想学的东西,就会觉得时间好少的感觉。。。尤其是看到校内还要学一些奇怪的、占时间的课。。。