1 #include <algorithm>
2 #include <cassert>
3 #include <climits>
4 #include <cmath>
5 #include <ctime>
6 #include <iostream>
7 #include <vector>
8
9 #define NOMINMAX
10 #include "tiny_obj_loader.h"
11
12 #define STB_IMAGE_WRITE_IMPLEMENTATION
13 #include "stb_image_write.h"
14
15 #define TINYEXR_IMPLEMENTATION
16 #include "tinyexr.h"
17
18 #include "nanort.h"
19
20 #define USE_MULTIHIT_RAY_TRAVERSAL (0)
21
22 #ifndef M_PI
23 #define M_PI 3.141592683
24 #endif
25
26 const int uMaxBounces = 10;
27 // const int SPP = 10000;
28 const int SPP = 100;
29
30 #ifdef _WIN32
31 #ifdef __cplusplus
32 extern "C" {
33 #endif
34 #include <windows.h>
35 #ifdef __cplusplus
36 }
37 #endif
38 #pragma comment(lib, "winmm.lib")
39 #else
40 #if defined(__unix__) || defined(__APPLE__)
41 #include <sys/time.h>
42 #else
43 #include <ctime>
44 #endif
45 #endif
46
47 namespace {
48
49 // This class is NOT thread-safe timer!
50 class timerutil {
51 public:
52 #ifdef _WIN32
53 typedef DWORD time_t;
54
timerutil()55 timerutil() { ::timeBeginPeriod(1); }
~timerutil()56 ~timerutil() { ::timeEndPeriod(1); }
57
start()58 void start() { t_[0] = ::timeGetTime(); }
end()59 void end() { t_[1] = ::timeGetTime(); }
60
sec()61 time_t sec() { return (time_t)((t_[1] - t_[0]) / 1000); }
msec()62 time_t msec() { return (time_t)((t_[1] - t_[0])); }
usec()63 time_t usec() { return (time_t)((t_[1] - t_[0]) * 1000); }
current()64 time_t current() { return ::timeGetTime(); }
65
66 #else
67 #if defined(__unix__) || defined(__APPLE__)
68 typedef unsigned long int time_t;
69
70 void start() { gettimeofday(tv + 0, &tz); }
71 void end() { gettimeofday(tv + 1, &tz); }
72
73 time_t sec() { return (time_t)(tv[1].tv_sec - tv[0].tv_sec); }
74 time_t msec() {
75 return this->sec() * 1000 +
76 (time_t)((tv[1].tv_usec - tv[0].tv_usec) / 1000);
77 }
78 time_t usec() {
79 return this->sec() * 1000000 + (time_t)(tv[1].tv_usec - tv[0].tv_usec);
80 }
81 time_t current() {
82 struct timeval t;
83 gettimeofday(&t, NULL);
84 return (time_t)(t.tv_sec * 1000 + t.tv_usec);
85 }
86
87 #else // C timer
88 // using namespace std;
89 typedef clock_t time_t;
90
91 void start() { t_[0] = clock(); }
92 void end() { t_[1] = clock(); }
93
94 time_t sec() { return (time_t)((t_[1] - t_[0]) / CLOCKS_PER_SEC); }
95 time_t msec() { return (time_t)((t_[1] - t_[0]) * 1000 / CLOCKS_PER_SEC); }
96 time_t usec() { return (time_t)((t_[1] - t_[0]) * 1000000 / CLOCKS_PER_SEC); }
97 time_t current() { return (time_t)clock(); }
98
99 #endif
100 #endif
101
102 private:
103 #ifdef _WIN32
104 DWORD t_[2];
105 #else
106 #if defined(__unix__) || defined(__APPLE__)
107 struct timeval tv[2];
108 struct timezone tz;
109 #else
110 time_t t_[2];
111 #endif
112 #endif
113 };
114
115 struct float3 {
float3__anone3d206b80111::float3116 float3() {}
float3__anone3d206b80111::float3117 float3(float xx, float yy, float zz) {
118 x = xx;
119 y = yy;
120 z = zz;
121 }
float3__anone3d206b80111::float3122 float3(const float *p) {
123 x = p[0];
124 y = p[1];
125 z = p[2];
126 }
127
operator *__anone3d206b80111::float3128 float3 operator*(float f) const { return float3(x * f, y * f, z * f); }
operator -__anone3d206b80111::float3129 float3 operator-(const float3 &f2) const {
130 return float3(x - f2.x, y - f2.y, z - f2.z);
131 }
operator -__anone3d206b80111::float3132 float3 operator-() const { return float3(-x, -y, -z); }
operator *__anone3d206b80111::float3133 float3 operator*(const float3 &f2) const {
134 return float3(x * f2.x, y * f2.y, z * f2.z);
135 }
operator +__anone3d206b80111::float3136 float3 operator+(const float3 &f2) const {
137 return float3(x + f2.x, y + f2.y, z + f2.z);
138 }
operator +=__anone3d206b80111::float3139 float3 &operator+=(const float3 &f2) {
140 x += f2.x;
141 y += f2.y;
142 z += f2.z;
143 return (*this);
144 }
operator *=__anone3d206b80111::float3145 float3 &operator*=(const float3 &f2) {
146 x *= f2.x;
147 y *= f2.y;
148 z *= f2.z;
149 return (*this);
150 }
operator *=__anone3d206b80111::float3151 float3 &operator*=(const float &f2) {
152 x *= f2;
153 y *= f2;
154 z *= f2;
155 return (*this);
156 }
operator /__anone3d206b80111::float3157 float3 operator/(const float3 &f2) const {
158 return float3(x / f2.x, y / f2.y, z / f2.z);
159 }
operator /__anone3d206b80111::float3160 float3 operator/(const float &f2) const {
161 return float3(x / f2, y / f2, z / f2);
162 }
operator []__anone3d206b80111::float3163 float operator[](int i) const { return (&x)[i]; }
operator []__anone3d206b80111::float3164 float &operator[](int i) { return (&x)[i]; }
165
neg__anone3d206b80111::float3166 float3 neg() { return float3(-x, -y, -z); }
167
length__anone3d206b80111::float3168 float length() { return sqrtf(x * x + y * y + z * z); }
169
normalize__anone3d206b80111::float3170 void normalize() {
171 float len = length();
172 if (fabs(len) > 1.0e-6) {
173 float inv_len = 1.0 / len;
174 x *= inv_len;
175 y *= inv_len;
176 z *= inv_len;
177 }
178 }
179
180 float x, y, z;
181 // float pad; // for alignment
182 };
183
normalize(float3 v)184 inline float3 normalize(float3 v) {
185 v.normalize();
186 return v;
187 }
188
operator *(float f,const float3 & v)189 inline float3 operator*(float f, const float3 &v) {
190 return float3(v.x * f, v.y * f, v.z * f);
191 }
192
vcross(float3 a,float3 b)193 inline float3 vcross(float3 a, float3 b) {
194 float3 c;
195 c[0] = a[1] * b[2] - a[2] * b[1];
196 c[1] = a[2] * b[0] - a[0] * b[2];
197 c[2] = a[0] * b[1] - a[1] * b[0];
198 return c;
199 }
200
vdot(float3 a,float3 b)201 inline float vdot(float3 a, float3 b) {
202 return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
203 }
204
uniformFloat(float min,float max)205 float uniformFloat(float min, float max) {
206 return min + float(rand()) / RAND_MAX * (max - min);
207 }
208
209 // Building an Orthonormal Basis, Revisited
210 // http://jcgt.org/published/0006/01/01/
revisedONB(const float3 & n,float3 & b1,float3 & b2)211 void revisedONB(const float3 &n, float3 &b1, float3 &b2) {
212 if (n.z < 0.0f) {
213 const float a = 1.0f / (1.0f - n.z);
214 const float b = n.x * n.y * a;
215 b1 = float3(1.0f - n.x * n.x * a, -b, n.x);
216 b2 = float3(b, n.y * n.y * a - 1.0f, -n.y);
217 } else {
218 const float a = 1.0f / (1.0f + n.z);
219 const float b = -n.x * n.y * a;
220 b1 = float3(1.0f - n.x * n.x * a, b, -n.x);
221 b2 = float3(b, 1.0f - n.y * n.y * a, -n.y);
222 }
223 }
224
directionCosTheta(float3 normal)225 float3 directionCosTheta(float3 normal) {
226 float u1 = uniformFloat(0, 1);
227 float phi = uniformFloat(0, 2 * M_PI);
228
229 float r = sqrt(u1);
230
231 float x = r * cosf(phi);
232 float y = r * sinf(phi);
233 float z = sqrtf(1.0 - u1);
234
235 #if 0 // simpler way
236 float3 xDir =
237 fabsf(normal.x) < fabsf(normal.y) ? float3(1, 0, 0) : float3(0, 1, 0);
238 float3 yDir = normalize(vcross(normal, xDir));
239 xDir = vcross(yDir, normal);
240 #else // better way
241 float3 xDir, yDir;
242 revisedONB(normal, xDir, yDir);
243 #endif
244 return xDir * x + yDir * y + z * normal;
245 }
246
PdfAtoW(const float aPdfA,const float aDist,const float aCosThere)247 inline float PdfAtoW(const float aPdfA, const float aDist,
248 const float aCosThere) {
249 return aPdfA * (aDist * aDist) / std::abs(aCosThere);
250 }
251
252 typedef struct {
253 size_t num_vertices;
254 size_t num_faces;
255 float *vertices; /// [xyz] * num_vertices
256 float *facevarying_normals; /// [xyz] * 3(triangle) * num_faces
257 float *facevarying_tangents; /// [xyz] * 3(triangle) * num_faces
258 float *facevarying_binormals; /// [xyz] * 3(triangle) * num_faces
259 float *facevarying_uvs; /// [xyz] * 3(triangle) * num_faces
260 float *facevarying_vertex_colors; /// [xyz] * 3(triangle) * num_faces
261 unsigned int *faces; /// triangle x num_faces
262 unsigned int *material_ids; /// index x num_faces
263 } Mesh;
264
265 struct Material {
266 float ambient[3];
267 float diffuse[3];
268 float reflection[3];
269 float refraction[3];
270 int id;
271 int diffuse_texid;
272 int reflection_texid;
273 int transparency_texid;
274 int bump_texid;
275 int normal_texid; // normal map
276 int alpha_texid; // alpha map
277
Material__anone3d206b80111::Material278 Material() {
279 ambient[0] = 0.0;
280 ambient[1] = 0.0;
281 ambient[2] = 0.0;
282 diffuse[0] = 0.5;
283 diffuse[1] = 0.5;
284 diffuse[2] = 0.5;
285 reflection[0] = 0.0;
286 reflection[1] = 0.0;
287 reflection[2] = 0.0;
288 refraction[0] = 0.0;
289 refraction[1] = 0.0;
290 refraction[2] = 0.0;
291 id = -1;
292 diffuse_texid = -1;
293 reflection_texid = -1;
294 transparency_texid = -1;
295 bump_texid = -1;
296 normal_texid = -1;
297 alpha_texid = -1;
298 }
299 };
300
calcNormal(float3 & N,float3 v0,float3 v1,float3 v2)301 void calcNormal(float3 &N, float3 v0, float3 v1, float3 v2) {
302 float3 v10 = v1 - v0;
303 float3 v20 = v2 - v0;
304
305 N = vcross(v20, v10);
306 N.normalize();
307 }
308
309 class EmissiveFace {
310 public:
EmissiveFace(unsigned int f=0,unsigned int m=0)311 EmissiveFace(unsigned int f = 0, unsigned int m = 0) : face_(f), mtl_(m) {}
312 unsigned int face_;
313 unsigned int mtl_;
314 };
315
316 class MeshLight {
317 public:
MeshLight(const Mesh & mesh,const std::vector<tinyobj::material_t> & materials)318 MeshLight(const Mesh &mesh, const std::vector<tinyobj::material_t> &materials)
319 : mesh_(mesh), materials_(materials) {
320 for (unsigned int face = 0; face < mesh_.num_faces; face++) {
321 unsigned int mtl_id = mesh_.material_ids[face];
322 const tinyobj::material_t &faceMtl = materials_[mtl_id];
323
324 if (faceMtl.emission[0] > 0.0f || faceMtl.emission[1] > 0.0f ||
325 faceMtl.emission[2] > 0.0f) {
326 EmissiveFace ef(face, mtl_id);
327 emissive_faces_.push_back(ef);
328 }
329 }
330 }
331
sampleDirect(const float3 & x,float Xi1,float Xi2,float3 & dstDir,float & dstDist,float & dstPdf,float3 & dstRadiance) const332 void sampleDirect(const float3 &x, float Xi1, float Xi2, float3 &dstDir,
333 float &dstDist, float &dstPdf, float3 &dstRadiance) const {
334 unsigned int num_faces = emissive_faces_.size();
335 unsigned int face = std::min(
336 static_cast<unsigned int>(floor(Xi1 * num_faces)), num_faces - 1);
337 float lightPickPdf = 1.0f / float(num_faces);
338
339 // normalize random number
340 Xi1 = Xi1 * num_faces - face;
341
342 unsigned int fid = emissive_faces_[face].face_;
343 unsigned int mtlid = emissive_faces_[face].mtl_;
344
345 unsigned int f0 = mesh_.faces[3 * fid + 0];
346 unsigned int f1 = mesh_.faces[3 * fid + 1];
347 unsigned int f2 = mesh_.faces[3 * fid + 2];
348
349 float3 v0, v1, v2;
350
351 v0[0] = mesh_.vertices[3 * f0 + 0];
352 v0[1] = mesh_.vertices[3 * f0 + 1];
353 v0[2] = mesh_.vertices[3 * f0 + 2];
354
355 v1[0] = mesh_.vertices[3 * f1 + 0];
356 v1[1] = mesh_.vertices[3 * f1 + 1];
357 v1[2] = mesh_.vertices[3 * f1 + 2];
358
359 v2[0] = mesh_.vertices[3 * f2 + 0];
360 v2[1] = mesh_.vertices[3 * f2 + 1];
361 v2[2] = mesh_.vertices[3 * f2 + 2];
362
363 float Xi1_ = std::sqrt(Xi1);
364 float c0 = 1.0f - Xi1_;
365 float c1 = Xi1_ * (1.0 - Xi2);
366 float c2 = Xi1_ * Xi2;
367
368 float3 tmp = vcross(v1 - v0, v2 - v0);
369 float3 norm = normalize(tmp);
370 float area = tmp.length() / 2.0f;
371
372 dstPdf = 0.0;
373 float3 lp = c0 * v0 + c1 * v1 + c2 * v2;
374 float areaPdf = lightPickPdf * (1.0f / area);
375 dstDir = lp - x;
376 dstDist = dstDir.length();
377
378 float3 ll = materials_[mtlid].emission;
379 if (dstDist > 0.000001f) {
380 dstDir.normalize();
381 float cosAtLight = std::max(vdot(-dstDir, norm), 0.0f);
382 dstRadiance = ll * cosAtLight; // light has cosine edf
383
384 // convert pdf to solid angle measure
385 dstPdf = PdfAtoW(areaPdf, dstDist, cosAtLight);
386 }
387 }
388
389 std::vector<EmissiveFace> emissive_faces_;
390 const Mesh &mesh_;
391 const std::vector<tinyobj::material_t> &materials_;
392 };
393
394 // Save in RAW headerless format, for use when exr tools are not available in
395 // system
SaveImageRaw(const char * filename,const float * rgb,int width,int height)396 void SaveImageRaw(const char *filename, const float *rgb, int width,
397 int height) {
398 std::vector<unsigned char> rawbuf;
399 rawbuf.resize(3 * width * height);
400 unsigned char *raw = &rawbuf.at(0);
401
402 // @note { Apply gamma correction would be nice? }
403 for (int i = 0; i < width * height; i++) {
404 raw[i * 3] = (char)(rgb[3 * i + 0] * 255.0);
405 raw[i * 3 + 1] = (char)(rgb[3 * i + 1] * 255.0);
406 raw[i * 3 + 2] = (char)(rgb[3 * i + 2] * 255.0);
407 }
408 FILE *f = fopen(filename, "wb");
409 if (!f) {
410 printf("Error: Couldnt open output binary file %s\n", filename);
411 return;
412 }
413 fwrite(raw, 3 * width * height, 1, f);
414 fclose(f);
415 printf("Info: Saved RAW RGB image of [%dx%d] dimensions to [ %s ]\n", width,
416 height, filename);
417 }
418
SaveImagePNG(const char * filename,const float * rgb,int width,int height)419 void SaveImagePNG(const char *filename, const float *rgb, int width,
420 int height) {
421 unsigned char *bytes = new unsigned char[width * height * 3];
422 for (int y = 0; y < height; y++) {
423 for (int x = 0; x < width; x++) {
424 const int index = y * width + x;
425 bytes[index * 3 + 0] = (unsigned char)std::max(
426 0.0f, std::min(rgb[index * 3 + 0] * 255.0f, 255.0f));
427 bytes[index * 3 + 1] = (unsigned char)std::max(
428 0.0f, std::min(rgb[index * 3 + 1] * 255.0f, 255.0f));
429 bytes[index * 3 + 2] = (unsigned char)std::max(
430 0.0f, std::min(rgb[index * 3 + 2] * 255.0f, 255.0f));
431 }
432 }
433 stbi_write_png(filename, width, height, 3, bytes, width * 3);
434 delete[] bytes;
435 }
436
SaveImage(const char * filename,const float * rgb,int width,int height)437 void SaveImage(const char *filename, const float *rgb, int width, int height) {
438 const char *err = NULL;
439 int ret = SaveEXR(rgb, width, height, /* RGB */ 3, /* fp16 */0, filename, &err);
440 if (ret != TINYEXR_SUCCESS) {
441 if (err) {
442 fprintf(stderr, "EXR save error: %s(%d)\n", err, ret);
443 FreeEXRErrorMessage(err);
444 } else {
445 fprintf(stderr, "EXR save error: %d\n", ret);
446 }
447 } else {
448 printf("Saved image to [ %s ]\n", filename);
449 }
450 }
451
LoadObj(Mesh & mesh,std::vector<tinyobj::material_t> & materials,const char * filename,float scale,const char * mtl_path)452 bool LoadObj(Mesh &mesh, std::vector<tinyobj::material_t> &materials,
453 const char *filename, float scale, const char *mtl_path) {
454 std::vector<tinyobj::shape_t> shapes;
455
456 std::string err = tinyobj::LoadObj(shapes, materials, filename, mtl_path);
457
458 if (!err.empty()) {
459 std::cerr << err << std::endl;
460 return false;
461 }
462
463 std::cout << "[LoadOBJ] # of shapes in .obj : " << shapes.size() << std::endl;
464 std::cout << "[LoadOBJ] # of materials in .obj : " << materials.size()
465 << std::endl;
466
467 size_t num_vertices = 0;
468 size_t num_faces = 0;
469 for (size_t i = 0; i < shapes.size(); i++) {
470 printf(" shape[%ld].name = %s\n", i, shapes[i].name.c_str());
471 printf(" shape[%ld].indices: %ld\n", i, shapes[i].mesh.indices.size());
472 assert((shapes[i].mesh.indices.size() % 3) == 0);
473 printf(" shape[%ld].vertices: %ld\n", i, shapes[i].mesh.positions.size());
474 assert((shapes[i].mesh.positions.size() % 3) == 0);
475 printf(" shape[%ld].normals: %ld\n", i, shapes[i].mesh.normals.size());
476 assert((shapes[i].mesh.normals.size() % 3) == 0);
477
478 num_vertices += shapes[i].mesh.positions.size() / 3;
479 num_faces += shapes[i].mesh.indices.size() / 3;
480 }
481 std::cout << "[LoadOBJ] # of faces: " << num_faces << std::endl;
482 std::cout << "[LoadOBJ] # of vertices: " << num_vertices << std::endl;
483
484 // @todo { material and texture. }
485
486 // Shape -> Mesh
487 mesh.num_faces = num_faces;
488 mesh.num_vertices = num_vertices;
489 mesh.vertices = new float[num_vertices * 3];
490 mesh.faces = new unsigned int[num_faces * 3];
491 mesh.material_ids = new unsigned int[num_faces];
492 memset(mesh.material_ids, 0, sizeof(int) * num_faces);
493 mesh.facevarying_normals = new float[num_faces * 3 * 3];
494 mesh.facevarying_uvs = new float[num_faces * 3 * 2];
495 memset(mesh.facevarying_uvs, 0, sizeof(float) * 2 * 3 * num_faces);
496
497 // @todo {}
498 mesh.facevarying_tangents = NULL;
499 mesh.facevarying_binormals = NULL;
500
501 size_t vertexIdxOffset = 0;
502 size_t faceIdxOffset = 0;
503 for (size_t i = 0; i < shapes.size(); i++) {
504 for (size_t f = 0; f < shapes[i].mesh.indices.size() / 3; f++) {
505 mesh.faces[3 * (faceIdxOffset + f) + 0] =
506 shapes[i].mesh.indices[3 * f + 0];
507 mesh.faces[3 * (faceIdxOffset + f) + 1] =
508 shapes[i].mesh.indices[3 * f + 1];
509 mesh.faces[3 * (faceIdxOffset + f) + 2] =
510 shapes[i].mesh.indices[3 * f + 2];
511
512 mesh.faces[3 * (faceIdxOffset + f) + 0] += vertexIdxOffset;
513 mesh.faces[3 * (faceIdxOffset + f) + 1] += vertexIdxOffset;
514 mesh.faces[3 * (faceIdxOffset + f) + 2] += vertexIdxOffset;
515
516 mesh.material_ids[faceIdxOffset + f] = shapes[i].mesh.material_ids[f];
517 }
518
519 for (size_t v = 0; v < shapes[i].mesh.positions.size() / 3; v++) {
520 mesh.vertices[3 * (vertexIdxOffset + v) + 0] =
521 scale * shapes[i].mesh.positions[3 * v + 0];
522 mesh.vertices[3 * (vertexIdxOffset + v) + 1] =
523 scale * shapes[i].mesh.positions[3 * v + 1];
524 mesh.vertices[3 * (vertexIdxOffset + v) + 2] =
525 scale * shapes[i].mesh.positions[3 * v + 2];
526 }
527
528 if (shapes[i].mesh.normals.size() > 0) {
529 for (size_t f = 0; f < shapes[i].mesh.indices.size() / 3; f++) {
530 int f0, f1, f2;
531
532 f0 = shapes[i].mesh.indices[3 * f + 0];
533 f1 = shapes[i].mesh.indices[3 * f + 1];
534 f2 = shapes[i].mesh.indices[3 * f + 2];
535
536 float3 n0, n1, n2;
537
538 n0[0] = shapes[i].mesh.normals[3 * f0 + 0];
539 n0[1] = shapes[i].mesh.normals[3 * f0 + 1];
540 n0[2] = shapes[i].mesh.normals[3 * f0 + 2];
541
542 n1[0] = shapes[i].mesh.normals[3 * f1 + 0];
543 n1[1] = shapes[i].mesh.normals[3 * f1 + 1];
544 n1[2] = shapes[i].mesh.normals[3 * f1 + 2];
545
546 n2[0] = shapes[i].mesh.normals[3 * f2 + 0];
547 n2[1] = shapes[i].mesh.normals[3 * f2 + 1];
548 n2[2] = shapes[i].mesh.normals[3 * f2 + 2];
549
550 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 0) + 0] = n0[0];
551 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 0) + 1] = n0[1];
552 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 0) + 2] = n0[2];
553
554 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 1) + 0] = n1[0];
555 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 1) + 1] = n1[1];
556 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 1) + 2] = n1[2];
557
558 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 2) + 0] = n2[0];
559 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 2) + 1] = n2[1];
560 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 2) + 2] = n2[2];
561 }
562 } else {
563 // calc geometric normal
564 for (size_t f = 0; f < shapes[i].mesh.indices.size() / 3; f++) {
565 int f0, f1, f2;
566
567 f0 = shapes[i].mesh.indices[3 * f + 0];
568 f1 = shapes[i].mesh.indices[3 * f + 1];
569 f2 = shapes[i].mesh.indices[3 * f + 2];
570
571 float3 v0, v1, v2;
572
573 v0[0] = shapes[i].mesh.positions[3 * f0 + 0];
574 v0[1] = shapes[i].mesh.positions[3 * f0 + 1];
575 v0[2] = shapes[i].mesh.positions[3 * f0 + 2];
576
577 v1[0] = shapes[i].mesh.positions[3 * f1 + 0];
578 v1[1] = shapes[i].mesh.positions[3 * f1 + 1];
579 v1[2] = shapes[i].mesh.positions[3 * f1 + 2];
580
581 v2[0] = shapes[i].mesh.positions[3 * f2 + 0];
582 v2[1] = shapes[i].mesh.positions[3 * f2 + 1];
583 v2[2] = shapes[i].mesh.positions[3 * f2 + 2];
584
585 float3 N;
586 calcNormal(N, v0, v1, v2);
587
588 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 0) + 0] = N[0];
589 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 0) + 1] = N[1];
590 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 0) + 2] = N[2];
591
592 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 1) + 0] = N[0];
593 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 1) + 1] = N[1];
594 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 1) + 2] = N[2];
595
596 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 2) + 0] = N[0];
597 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 2) + 1] = N[1];
598 mesh.facevarying_normals[3 * (3 * (faceIdxOffset + f) + 2) + 2] = N[2];
599 }
600 }
601
602 if (shapes[i].mesh.texcoords.size() > 0) {
603 for (size_t f = 0; f < shapes[i].mesh.indices.size() / 3; f++) {
604 int f0, f1, f2;
605
606 f0 = shapes[i].mesh.indices[3 * f + 0];
607 f1 = shapes[i].mesh.indices[3 * f + 1];
608 f2 = shapes[i].mesh.indices[3 * f + 2];
609
610 float3 n0, n1, n2;
611
612 n0[0] = shapes[i].mesh.texcoords[2 * f0 + 0];
613 n0[1] = shapes[i].mesh.texcoords[2 * f0 + 1];
614
615 n1[0] = shapes[i].mesh.texcoords[2 * f1 + 0];
616 n1[1] = shapes[i].mesh.texcoords[2 * f1 + 1];
617
618 n2[0] = shapes[i].mesh.texcoords[2 * f2 + 0];
619 n2[1] = shapes[i].mesh.texcoords[2 * f2 + 1];
620
621 mesh.facevarying_uvs[2 * (3 * (faceIdxOffset + f) + 0) + 0] = n0[0];
622 mesh.facevarying_uvs[2 * (3 * (faceIdxOffset + f) + 0) + 1] = n0[1];
623
624 mesh.facevarying_uvs[2 * (3 * (faceIdxOffset + f) + 1) + 0] = n1[0];
625 mesh.facevarying_uvs[2 * (3 * (faceIdxOffset + f) + 1) + 1] = n1[1];
626
627 mesh.facevarying_uvs[2 * (3 * (faceIdxOffset + f) + 2) + 0] = n2[0];
628 mesh.facevarying_uvs[2 * (3 * (faceIdxOffset + f) + 2) + 1] = n2[1];
629 }
630 }
631
632 vertexIdxOffset += shapes[i].mesh.positions.size() / 3;
633 faceIdxOffset += shapes[i].mesh.indices.size() / 3;
634 }
635
636 return true;
637 }
638 } // namespace
639
sign(float f)640 inline float sign(float f) { return f < 0 ? -1 : 1; }
641
reflect(float3 I,float3 N)642 inline float3 reflect(float3 I, float3 N) { return I - 2 * vdot(I, N) * N; }
643
refract(float3 I,float3 N,float eta)644 inline float3 refract(float3 I, float3 N, float eta) {
645 float NdotI = vdot(N, I);
646 float k = 1.0f - eta * eta * (1.0f - NdotI * NdotI);
647 if (k < 0.0f)
648 return float3(0, 0, 0);
649 else
650 return eta * I - (eta * NdotI + sqrtf(k)) * N;
651 }
652
pow5(float val)653 inline float pow5(float val) { return val * val * val * val * val; }
654
fresnel_schlick(float3 H,float3 norm,float n1)655 inline float fresnel_schlick(float3 H, float3 norm, float n1) {
656 float r0 = n1 * n1;
657 return r0 + (1 - r0) * pow5(1 - vdot(H, norm));
658 }
659
progressBar(int tick,int total,int width=50)660 void progressBar(int tick, int total, int width = 50) {
661 float ratio = 100.0f * tick / total;
662 float count = width * tick / total;
663 std::string bar(width, ' ');
664 std::fill(bar.begin(), bar.begin() + count, '+');
665 printf("[ %6.2f %% ] [ %s ]%c", ratio, bar.c_str(),
666 tick == total ? '\n' : '\r');
667 std::fflush(stdout);
668 }
669
CheckForOccluder(float3 p1,float3 p2,const Mesh & mesh,const nanort::BVHAccel<float> & accel)670 bool CheckForOccluder(float3 p1, float3 p2, const Mesh &mesh,
671 const nanort::BVHAccel<float> &accel) {
672 static const float ray_eps = 0.00001f;
673
674 float3 dir = p2 - p1;
675 float dist = dir.length();
676 dir.normalize();
677
678 nanort::Ray<float> shadow_ray;
679 shadow_ray.min_t = ray_eps;
680 shadow_ray.max_t = dist - ray_eps;
681 shadow_ray.dir[0] = dir[0];
682 shadow_ray.dir[1] = dir[1];
683 shadow_ray.dir[2] = dir[2];
684 shadow_ray.org[0] = p1[0];
685 shadow_ray.org[1] = p1[1];
686 shadow_ray.org[2] = p1[2];
687
688 nanort::TriangleIntersector<> triangle_intersector(mesh.vertices, mesh.faces,
689 sizeof(float) * 3);
690 nanort::TriangleIntersection<> isect;
691 if (!accel.Traverse(shadow_ray, triangle_intersector, &isect)) {
692 return false;
693 }
694
695 return true;
696 }
697
main(int argc,char ** argv)698 int main(int argc, char **argv) {
699 int width = 512;
700 int height = 512;
701
702 float scale = 1.0f;
703
704 std::string objFilename = "../common/cornellbox_suzanne_lucy.obj";
705 std::string mtlPath = "../common/";
706
707 if (argc > 1) {
708 objFilename = std::string(argv[1]);
709 }
710
711 if (argc > 2) {
712 scale = atof(argv[2]);
713 }
714
715 if (argc > 3) {
716 mtlPath = std::string(argv[3]);
717 }
718
719 #ifdef _OPENMP
720 printf("Using OpenMP: yes!\n");
721 #else
722 printf("Using OpenMP: no!\n");
723 #endif
724
725 bool ret = false;
726
727 Mesh mesh;
728 std::vector<tinyobj::material_t> materials;
729 ret = LoadObj(mesh, materials, objFilename.c_str(), scale, mtlPath.c_str());
730 if (!ret) {
731 fprintf(stderr, "Failed to load [ %s ]\n", objFilename.c_str());
732 return -1;
733 }
734
735 MeshLight lights(mesh, materials);
736
737 nanort::BVHBuildOptions<float> build_options; // Use default option
738 build_options.cache_bbox = false;
739
740 printf(" BVH build option:\n");
741 printf(" # of leaf primitives: %d\n", build_options.min_leaf_primitives);
742 printf(" SAH binsize : %d\n", build_options.bin_size);
743
744 timerutil t;
745 t.start();
746
747 nanort::TriangleMesh<float> triangle_mesh(mesh.vertices, mesh.faces,
748 sizeof(float) * 3);
749 nanort::TriangleSAHPred<float> triangle_pred(mesh.vertices, mesh.faces,
750 sizeof(float) * 3);
751
752 printf("num_triangles = %lu\n", mesh.num_faces);
753 printf("faces = %p\n", mesh.faces);
754
755 nanort::BVHAccel<float> accel;
756 ret =
757 accel.Build(mesh.num_faces, triangle_mesh, triangle_pred, build_options);
758 assert(ret);
759
760 t.end();
761 printf(" BVH build time: %f secs\n", t.msec() / 1000.0);
762
763 nanort::BVHBuildStatistics stats = accel.GetStatistics();
764
765 printf(" BVH statistics:\n");
766 printf(" # of leaf nodes: %d\n", stats.num_leaf_nodes);
767 printf(" # of branch nodes: %d\n", stats.num_branch_nodes);
768 printf(" Max tree depth : %d\n", stats.max_tree_depth);
769 float bmin[3], bmax[3];
770 accel.BoundingBox(bmin, bmax);
771 printf(" Bmin : %f, %f, %f\n", bmin[0], bmin[1], bmin[2]);
772 printf(" Bmax : %f, %f, %f\n", bmax[0], bmax[1], bmax[2]);
773
774 std::vector<float> rgb(width * height * 3, 0.0f);
775
776 srand(0);
777
778 // Shoot rays.
779 #ifdef _OPENMP
780 #pragma omp parallel for schedule(dynamic, 1)
781 #endif
782 for (int y = 0; y < height; y++) {
783 for (int x = 0; x < width; x++) {
784 float3 finalColor = float3(0, 0, 0);
785 for (int i = 0; i < SPP; ++i) {
786 float px = x + uniformFloat(-0.5, 0.5);
787 float py = y + uniformFloat(-0.5, 0.5);
788 // Simple camera. change eye pos and direction fit to .obj model.
789
790 float3 rayDir = float3((px / (float)width) - 0.5f,
791 (py / (float)height) - 0.5f, -1.0f);
792 rayDir.normalize();
793
794 float3 rayOrg = float3(0.0f, 5.0f, 20.0f);
795
796 float3 color = float3(0, 0, 0);
797 float3 weight = float3(1, 1, 1);
798
799 int b;
800
801 bool do_emmition = true; // just skit emmition if light sampling was
802 // done on previous event (No MIS)
803 for (b = 0; b < uMaxBounces; ++b) {
804 // Russian Roulette
805 float rr_fac = 1.0f;
806 if (b > 3) {
807 float rr_rand = uniformFloat(0, 1);
808 float termination_probability = 0.2f;
809 if (rr_rand < termination_probability) {
810 break;
811 }
812 rr_fac = 1.0 - termination_probability;
813 }
814 weight *= 1.0 / rr_fac;
815
816 nanort::Ray<float> ray;
817 float kFar = 1.0e+30f;
818 ray.min_t = 0.001f;
819 ray.max_t = kFar;
820
821 ray.dir[0] = rayDir[0];
822 ray.dir[1] = rayDir[1];
823 ray.dir[2] = rayDir[2];
824 ray.org[0] = rayOrg[0];
825 ray.org[1] = rayOrg[1];
826 ray.org[2] = rayOrg[2];
827
828 nanort::TriangleIntersector<> triangle_intersector(
829 mesh.vertices, mesh.faces, sizeof(float) * 3);
830 nanort::TriangleIntersection<> isect;
831 bool hit = accel.Traverse(ray, triangle_intersector, &isect);
832
833 if (!hit) {
834 break;
835 }
836
837 rayOrg += rayDir * isect.t;
838
839 unsigned int fid = isect.prim_id;
840 float3 norm(0, 0, 0);
841 if (mesh.facevarying_normals) {
842 float3 normals[3];
843 for (int vId = 0; vId < 3; vId++) {
844 normals[vId][0] = mesh.facevarying_normals[9 * fid + 3 * vId + 0];
845 normals[vId][1] = mesh.facevarying_normals[9 * fid + 3 * vId + 1];
846 normals[vId][2] = mesh.facevarying_normals[9 * fid + 3 * vId + 2];
847 }
848 float u = isect.u;
849 float v = isect.v;
850 norm = (1.0 - u - v) * normals[0] + u * normals[1] + v * normals[2];
851 norm.normalize();
852 }
853
854 // Flip normal torwards incoming ray for backface shading
855 float3 originalNorm = norm;
856 if (vdot(norm, rayDir) > 0) {
857 norm *= -1;
858 }
859
860 // Get properties from the material of the hit primitive
861 unsigned int matId = mesh.material_ids[fid];
862 tinyobj::material_t mat = materials[matId];
863
864 float3 diffuseColor(mat.diffuse);
865 float3 emissiveColor(mat.emission);
866 float3 specularColor(mat.specular);
867 float3 refractionColor(mat.transmittance);
868 float ior = mat.ior;
869
870 // Calculate fresnel factor based on ior.
871 float inside =
872 sign(vdot(rayDir, originalNorm)); // 1 for inside, -1 for outside
873 // Assume ior of medium outside of objects = 1.0
874 float n1 = inside < 0 ? 1.0 / ior : ior;
875 float n2 = 1.0 / n1;
876
877 float fresnel = fresnel_schlick(-rayDir, norm, (n1 - n2) / (n1 + n2));
878
879 // Compute probabilities for each surface interaction.
880 // Specular is just regular reflectiveness * fresnel.
881 float rhoS = vdot(float3(1, 1, 1) / 3.0f, specularColor) * fresnel;
882 // If we don't have a specular reflection, choose either diffuse or
883 // transmissive
884 // Mix them based on the dissolve value of the material
885 float rhoD = vdot(float3(1, 1, 1) / 3.0f, diffuseColor) *
886 (1.0 - fresnel) * (1.0 - mat.dissolve);
887 float rhoR = vdot(float3(1, 1, 1) / 3.0f, refractionColor) *
888 (1.0 - fresnel) * mat.dissolve;
889
890 float rhoE = vdot(float3(1, 1, 1) / 3.0f, emissiveColor);
891
892 // Normalize probabilities so they sum to 1.0
893 float totalrho = rhoS + rhoD + rhoR + rhoE;
894 // No scattering event is likely, just stop here
895 if (totalrho < 0.0001) {
896 break;
897 }
898
899 rhoS /= totalrho;
900 rhoD /= totalrho;
901 rhoR /= totalrho;
902 rhoE /= totalrho;
903
904 // Choose an interaction based on the calculated probabilities
905 float rand = uniformFloat(0, 1);
906 float3 outDir;
907 // REFLECT glossy
908 if (rand < rhoS) {
909 outDir = reflect(rayDir, norm);
910 weight *= specularColor;
911 do_emmition = true;
912 }
913 // REFLECT diffuse
914 else if (rand < rhoS + rhoD) {
915 float3 brdfEval = (1.0f / M_PI) * diffuseColor;
916 float3 dl = float3(0.0, 0.0, 0.0), ldir, ll;
917 float lpdf, ldist;
918 lights.sampleDirect(rayOrg, uniformFloat(0, 1), uniformFloat(0, 1),
919 ldir, ldist, lpdf, ll);
920
921 if (lpdf > 0.0f) {
922 float cosTheta = std::abs(vdot(ldir, norm));
923 float3 directLight = (brdfEval * ll * cosTheta) / lpdf;
924 bool visible =
925 !CheckForOccluder(rayOrg, rayOrg + ldir * ldist, mesh, accel);
926
927 color += directLight * visible * weight;
928 }
929
930 // Sample cosine weighted hemisphere
931 outDir = directionCosTheta(norm);
932 weight *= diffuseColor;
933 do_emmition = false;
934 }
935 // REFRACT
936 else if (rand < rhoD + rhoS + rhoR) {
937 outDir = refract(rayDir, -inside * originalNorm, n1);
938 weight *= refractionColor;
939 do_emmition = true;
940 }
941 // EMIT
942 else {
943 // Weight light by cosine factor (surface emits most light in normal
944 // direction)
945 if (do_emmition) {
946 color += std::max(vdot(originalNorm, -rayDir), 0.0f) *
947 emissiveColor * weight;
948 }
949 break;
950 }
951
952 // Calculate new ray start position and set outgoing direction.
953 rayDir = outDir;
954 }
955
956 finalColor += color;
957 }
958
959 finalColor *= 1.0 / SPP;
960
961 // Gamme Correct
962 finalColor[0] = pow(finalColor[0], 1.0 / 2.2);
963 finalColor[1] = pow(finalColor[1], 1.0 / 2.2);
964 finalColor[2] = pow(finalColor[2], 1.0 / 2.2);
965
966 rgb[3 * ((height - y - 1) * width + x) + 0] = finalColor[0];
967 rgb[3 * ((height - y - 1) * width + x) + 1] = finalColor[1];
968 rgb[3 * ((height - y - 1) * width + x) + 2] = finalColor[2];
969 }
970
971 progressBar(y + 1, height);
972 }
973
974 // Save image.
975 SaveImage("render.exr", &rgb.at(0), width, height);
976 // Save Raw Image that can be opened by tools like GIMP
977 SaveImageRaw("render.data", &rgb.at(0), width, height);
978 SaveImagePNG("render.png", &rgb.at(0), width, height);
979
980 return 0;
981 }
982