1 /*
2 The MIT License (MIT)
3
4 Copyright (c) 2017 Light Transport Entertainment, Inc.
5
6 Permission is hereby granted, free of charge, to any person obtaining a copy
7 of this software and associated documentation files (the "Software"), to deal
8 in the Software without restriction, including without limitation the rights
9 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10 copies of the Software, and to permit persons to whom the Software is
11 furnished to do so, subject to the following conditions:
12
13 The above copyright notice and this permission notice shall be included in
14 all copies or substantial portions of the Software.
15
16 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
22 THE SOFTWARE.
23 */
24
25 #ifndef NANOSG_H_
26 #define NANOSG_H_
27
28 #ifdef _WIN32
29 #ifndef NOMINMAX
30 #define NOMINMAX
31 #endif
32 #endif
33
34 #include <iostream>
35 #include <limits>
36 #include <vector>
37
38 #include "nanort.h"
39
40 namespace nanosg {
41
42 template <class T>
43 class PrimitiveInterface;
44
45 template <class T>
46 class PrimitiveInterface {
47 public:
print()48 void print() { static_cast<T &>(this)->print(); }
49 };
50
51 class SpherePrimitive : PrimitiveInterface<SpherePrimitive> {
52 public:
print()53 void print() { std::cout << "Sphere" << std::endl; }
54 };
55
56 // 4x4 matrix
57 template <typename T>
58 class Matrix {
59 public:
60 Matrix();
61 ~Matrix();
62
Print(const T m[4][4])63 static void Print(const T m[4][4]) {
64 for (int i = 0; i < 4; i++) {
65 printf("m[%d] = %f, %f, %f, %f\n", i, m[i][0], m[i][1], m[i][2], m[i][3]);
66 }
67 }
68
Identity(T m[4][4])69 static void Identity(T m[4][4]) {
70 m[0][0] = static_cast<T>(1);
71 m[0][1] = static_cast<T>(0);
72 m[0][2] = static_cast<T>(0);
73 m[0][3] = static_cast<T>(0);
74 m[1][0] = static_cast<T>(0);
75 m[1][1] = static_cast<T>(1);
76 m[1][2] = static_cast<T>(0);
77 m[1][3] = static_cast<T>(0);
78 m[2][0] = static_cast<T>(0);
79 m[2][1] = static_cast<T>(0);
80 m[2][2] = static_cast<T>(1);
81 m[2][3] = static_cast<T>(0);
82 m[3][0] = static_cast<T>(0);
83 m[3][1] = static_cast<T>(0);
84 m[3][2] = static_cast<T>(0);
85 m[3][3] = static_cast<T>(1);
86 }
87
Copy(T dst[4][4],const T src[4][4])88 static void Copy(T dst[4][4], const T src[4][4]) {
89 memcpy(dst, src, sizeof(T) * 16);
90 }
91
Inverse(T m[4][4])92 static void Inverse(T m[4][4]) {
93 /*
94 * codes from intel web
95 * cramer's rule version
96 */
97 int i, j;
98 T tmp[12]; /* tmp array for pairs */
99 T tsrc[16]; /* array of transpose source matrix */
100 T det; /* determinant */
101
102 /* transpose matrix */
103 for (i = 0; i < 4; i++) {
104 tsrc[i] = m[i][0];
105 tsrc[i + 4] = m[i][1];
106 tsrc[i + 8] = m[i][2];
107 tsrc[i + 12] = m[i][3];
108 }
109
110 /* calculate pair for first 8 elements(cofactors) */
111 tmp[0] = tsrc[10] * tsrc[15];
112 tmp[1] = tsrc[11] * tsrc[14];
113 tmp[2] = tsrc[9] * tsrc[15];
114 tmp[3] = tsrc[11] * tsrc[13];
115 tmp[4] = tsrc[9] * tsrc[14];
116 tmp[5] = tsrc[10] * tsrc[13];
117 tmp[6] = tsrc[8] * tsrc[15];
118 tmp[7] = tsrc[11] * tsrc[12];
119 tmp[8] = tsrc[8] * tsrc[14];
120 tmp[9] = tsrc[10] * tsrc[12];
121 tmp[10] = tsrc[8] * tsrc[13];
122 tmp[11] = tsrc[9] * tsrc[12];
123
124 /* calculate first 8 elements(cofactors) */
125 m[0][0] = tmp[0] * tsrc[5] + tmp[3] * tsrc[6] + tmp[4] * tsrc[7];
126 m[0][0] -= tmp[1] * tsrc[5] + tmp[2] * tsrc[6] + tmp[5] * tsrc[7];
127 m[0][1] = tmp[1] * tsrc[4] + tmp[6] * tsrc[6] + tmp[9] * tsrc[7];
128 m[0][1] -= tmp[0] * tsrc[4] + tmp[7] * tsrc[6] + tmp[8] * tsrc[7];
129 m[0][2] = tmp[2] * tsrc[4] + tmp[7] * tsrc[5] + tmp[10] * tsrc[7];
130 m[0][2] -= tmp[3] * tsrc[4] + tmp[6] * tsrc[5] + tmp[11] * tsrc[7];
131 m[0][3] = tmp[5] * tsrc[4] + tmp[8] * tsrc[5] + tmp[11] * tsrc[6];
132 m[0][3] -= tmp[4] * tsrc[4] + tmp[9] * tsrc[5] + tmp[10] * tsrc[6];
133 m[1][0] = tmp[1] * tsrc[1] + tmp[2] * tsrc[2] + tmp[5] * tsrc[3];
134 m[1][0] -= tmp[0] * tsrc[1] + tmp[3] * tsrc[2] + tmp[4] * tsrc[3];
135 m[1][1] = tmp[0] * tsrc[0] + tmp[7] * tsrc[2] + tmp[8] * tsrc[3];
136 m[1][1] -= tmp[1] * tsrc[0] + tmp[6] * tsrc[2] + tmp[9] * tsrc[3];
137 m[1][2] = tmp[3] * tsrc[0] + tmp[6] * tsrc[1] + tmp[11] * tsrc[3];
138 m[1][2] -= tmp[2] * tsrc[0] + tmp[7] * tsrc[1] + tmp[10] * tsrc[3];
139 m[1][3] = tmp[4] * tsrc[0] + tmp[9] * tsrc[1] + tmp[10] * tsrc[2];
140 m[1][3] -= tmp[5] * tsrc[0] + tmp[8] * tsrc[1] + tmp[11] * tsrc[2];
141
142 /* calculate pairs for second 8 elements(cofactors) */
143 tmp[0] = tsrc[2] * tsrc[7];
144 tmp[1] = tsrc[3] * tsrc[6];
145 tmp[2] = tsrc[1] * tsrc[7];
146 tmp[3] = tsrc[3] * tsrc[5];
147 tmp[4] = tsrc[1] * tsrc[6];
148 tmp[5] = tsrc[2] * tsrc[5];
149 tmp[6] = tsrc[0] * tsrc[7];
150 tmp[7] = tsrc[3] * tsrc[4];
151 tmp[8] = tsrc[0] * tsrc[6];
152 tmp[9] = tsrc[2] * tsrc[4];
153 tmp[10] = tsrc[0] * tsrc[5];
154 tmp[11] = tsrc[1] * tsrc[4];
155
156 /* calculate second 8 elements(cofactors) */
157 m[2][0] = tmp[0] * tsrc[13] + tmp[3] * tsrc[14] + tmp[4] * tsrc[15];
158 m[2][0] -= tmp[1] * tsrc[13] + tmp[2] * tsrc[14] + tmp[5] * tsrc[15];
159 m[2][1] = tmp[1] * tsrc[12] + tmp[6] * tsrc[14] + tmp[9] * tsrc[15];
160 m[2][1] -= tmp[0] * tsrc[12] + tmp[7] * tsrc[14] + tmp[8] * tsrc[15];
161 m[2][2] = tmp[2] * tsrc[12] + tmp[7] * tsrc[13] + tmp[10] * tsrc[15];
162 m[2][2] -= tmp[3] * tsrc[12] + tmp[6] * tsrc[13] + tmp[11] * tsrc[15];
163 m[2][3] = tmp[5] * tsrc[12] + tmp[8] * tsrc[13] + tmp[11] * tsrc[14];
164 m[2][3] -= tmp[4] * tsrc[12] + tmp[9] * tsrc[13] + tmp[10] * tsrc[14];
165 m[3][0] = tmp[2] * tsrc[10] + tmp[5] * tsrc[11] + tmp[1] * tsrc[9];
166 m[3][0] -= tmp[4] * tsrc[11] + tmp[0] * tsrc[9] + tmp[3] * tsrc[10];
167 m[3][1] = tmp[8] * tsrc[11] + tmp[0] * tsrc[8] + tmp[7] * tsrc[10];
168 m[3][1] -= tmp[6] * tsrc[10] + tmp[9] * tsrc[11] + tmp[1] * tsrc[8];
169 m[3][2] = tmp[6] * tsrc[9] + tmp[11] * tsrc[11] + tmp[3] * tsrc[8];
170 m[3][2] -= tmp[10] * tsrc[11] + tmp[2] * tsrc[8] + tmp[7] * tsrc[9];
171 m[3][3] = tmp[10] * tsrc[10] + tmp[4] * tsrc[8] + tmp[9] * tsrc[9];
172 m[3][3] -= tmp[8] * tsrc[9] + tmp[11] * tsrc[0] + tmp[5] * tsrc[8];
173
174 /* calculate determinant */
175 det = tsrc[0] * m[0][0] + tsrc[1] * m[0][1] + tsrc[2] * m[0][2] +
176 tsrc[3] * m[0][3];
177
178 /* calculate matrix inverse */
179 det = static_cast<T>(1.0) / det;
180
181 for (j = 0; j < 4; j++) {
182 for (i = 0; i < 4; i++) {
183 m[j][i] *= det;
184 }
185 }
186 }
187
Transpose(T m[4][4])188 static void Transpose(T m[4][4]) {
189 T t[4][4];
190
191 // Transpose
192 for (int j = 0; j < 4; j++) {
193 for (int i = 0; i < 4; i++) {
194 t[j][i] = m[i][j];
195 }
196 }
197
198 // Copy
199 for (int j = 0; j < 4; j++) {
200 for (int i = 0; i < 4; i++) {
201 m[j][i] = t[j][i];
202 }
203 }
204 }
205
Mult(T dst[4][4],const T m0[4][4],const T m1[4][4])206 static void Mult(T dst[4][4], const T m0[4][4], const T m1[4][4]) {
207 for (int i = 0; i < 4; ++i) {
208 for (int j = 0; j < 4; ++j) {
209 dst[i][j] = 0;
210 for (int k = 0; k < 4; ++k) {
211 dst[i][j] += m0[k][j] * m1[i][k];
212 }
213 }
214 }
215 }
216
MultV(T dst[3],const T m[4][4],const T v[3])217 static void MultV(T dst[3], const T m[4][4], const T v[3]) {
218 T tmp[3];
219 tmp[0] = m[0][0] * v[0] + m[1][0] * v[1] + m[2][0] * v[2] + m[3][0];
220 tmp[1] = m[0][1] * v[0] + m[1][1] * v[1] + m[2][1] * v[2] + m[3][1];
221 tmp[2] = m[0][2] * v[0] + m[1][2] * v[1] + m[2][2] * v[2] + m[3][2];
222 dst[0] = tmp[0];
223 dst[1] = tmp[1];
224 dst[2] = tmp[2];
225 }
226
MultV(nanort::real3<T> & dst,const T m[4][4],const T v[3])227 static void MultV(nanort::real3<T> &dst, const T m[4][4], const T v[3]) {
228 T tmp[3];
229 tmp[0] = m[0][0] * v[0] + m[1][0] * v[1] + m[2][0] * v[2] + m[3][0];
230 tmp[1] = m[0][1] * v[0] + m[1][1] * v[1] + m[2][1] * v[2] + m[3][1];
231 tmp[2] = m[0][2] * v[0] + m[1][2] * v[1] + m[2][2] * v[2] + m[3][2];
232 dst[0] = tmp[0];
233 dst[1] = tmp[1];
234 dst[2] = tmp[2];
235 }
236 };
237
238 // typedef Matrix<float> Matrixf;
239 // typedef Matrix<double> Matrixd;
240
241 template <typename T>
XformBoundingBox(T xbmin[3],T xbmax[3],T bmin[3],T bmax[3],T m[4][4])242 static void XformBoundingBox(T xbmin[3], // out
243 T xbmax[3], // out
244 T bmin[3], T bmax[3], T m[4][4]) {
245 // create bounding vertex from (bmin, bmax)
246 T b[8][3];
247
248 b[0][0] = bmin[0];
249 b[0][1] = bmin[1];
250 b[0][2] = bmin[2];
251 b[1][0] = bmax[0];
252 b[1][1] = bmin[1];
253 b[1][2] = bmin[2];
254 b[2][0] = bmin[0];
255 b[2][1] = bmax[1];
256 b[2][2] = bmin[2];
257 b[3][0] = bmax[0];
258 b[3][1] = bmax[1];
259 b[3][2] = bmin[2];
260
261 b[4][0] = bmin[0];
262 b[4][1] = bmin[1];
263 b[4][2] = bmax[2];
264 b[5][0] = bmax[0];
265 b[5][1] = bmin[1];
266 b[5][2] = bmax[2];
267 b[6][0] = bmin[0];
268 b[6][1] = bmax[1];
269 b[6][2] = bmax[2];
270 b[7][0] = bmax[0];
271 b[7][1] = bmax[1];
272 b[7][2] = bmax[2];
273
274 T xb[8][3];
275 for (int i = 0; i < 8; i++) {
276 Matrix<T>::MultV(xb[i], m, b[i]);
277 }
278
279 xbmin[0] = xb[0][0];
280 xbmin[1] = xb[0][1];
281 xbmin[2] = xb[0][2];
282 xbmax[0] = xb[0][0];
283 xbmax[1] = xb[0][1];
284 xbmax[2] = xb[0][2];
285
286 for (int i = 1; i < 8; i++) {
287 xbmin[0] = std::min(xb[i][0], xbmin[0]);
288 xbmin[1] = std::min(xb[i][1], xbmin[1]);
289 xbmin[2] = std::min(xb[i][2], xbmin[2]);
290
291 xbmax[0] = std::max(xb[i][0], xbmax[0]);
292 xbmax[1] = std::max(xb[i][1], xbmax[1]);
293 xbmax[2] = std::max(xb[i][2], xbmax[2]);
294 }
295 }
296
297 ///
298 /// Intersection info struct for Node intersector
299 ///
300 /// @tparam T Precision(usually `float` or `double`)
301 ///
302 template <typename T>
303 struct Intersection {
304 // required fields.
305 T t; // hit distance
306 unsigned int prim_id; // primitive ID of the hit
307 T u;
308 T v;
309
310 unsigned int node_id; // node ID of the hit.
311 nanort::real3<T> P; // intersection point
312 nanort::real3<T> Ns; // shading normal
313 nanort::real3<T> Ng; // geometric normal
314 };
315
316 ///
317 /// Renderable node
318 ///
319 /// @tparam T Type of xform and bounding box(usually `float` or `double`).
320 /// @tparam M Mesh class
321 ///
322 template <typename T, class M>
323 class Node {
324 public:
325 typedef Node<T, M> type;
326
Node(const M * mesh)327 explicit Node(const M *mesh) : mesh_(mesh) {
328 xbmin_[0] = xbmin_[1] = xbmin_[2] = std::numeric_limits<T>::max();
329 xbmax_[0] = xbmax_[1] = xbmax_[2] = -std::numeric_limits<T>::max();
330
331 lbmin_[0] = lbmin_[1] = lbmin_[2] = std::numeric_limits<T>::max();
332 lbmax_[0] = lbmax_[1] = lbmax_[2] = -std::numeric_limits<T>::max();
333
334 Matrix<T>::Identity(local_xform_);
335 Matrix<T>::Identity(xform_);
336 Matrix<T>::Identity(inv_xform_);
337 Matrix<T>::Identity(inv_xform33_);
338 inv_xform33_[3][3] = static_cast<T>(0.0);
339 Matrix<T>::Identity(inv_transpose_xform33_);
340 inv_transpose_xform33_[3][3] = static_cast<T>(0.0);
341 }
342
~Node()343 ~Node() {}
344
Copy(const type & rhs)345 void Copy(const type &rhs) {
346 Matrix<T>::Copy(local_xform_, rhs.local_xform_);
347 Matrix<T>::Copy(xform_, rhs.xform_);
348 Matrix<T>::Copy(inv_xform_, rhs.inv_xform_);
349 Matrix<T>::Copy(inv_xform33_, rhs.inv_xform33_);
350 Matrix<T>::Copy(inv_transpose_xform33_, rhs.inv_transpose_xform33_);
351
352 lbmin_[0] = rhs.lbmin_[0];
353 lbmin_[1] = rhs.lbmin_[1];
354 lbmin_[2] = rhs.lbmin_[2];
355
356 lbmax_[0] = rhs.lbmax_[0];
357 lbmax_[1] = rhs.lbmax_[1];
358 lbmax_[2] = rhs.lbmax_[2];
359
360 xbmin_[0] = rhs.xbmin_[0];
361 xbmin_[1] = rhs.xbmin_[1];
362 xbmin_[2] = rhs.xbmin_[2];
363
364 xbmax_[0] = rhs.xbmax_[0];
365 xbmax_[1] = rhs.xbmax_[1];
366 xbmax_[2] = rhs.xbmax_[2];
367
368 mesh_ = rhs.mesh_;
369 name_ = rhs.name_;
370
371 children_ = rhs.children_;
372 }
373
Node(const type & rhs)374 Node(const type &rhs) { Copy(rhs); }
375
376 const type &operator=(const type &rhs) {
377 Copy(rhs);
378 return (*this);
379 }
380
SetName(const std::string & name)381 void SetName(const std::string &name) { name_ = name; }
382
GetName()383 const std::string &GetName() const { return name_; }
384
385 ///
386 /// Add child node.
387 ///
AddChild(const type & child)388 void AddChild(const type &child) { children_.push_back(child); }
389
390 ///
391 /// Get chidren
392 ///
GetChildren()393 const std::vector<type> &GetChildren() const { return children_; }
394
GetChildren()395 std::vector<type> &GetChildren() { return children_; }
396
397 ///
398 /// Update internal state.
399 ///
Update(const T parent_xform[4][4])400 void Update(const T parent_xform[4][4]) {
401 if (!accel_.IsValid() && mesh_ && (mesh_->vertices.size() > 3) &&
402 (mesh_->faces.size() >= 3)) {
403 // Assume mesh is composed of triangle faces only.
404 nanort::TriangleMesh<float> triangle_mesh(
405 mesh_->vertices.data(), mesh_->faces.data(), mesh_->stride);
406 nanort::TriangleSAHPred<float> triangle_pred(
407 mesh_->vertices.data(), mesh_->faces.data(), mesh_->stride);
408
409 bool ret =
410 accel_.Build(static_cast<unsigned int>(mesh_->faces.size()) / 3,
411 triangle_mesh, triangle_pred);
412
413 // Update local bbox.
414 if (ret) {
415 accel_.BoundingBox(lbmin_, lbmax_);
416 }
417 }
418
419 // xform = parent_xform x local_xform
420 Matrix<T>::Mult(xform_, parent_xform, local_xform_);
421
422 // Compute the bounding box in world coordinate.
423 XformBoundingBox(xbmin_, xbmax_, lbmin_, lbmax_, xform_);
424
425 // Inverse(xform)
426 Matrix<T>::Copy(inv_xform_, xform_);
427 Matrix<T>::Inverse(inv_xform_);
428
429 // Clear translation, then inverse(xform)
430 Matrix<T>::Copy(inv_xform33_, xform_);
431 inv_xform33_[3][0] = static_cast<T>(0.0);
432 inv_xform33_[3][1] = static_cast<T>(0.0);
433 inv_xform33_[3][2] = static_cast<T>(0.0);
434 Matrix<T>::Inverse(inv_xform33_);
435
436 // Inverse transpose of xform33
437 Matrix<T>::Copy(inv_transpose_xform33_, inv_xform33_);
438 Matrix<T>::Transpose(inv_transpose_xform33_);
439
440 // Update children nodes
441 for (size_t i = 0; i < children_.size(); i++) {
442 children_[i].Update(xform_);
443 }
444 }
445
446 ///
447 /// Set local transformation.
448 ///
SetLocalXform(const T xform[4][4])449 void SetLocalXform(const T xform[4][4]) {
450 memcpy(local_xform_, xform, sizeof(float) * 16);
451 }
452
GetLocalXformPtr()453 const T *GetLocalXformPtr() const { return &local_xform_[0][0]; }
454
GetXformPtr()455 const T *GetXformPtr() const { return &xform_[0][0]; }
456
GetMesh()457 const M *GetMesh() const { return mesh_; }
458
GetAccel()459 const nanort::BVHAccel<T> &GetAccel() const { return accel_; }
460
GetWorldBoundingBox(T bmin[3],T bmax[3])461 inline void GetWorldBoundingBox(T bmin[3], T bmax[3]) const {
462 bmin[0] = xbmin_[0];
463 bmin[1] = xbmin_[1];
464 bmin[2] = xbmin_[2];
465
466 bmax[0] = xbmax_[0];
467 bmax[1] = xbmax_[1];
468 bmax[2] = xbmax_[2];
469 }
470
GetLocalBoundingBox(T bmin[3],T bmax[3])471 inline void GetLocalBoundingBox(T bmin[3], T bmax[3]) const {
472 bmin[0] = lbmin_[0];
473 bmin[1] = lbmin_[1];
474 bmin[2] = lbmin_[2];
475
476 bmax[0] = lbmax_[0];
477 bmax[1] = lbmax_[1];
478 bmax[2] = lbmax_[2];
479 }
480
481 T local_xform_[4][4]; // Node's local transformation matrix.
482 T xform_[4][4]; // Parent xform x local_xform.
483 T inv_xform_[4][4]; // inverse(xform). world -> local
484 T inv_xform33_[4][4]; // inverse(xform0 with upper-left 3x3 elemets only(for
485 // transforming direction vector)
486 T inv_transpose_xform33_[4][4]; // inverse(transpose(xform)) with upper-left
487 // 3x3 elements only(for transforming normal
488 // vector)
489
490 private:
491 // bounding box(local space)
492 T lbmin_[3];
493 T lbmax_[3];
494
495 // bounding box after xform(world space)
496 T xbmin_[3];
497 T xbmax_[3];
498
499 nanort::BVHAccel<T> accel_;
500
501 std::string name_;
502
503 const M *mesh_{nullptr};
504
505 std::vector<type> children_;
506 };
507
508 // -------------------------------------------------
509
510 // Predefined SAH predicator for cube.
511 template <typename T, class M>
512 class NodeBBoxPred {
513 public:
NodeBBoxPred(const std::vector<Node<T,M>> * nodes)514 NodeBBoxPred(const std::vector<Node<T, M> > *nodes)
515 : axis_(0), pos_(0.0f), nodes_(nodes) {}
516
Set(int axis,float pos)517 void Set(int axis, float pos) const {
518 axis_ = axis;
519 pos_ = pos;
520 }
521
operator()522 bool operator()(unsigned int i) const {
523 int axis = axis_;
524 float pos = pos_;
525
526 T bmin[3], bmax[3];
527
528 (*nodes_)[i].GetWorldBoundingBox(bmin, bmax);
529
530 T center = bmax[axis] - bmin[axis];
531
532 return (center < pos);
533 }
534
535 private:
536 mutable int axis_;
537 mutable float pos_;
538 const std::vector<Node<T, M> > *nodes_;
539 };
540
541 template <typename T, class M>
542 class NodeBBoxGeometry {
543 public:
NodeBBoxGeometry(const std::vector<Node<T,M>> * nodes)544 NodeBBoxGeometry(const std::vector<Node<T, M> > *nodes) : nodes_(nodes) {}
545
546 /// Compute bounding box for `prim_index`th cube.
547 /// This function is called for each primitive in BVH build.
BoundingBox(nanort::real3<T> * bmin,nanort::real3<T> * bmax,unsigned int prim_index)548 void BoundingBox(nanort::real3<T> *bmin, nanort::real3<T> *bmax,
549 unsigned int prim_index) const {
550 T a[3], b[3];
551 (*nodes_)[prim_index].GetWorldBoundingBox(a, b);
552 (*bmin)[0] = a[0];
553 (*bmin)[1] = a[1];
554 (*bmin)[2] = a[2];
555 (*bmax)[0] = b[0];
556 (*bmax)[1] = b[1];
557 (*bmax)[2] = b[2];
558 }
559
560 const std::vector<Node<T, M> > *nodes_;
561 mutable nanort::real3<T> ray_org_;
562 mutable nanort::real3<T> ray_dir_;
563 mutable nanort::BVHTraceOptions trace_options_;
564 int _pad_;
565 };
566
567 class NodeBBoxIntersection {
568 public:
NodeBBoxIntersection()569 NodeBBoxIntersection() {}
570
571 float normal[3];
572
573 // Required member variables.
574 float t;
575 unsigned int prim_id;
576 };
577
578 template <typename T, class M>
579 class NodeBBoxIntersector {
580 public:
NodeBBoxIntersector(const std::vector<Node<T,M>> * nodes)581 NodeBBoxIntersector(const std::vector<Node<T, M> > *nodes) : nodes_(nodes) {}
582
Intersect(float * out_t_min,float * out_t_max,unsigned int prim_index)583 bool Intersect(float *out_t_min, float *out_t_max,
584 unsigned int prim_index) const {
585 T bmin[3], bmax[3];
586
587 (*nodes_)[prim_index].GetWorldBoundingBox(bmin, bmax);
588
589 float tmin, tmax;
590
591 const float min_x = ray_dir_sign_[0] ? bmax[0] : bmin[0];
592 const float min_y = ray_dir_sign_[1] ? bmax[1] : bmin[1];
593 const float min_z = ray_dir_sign_[2] ? bmax[2] : bmin[2];
594 const float max_x = ray_dir_sign_[0] ? bmin[0] : bmax[0];
595 const float max_y = ray_dir_sign_[1] ? bmin[1] : bmax[1];
596 const float max_z = ray_dir_sign_[2] ? bmin[2] : bmax[2];
597
598 // X
599 const float tmin_x = (min_x - ray_org_[0]) * ray_inv_dir_[0];
600 const float tmax_x = (max_x - ray_org_[0]) * ray_inv_dir_[0];
601
602 // Y
603 const float tmin_y = (min_y - ray_org_[1]) * ray_inv_dir_[1];
604 const float tmax_y = (max_y - ray_org_[1]) * ray_inv_dir_[1];
605
606 // Z
607 const float tmin_z = (min_z - ray_org_[2]) * ray_inv_dir_[2];
608 const float tmax_z = (max_z - ray_org_[2]) * ray_inv_dir_[2];
609
610 tmin = nanort::safemax(tmin_z, nanort::safemax(tmin_y, tmin_x));
611 tmax = nanort::safemin(tmax_z, nanort::safemin(tmax_y, tmax_x));
612
613 if (tmin <= tmax) {
614 (*out_t_min) = tmin;
615 (*out_t_max) = tmax;
616 return true;
617 }
618
619 return false;
620 }
621
622 /// Prepare BVH traversal(e.g. compute inverse ray direction)
623 /// This function is called only once in BVH traversal.
PrepareTraversal(const nanort::Ray<float> & ray)624 void PrepareTraversal(const nanort::Ray<float> &ray) const {
625 ray_org_[0] = ray.org[0];
626 ray_org_[1] = ray.org[1];
627 ray_org_[2] = ray.org[2];
628
629 ray_dir_[0] = ray.dir[0];
630 ray_dir_[1] = ray.dir[1];
631 ray_dir_[2] = ray.dir[2];
632
633 // FIXME(syoyo): Consider zero div case.
634 ray_inv_dir_[0] = static_cast<T>(1.0) / ray.dir[0];
635 ray_inv_dir_[1] = static_cast<T>(1.0) / ray.dir[1];
636 ray_inv_dir_[2] = static_cast<T>(1.0) / ray.dir[2];
637
638 ray_dir_sign_[0] = ray.dir[0] < static_cast<T>(0.0) ? 1 : 0;
639 ray_dir_sign_[1] = ray.dir[1] < static_cast<T>(0.0) ? 1 : 0;
640 ray_dir_sign_[2] = ray.dir[2] < static_cast<T>(0.0) ? 1 : 0;
641 }
642
643 const std::vector<Node<T, M> > *nodes_;
644 mutable nanort::real3<T> ray_org_;
645 mutable nanort::real3<T> ray_dir_;
646 mutable nanort::real3<T> ray_inv_dir_;
647 mutable int ray_dir_sign_[3];
648 };
649
650 template <typename T, class M>
651 class Scene {
652 public:
Scene()653 Scene() {
654 bmin_[0] = bmin_[1] = bmin_[2] = std::numeric_limits<T>::max();
655 bmax_[0] = bmax_[1] = bmax_[2] = -std::numeric_limits<T>::max();
656 }
657
~Scene()658 ~Scene() {}
659
660 ///
661 /// Add intersectable node to the scene.
662 ///
AddNode(const Node<T,M> & node)663 bool AddNode(const Node<T, M> &node) {
664 nodes_.push_back(node);
665 return true;
666 }
667
GetNodes()668 const std::vector<Node<T, M> > &GetNodes() const { return nodes_; }
669
FindNode(const std::string & name,Node<T,M> ** found_node)670 bool FindNode(const std::string &name, Node<T, M> **found_node) {
671 if (!found_node) {
672 return false;
673 }
674
675 if (name.empty()) {
676 return false;
677 }
678
679 // Simple exhaustive search.
680 for (size_t i = 0; i < nodes_.size(); i++) {
681 if (FindNodeRecursive(name, &(nodes_[i]), found_node)) {
682 return true;
683 }
684 }
685
686 return false;
687 }
688
689 ///
690 /// Commit the scene. Must be called before tracing rays into the scene.
691 ///
Commit()692 bool Commit() {
693 // the scene should contains something
694 if (nodes_.size() == 0) {
695 std::cerr << "You are attempting to commit an empty scene!\n";
696 return false;
697 }
698
699 // Update nodes.
700 for (size_t i = 0; i < nodes_.size(); i++) {
701 T ident[4][4];
702 Matrix<T>::Identity(ident);
703
704 nodes_[i].Update(ident);
705 }
706
707 // Build toplevel BVH.
708 NodeBBoxGeometry<T, M> geom(&nodes_);
709 NodeBBoxPred<T, M> pred(&nodes_);
710
711 // FIXME(LTE): Limit one leaf contains one node bbox primitive. This would
712 // work, but would be inefficient.
713 // e.g. will miss some node when constructed BVH depth is larger than the
714 // value of BVHBuildOptions.
715 // Implement more better and efficient BVH build and traverse for Toplevel
716 // BVH.
717 nanort::BVHBuildOptions<T> build_options;
718 build_options.min_leaf_primitives = 1;
719
720 bool ret = toplevel_accel_.Build(static_cast<unsigned int>(nodes_.size()),
721 geom, pred, build_options);
722
723 nanort::BVHBuildStatistics stats = toplevel_accel_.GetStatistics();
724 (void)stats;
725
726 // toplevel_accel_.Debug();
727
728 if (ret) {
729 toplevel_accel_.BoundingBox(bmin_, bmax_);
730 } else {
731 // Set invalid bbox value.
732 bmin_[0] = std::numeric_limits<T>::max();
733 bmin_[1] = std::numeric_limits<T>::max();
734 bmin_[2] = std::numeric_limits<T>::max();
735
736 bmax_[0] = -std::numeric_limits<T>::max();
737 bmax_[1] = -std::numeric_limits<T>::max();
738 bmax_[2] = -std::numeric_limits<T>::max();
739 }
740
741 return ret;
742 }
743
744 ///
745 /// Get the scene bounding box.
746 ///
GetBoundingBox(T bmin[3],T bmax[3])747 void GetBoundingBox(T bmin[3], T bmax[3]) const {
748 bmin[0] = bmin_[0];
749 bmin[1] = bmin_[1];
750 bmin[2] = bmin_[2];
751
752 bmax[0] = bmax_[0];
753 bmax[1] = bmax_[1];
754 bmax[2] = bmax_[2];
755 }
756
757 ///
758 /// Trace the ray into the scene.
759 /// First find the intersection of nodes' bounding box using toplevel BVH.
760 /// Then, trace into the hit node to find the intersection of the primitive.
761 ///
762 /// @tparam H Hit intersection info class
763 /// @tparam I Intersector class
764 ///
765 template <class H, class I>
766 bool Traverse(nanort::Ray<T> &ray, H *isect,
767 const bool cull_back_face = false) const {
768 if (!toplevel_accel_.IsValid()) {
769 return false;
770 }
771
772 const int kMaxIntersections = 64;
773
774 bool has_hit = false;
775
776 NodeBBoxIntersector<T, M> isector(&nodes_);
777 nanort::StackVector<nanort::NodeHit<T>, 128> node_hits;
778 bool may_hit = toplevel_accel_.ListNodeIntersections(ray, kMaxIntersections,
779 isector, &node_hits);
780
781 if (may_hit) {
782 T t_max = std::numeric_limits<T>::max();
783 T t_nearest = t_max;
784
785 nanort::BVHTraceOptions trace_options;
786 trace_options.cull_back_face = cull_back_face;
787
788 // Find actual intersection point.
789 for (size_t i = 0; i < node_hits->size(); i++) {
790 // Early cull test.
791 if (t_nearest < node_hits[i].t_min) {
792 // printf("near: %f, t_min: %f, t_max: %f\n", t_nearest,
793 // node_hits[i].t_min, node_hits[i].t_max);
794 continue;
795 }
796
797 assert(node_hits[i].node_id < nodes_.size());
798 const Node<T, M> &node = nodes_[node_hits[i].node_id];
799
800 // Transform ray into node's local space
801 // TODO(LTE): Set ray tmin and tmax
802 nanort::Ray<T> local_ray;
803 Matrix<T>::MultV(local_ray.org, node.inv_xform_, ray.org);
804 Matrix<T>::MultV(local_ray.dir, node.inv_xform33_, ray.dir);
805
806 // TODO(LTE): Provide custom intersector
807 //nanort::TriangleIntersector<T, H> triangle_intersector(
808 // node.GetMesh()->vertices.data(), node.GetMesh()->faces.data(),
809 // node.GetMesh()->stride);
810 I intersector(node.GetMesh());
811 H local_isect;
812
813 bool hit = node.GetAccel().Traverse(local_ray, intersector,
814 &local_isect);
815
816 if (hit) {
817 // Calulcate hit distance in world coordiante.
818 T local_P[3];
819 local_P[0] = local_ray.org[0] + local_isect.t * local_ray.dir[0];
820 local_P[1] = local_ray.org[1] + local_isect.t * local_ray.dir[1];
821 local_P[2] = local_ray.org[2] + local_isect.t * local_ray.dir[2];
822
823 T world_P[3];
824 Matrix<T>::MultV(world_P, node.xform_, local_P);
825
826 nanort::real3<T> po;
827 po[0] = world_P[0] - ray.org[0];
828 po[1] = world_P[1] - ray.org[1];
829 po[2] = world_P[2] - ray.org[2];
830
831 float t_world = vlength(po);
832 // printf("tworld %f, tnear %f\n", t_world, t_nearest);
833
834 if (t_world < t_nearest) {
835 t_nearest = t_world;
836 has_hit = true;
837 //(*isect) = local_isect;
838 isect->node_id = node_hits[i].node_id;
839 isect->prim_id = local_isect.prim_id;
840 isect->u = local_isect.u;
841 isect->v = local_isect.v;
842
843 // TODO(LTE): Implement
844 T Ng[3], Ns[3]; // geometric normal, shading normal.
845
846 node.GetMesh()->GetNormal(Ng, Ns, isect->prim_id, isect->u,
847 isect->v);
848
849 // Convert position and normal into world coordinate.
850 isect->t = t_world;
851 Matrix<T>::MultV(isect->P, node.xform_, local_P);
852 Matrix<T>::MultV(isect->Ng, node.inv_transpose_xform33_, Ng);
853 Matrix<T>::MultV(isect->Ns, node.inv_transpose_xform33_, Ns);
854 }
855 }
856 }
857 }
858
859 return has_hit;
860 }
861
862 private:
863 ///
864 /// Find a node by name.
865 ///
FindNodeRecursive(const std::string & name,Node<T,M> * root,Node<T,M> ** found_node)866 bool FindNodeRecursive(const std::string &name, Node<T, M> *root,
867 Node<T, M> **found_node) {
868 if (root->GetName().compare(name) == 0) {
869 (*found_node) = root;
870 return true;
871 }
872
873 // Simple exhaustive search.
874 for (size_t i = 0; i < root->GetChildren().size(); i++) {
875 if (FindNodeRecursive(name, &(root->GetChildren()[i]), found_node)) {
876 return true;
877 }
878 }
879
880 return false;
881 }
882
883 // Scene bounding box.
884 // Valid after calling `Commit()`.
885 T bmin_[3];
886 T bmax_[3];
887
888 // Toplevel BVH accel.
889 nanort::BVHAccel<T> toplevel_accel_;
890 std::vector<Node<T, M> > nodes_;
891 };
892
893 } // namespace nanosg
894
895 #endif // NANOSG_H_
896