KiCad PCB EDA Suite
bvh_pbrt.cpp
Go to the documentation of this file.
1/*
2 * This program source code file is part of KiCad, a free EDA CAD application.
3 *
4 * Copyright (C) 2015-2016 Mario Luzeiro <[email protected]>
5 * Copyright (C) 2015-2020 KiCad Developers, see AUTHORS.txt for contributors.
6 *
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, you may find one here:
19 * http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
20 * or you may search the http://www.gnu.org website for the version 2 license,
21 * or you may write to the Free Software Foundation, Inc.,
22 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
23 */
24
69#include "bvh_pbrt.h"
70#include "../../../3d_fastmath.h"
71#include <macros.h>
72
73#include <boost/range/algorithm/nth_element.hpp>
74#include <boost/range/algorithm/partition.hpp>
75#include <cstdlib>
76#include <vector>
77
78#include <stack>
79#include <wx/debug.h>
80
81#ifdef PRINT_STATISTICS_3D_VIEWER
82#include <stdio.h>
83#endif
84
85// BVHAccel Local Declarations
87{
89 {
91 bounds.Reset();
92 centroid = SFVEC3F( 0.0f );
93 }
94
95 BVHPrimitiveInfo( int aPrimitiveNumber, const BBOX_3D& aBounds ) :
96 primitiveNumber( aPrimitiveNumber ),
97 bounds( aBounds ),
98 centroid( .5f * aBounds.Min() + .5f * aBounds.Max() ) {}
99
103};
104
105
107{
108 // BVHBuildNode Public Methods
109 void InitLeaf( int first, int n, const BBOX_3D& b)
110 {
111 firstPrimOffset = first;
112 nPrimitives = n;
113 bounds = b;
114 children[0] = children[1] = nullptr;
115 }
116
117 void InitInterior( int axis, BVHBuildNode* c0, BVHBuildNode* c1 )
118 {
119 children[0] = c0;
120 children[1] = c1;
121 bounds.Set( c0->bounds );
122 bounds.Union( c1->bounds );
123 splitAxis = axis;
124 nPrimitives = 0;
125 }
126
130};
131
132
134{
136 uint32_t mortonCode;
137};
138
139
141{
144};
145
146
147// BVHAccel Utility Functions
148inline uint32_t LeftShift3( uint32_t x )
149{
150 wxASSERT( x <= (1 << 10) );
151
152 if( x == ( 1 << 10 ) )
153 --x;
154
155 x = ( x | ( x << 16 ) ) & 0b00000011000000000000000011111111;
156 // x = ---- --98 ---- ---- ---- ---- 7654 3210
157 x = ( x | ( x << 8 ) ) & 0b00000011000000001111000000001111;
158 // x = ---- --98 ---- ---- 7654 ---- ---- 3210
159 x = ( x | ( x << 4 ) ) & 0b00000011000011000011000011000011;
160 // x = ---- --98 ---- 76-- --54 ---- 32-- --10
161 x = ( x | ( x << 2 ) ) & 0b00001001001001001001001001001001;
162 // x = ---- 9--8 --7- -6-- 5--4 --3- -2-- 1--0
163
164 return x;
165}
166
167
168inline uint32_t EncodeMorton3( const SFVEC3F& v )
169{
170 wxASSERT( v.x >= 0 && v.x <= ( 1 << 10 ) );
171 wxASSERT( v.y >= 0 && v.y <= ( 1 << 10 ) );
172 wxASSERT( v.z >= 0 && v.z <= ( 1 << 10 ) );
173
174 return ( LeftShift3( v.z ) << 2 ) | ( LeftShift3( v.y ) << 1 ) | LeftShift3( v.x );
175}
176
177
178static void RadixSort( std::vector<MortonPrimitive>* v )
179{
180 std::vector<MortonPrimitive> tempVector( v->size() );
181
182 const int bitsPerPass = 6;
183 const int nBits = 30;
184
185 wxASSERT( ( nBits % bitsPerPass ) == 0 );
186
187 const int nPasses = nBits / bitsPerPass;
188
189 for( int pass = 0; pass < nPasses; ++pass )
190 {
191 // Perform one pass of radix sort, sorting _bitsPerPass_ bits
192 const int lowBit = pass * bitsPerPass;
193
194 // Set in and out vector pointers for radix sort pass
195 std::vector<MortonPrimitive>& in = ( pass & 1 ) ? tempVector : *v;
196 std::vector<MortonPrimitive>& out = ( pass & 1 ) ? *v : tempVector;
197
198 // Count number of zero bits in array for current radix sort bit
199 const int nBuckets = 1 << bitsPerPass;
200 int bucketCount[nBuckets] = {0};
201 const int bitMask = ( 1 << bitsPerPass ) - 1;
202
203 for( uint32_t i = 0; i < in.size(); ++i )
204 {
205 const MortonPrimitive &mp = in[i];
206 int bucket = ( mp.mortonCode >> lowBit ) & bitMask;
207
208 wxASSERT( ( bucket >= 0 ) && ( bucket < nBuckets ) );
209
210 ++bucketCount[bucket];
211 }
212
213 // Compute starting index in output array for each bucket
214 int startIndex[nBuckets];
215 startIndex[0] = 0;
216
217 for( int i = 1; i < nBuckets; ++i )
218 startIndex[i] = startIndex[i - 1] + bucketCount[i - 1];
219
220 // Store sorted values in output array
221 for( uint32_t i = 0; i < in.size(); ++i )
222 {
223 const MortonPrimitive& mp = in[i];
224 int bucket = (mp.mortonCode >> lowBit) & bitMask;
225 out[startIndex[bucket]++] = mp;
226 }
227 }
228
229 // Copy final result from _tempVector_, if needed
230 if( nPasses & 1 )
231 std::swap( *v, tempVector );
232}
233
234
235BVH_PBRT::BVH_PBRT( const CONTAINER_3D_BASE& aObjectContainer, int aMaxPrimsInNode,
236 SPLITMETHOD aSplitMethod ) :
237 m_maxPrimsInNode( std::min( 255, aMaxPrimsInNode ) ),
238 m_splitMethod( aSplitMethod )
239{
240 if( aObjectContainer.GetList().empty() )
241 {
242 m_nodes = nullptr;
243
244 return;
245 }
246
247 // Initialize the indexes of ray packet for partition traversal
248 for( unsigned int i = 0; i < RAYPACKET_RAYS_PER_PACKET; ++i )
249 {
250 m_I[i] = i;
251 }
252
253 // Convert the objects list to vector of objects
254 aObjectContainer.ConvertTo( m_primitives );
255
256 wxASSERT( aObjectContainer.GetList().size() == m_primitives.size() );
257
258 // Initialize _primitiveInfo_ array for primitives
259 std::vector<BVHPrimitiveInfo> primitiveInfo( m_primitives.size() );
260
261 for( size_t i = 0; i < m_primitives.size(); ++i )
262 {
263 wxASSERT( m_primitives[i]->GetBBox().IsInitialized() );
264
265 primitiveInfo[i] = BVHPrimitiveInfo( i, m_primitives[i]->GetBBox() );
266 }
267
268 // Build BVH tree for primitives using _primitiveInfo_
269 int totalNodes = 0;
270
271 CONST_VECTOR_OBJECT orderedPrims;
272 orderedPrims.clear();
273 orderedPrims.reserve( m_primitives.size() );
274
275 BVHBuildNode *root;
276
278 root = HLBVHBuild( primitiveInfo, &totalNodes, orderedPrims );
279 else
280 root = recursiveBuild( primitiveInfo, 0, m_primitives.size(), &totalNodes, orderedPrims );
281
282 wxASSERT( m_primitives.size() == orderedPrims.size() );
283
284 m_primitives.swap( orderedPrims );
285
286 // Compute representation of depth-first traversal of BVH tree
287 m_nodes = static_cast<LinearBVHNode*>( malloc( sizeof( LinearBVHNode ) * totalNodes ) );
288 m_nodesToFree.push_back( m_nodes );
289
290 for( int i = 0; i < totalNodes; ++i )
291 {
292 m_nodes[i].bounds.Reset();
294 m_nodes[i].nPrimitives = 0;
295 m_nodes[i].axis = 0;
296 }
297
298 uint32_t offset = 0;
299
300 flattenBVHTree( root, &offset );
301
302 wxASSERT( offset == (unsigned int)totalNodes );
303}
304
305
307{
308 if( !m_nodesToFree.empty() )
309 {
310 for( std::list<void *>::iterator ii = m_nodesToFree.begin(); ii != m_nodesToFree.end();
311 ++ii )
312 {
313 free( *ii );
314 *ii = nullptr;
315 }
316 }
317
318 m_nodesToFree.clear();
319}
320
321
323{
324 explicit ComparePoints( int d ) { dim = d; }
325
326 int dim;
327
328 bool operator()( const BVHPrimitiveInfo& a, const BVHPrimitiveInfo& b ) const
329 {
330 return a.centroid[dim] < b.centroid[dim];
331 }
332};
333
334
336{
337 explicit CompareToMid( int d, float m ) { dim = d; mid = m; }
338
339 int dim;
340 float mid;
341
342 bool operator()( const BVHPrimitiveInfo& a ) const
343 {
344 return a.centroid[dim] < mid;
345 }
346};
347
348
350{
351 CompareToBucket( int split, int num, int d, const BBOX_3D& b )
352 : centroidBounds( b )
353 {
355 nBuckets = num;
356 dim = d;
357 }
358
359 bool operator()( const BVHPrimitiveInfo& p ) const;
360
362
364};
365
366
368{
369 const float centroid = p.centroid[dim];
370
371 int b = nBuckets *
372 // Computes the offset (0.0 - 1.0) for one axis
373 ( ( centroid - centroidBounds.Min()[dim] ) /
375
376 if( b == nBuckets )
377 b = nBuckets - 1;
378
379 wxASSERT( ( b >= 0 ) && ( b < nBuckets ) );
380
381 return b <= splitBucket;
382}
383
384
386{
387 HLBVH_SAH_Evaluator( int split, int num, int d, const BBOX_3D& b )
388 : centroidBounds(b)
389 {
391 nBuckets = num; dim = d;
392 }
393
394 bool operator()( const BVHBuildNode* node ) const;
395
398};
399
400
402{
403 const float centroid = node->bounds.GetCenter( dim );
404
405 int b = nBuckets *
406 // Computes the offset (0.0 - 1.0) for one axis
407 ( ( centroid - centroidBounds.Min()[dim] ) /
409
410 if( b == nBuckets )
411 b = nBuckets - 1;
412
413 wxASSERT( b >= 0 && b < nBuckets );
414
415 return b <= minCostSplitBucket;
416}
417
418
420{
421 int count;
423};
424
425
426BVHBuildNode *BVH_PBRT::recursiveBuild ( std::vector<BVHPrimitiveInfo>& primitiveInfo,
427 int start, int end, int* totalNodes,
428 CONST_VECTOR_OBJECT& orderedPrims )
429{
430 wxASSERT( totalNodes != nullptr );
431 wxASSERT( start >= 0 );
432 wxASSERT( end >= 0 );
433 wxASSERT( start != end );
434 wxASSERT( start < end );
435 wxASSERT( start <= (int)primitiveInfo.size() );
436 wxASSERT( end <= (int)primitiveInfo.size() );
437
438 (*totalNodes)++;
439
440 // !TODO: implement an memory Arena
441 BVHBuildNode *node = static_cast<BVHBuildNode *>( malloc( sizeof( BVHBuildNode ) ) );
442 m_nodesToFree.push_back( node );
443
444 node->bounds.Reset();
445 node->firstPrimOffset = 0;
446 node->nPrimitives = 0;
447 node->splitAxis = 0;
448 node->children[0] = nullptr;
449 node->children[1] = nullptr;
450
451 // Compute bounds of all primitives in BVH node
452 BBOX_3D bounds;
453 bounds.Reset();
454
455 for( int i = start; i < end; ++i )
456 bounds.Union( primitiveInfo[i].bounds );
457
458 int nPrimitives = end - start;
459
460 if( nPrimitives == 1 )
461 {
462 // Create leaf _BVHBuildNode_
463 int firstPrimOffset = orderedPrims.size();
464
465 for( int i = start; i < end; ++i )
466 {
467 int primitiveNr = primitiveInfo[i].primitiveNumber;
468 wxASSERT( primitiveNr < (int)m_primitives.size() );
469 orderedPrims.push_back( m_primitives[ primitiveNr ] );
470 }
471
472 node->InitLeaf( firstPrimOffset, nPrimitives, bounds );
473 }
474 else
475 {
476 // Compute bound of primitive centroids, choose split dimension _dim_
477 BBOX_3D centroidBounds;
478 centroidBounds.Reset();
479
480 for( int i = start; i < end; ++i )
481 centroidBounds.Union( primitiveInfo[i].centroid );
482
483 const int dim = centroidBounds.MaxDimension();
484
485 // Partition primitives into two sets and build children
486 int mid = (start + end) / 2;
487
488 if( fabs( centroidBounds.Max()[dim] -
489 centroidBounds.Min()[dim] ) < (FLT_EPSILON + FLT_EPSILON) )
490 {
491 // Create leaf _BVHBuildNode_
492 const int firstPrimOffset = orderedPrims.size();
493
494 for( int i = start; i < end; ++i )
495 {
496 int primitiveNr = primitiveInfo[i].primitiveNumber;
497
498 wxASSERT( ( primitiveNr >= 0 ) && ( primitiveNr < (int) m_primitives.size() ) );
499
500 const OBJECT_3D* obj = static_cast<const OBJECT_3D*>( m_primitives[ primitiveNr ] );
501
502 wxASSERT( obj != nullptr );
503
504 orderedPrims.push_back( obj );
505 }
506
507 node->InitLeaf( firstPrimOffset, nPrimitives, bounds );
508 }
509 else
510 {
511 // Partition primitives based on _splitMethod_
512 switch( m_splitMethod )
513 {
515 {
516 // Partition primitives through node's midpoint
517 float pmid = centroidBounds.GetCenter( dim );
518
519 BVHPrimitiveInfo *midPtr = std::partition( &primitiveInfo[start],
520 &primitiveInfo[end - 1] + 1,
521 CompareToMid( dim, pmid ) );
522 mid = midPtr - &primitiveInfo[0];
523
524 wxASSERT( ( mid >= start ) && ( mid <= end ) );
525
526 if( ( mid != start ) && ( mid != end ) )
527 break;
528 }
529
530 // Intentionally fall through to SPLITMETHOD::EQUAL_COUNTS since prims
531 // with large overlapping bounding boxes may fail to partition
533
535 {
536 // Partition primitives into equally-sized subsets
537 mid = ( start + end ) / 2;
538
539 std::nth_element( &primitiveInfo[start], &primitiveInfo[mid],
540 &primitiveInfo[end - 1] + 1, ComparePoints( dim ) );
541
542 break;
543 }
544
545 case SPLITMETHOD::SAH:
546 default:
547 {
548 // Partition primitives using approximate SAH
549 if( nPrimitives <= 2 )
550 {
551 // Partition primitives into equally-sized subsets
552 mid = ( start + end ) / 2;
553
554 std::nth_element( &primitiveInfo[start], &primitiveInfo[mid],
555 &primitiveInfo[end - 1] + 1, ComparePoints( dim ) );
556 }
557 else
558 {
559 // Allocate _BucketInfo_ for SAH partition buckets
560 const int nBuckets = 12;
561
562 BucketInfo buckets[nBuckets];
563
564 for( int i = 0; i < nBuckets; ++i )
565 {
566 buckets[i].count = 0;
567 buckets[i].bounds.Reset();
568 }
569
570 // Initialize _BucketInfo_ for SAH partition buckets
571 for( int i = start; i < end; ++i )
572 {
573 int b = nBuckets *
574 centroidBounds.Offset( primitiveInfo[i].centroid )[dim];
575
576 if( b == nBuckets )
577 b = nBuckets - 1;
578
579 wxASSERT( b >= 0 && b < nBuckets );
580
581 buckets[b].count++;
582 buckets[b].bounds.Union( primitiveInfo[i].bounds );
583 }
584
585 // Compute costs for splitting after each bucket
586 float cost[nBuckets - 1];
587
588 for( int i = 0; i < ( nBuckets - 1 ); ++i )
589 {
590 BBOX_3D b0, b1;
591
592 b0.Reset();
593 b1.Reset();
594
595 int count0 = 0;
596 int count1 = 0;
597
598 for( int j = 0; j <= i; ++j )
599 {
600 if( buckets[j].count )
601 {
602 count0 += buckets[j].count;
603 b0.Union( buckets[j].bounds );
604 }
605 }
606
607 for( int j = i + 1; j < nBuckets; ++j )
608 {
609 if( buckets[j].count )
610 {
611 count1 += buckets[j].count;
612 b1.Union( buckets[j].bounds );
613 }
614 }
615
616 cost[i] = 1.0f + ( count0 * b0.SurfaceArea() +
617 count1 * b1.SurfaceArea() ) / bounds.SurfaceArea();
618 }
619
620 // Find bucket to split at that minimizes SAH metric
621 float minCost = cost[0];
622 int minCostSplitBucket = 0;
623
624 for( int i = 1; i < ( nBuckets - 1 ); ++i )
625 {
626 if( cost[i] < minCost )
627 {
628 minCost = cost[i];
629 minCostSplitBucket = i;
630 }
631 }
632
633 // Either create leaf or split primitives at selected SAH bucket
634 if( ( nPrimitives > m_maxPrimsInNode ) || ( minCost < (float) nPrimitives ) )
635 {
636 BVHPrimitiveInfo *pmid =
637 std::partition( &primitiveInfo[start], &primitiveInfo[end - 1] + 1,
638 CompareToBucket( minCostSplitBucket, nBuckets,
639 dim, centroidBounds ) );
640 mid = pmid - &primitiveInfo[0];
641
642 wxASSERT( ( mid >= start ) && ( mid <= end ) );
643 }
644 else
645 {
646 // Create leaf _BVHBuildNode_
647 const int firstPrimOffset = orderedPrims.size();
648
649 for( int i = start; i < end; ++i )
650 {
651 const int primitiveNr = primitiveInfo[i].primitiveNumber;
652
653 wxASSERT( primitiveNr < (int)m_primitives.size() );
654
655 orderedPrims.push_back( m_primitives[ primitiveNr ] );
656 }
657
658 node->InitLeaf( firstPrimOffset, nPrimitives, bounds );
659
660 return node;
661 }
662 }
663 break;
664 }
665 }
666
667 node->InitInterior( dim, recursiveBuild( primitiveInfo, start, mid, totalNodes,
668 orderedPrims ),
669 recursiveBuild( primitiveInfo, mid, end, totalNodes,
670 orderedPrims) );
671 }
672 }
673
674 return node;
675}
676
677
678BVHBuildNode *BVH_PBRT::HLBVHBuild( const std::vector<BVHPrimitiveInfo>& primitiveInfo,
679 int* totalNodes, CONST_VECTOR_OBJECT& orderedPrims )
680{
681 // Compute bounding box of all primitive centroids
682 BBOX_3D bounds;
683 bounds.Reset();
684
685 for( unsigned int i = 0; i < primitiveInfo.size(); ++i )
686 bounds.Union( primitiveInfo[i].centroid );
687
688 // Compute Morton indices of primitives
689 std::vector<MortonPrimitive> mortonPrims( primitiveInfo.size() );
690
691 for( int i = 0; i < (int)primitiveInfo.size(); ++i )
692 {
693 // Initialize _mortonPrims[i]_ for _i_th primitive
694 const int mortonBits = 10;
695 const int mortonScale = 1 << mortonBits;
696
697 wxASSERT( primitiveInfo[i].primitiveNumber < (int)primitiveInfo.size() );
698
699 mortonPrims[i].primitiveIndex = primitiveInfo[i].primitiveNumber;
700
701 const SFVEC3F centroidOffset = bounds.Offset( primitiveInfo[i].centroid );
702
703 wxASSERT( ( centroidOffset.x >= 0.0f ) && ( centroidOffset.x <= 1.0f ) );
704 wxASSERT( ( centroidOffset.y >= 0.0f ) && ( centroidOffset.y <= 1.0f ) );
705 wxASSERT( ( centroidOffset.z >= 0.0f ) && ( centroidOffset.z <= 1.0f ) );
706
707 mortonPrims[i].mortonCode = EncodeMorton3( centroidOffset * SFVEC3F( (float)mortonScale ) );
708 }
709
710 // Radix sort primitive Morton indices
711 RadixSort( &mortonPrims );
712
713 // Create LBVH treelets at bottom of BVH
714
715 // Find intervals of primitives for each treelet
716 std::vector<LBVHTreelet> treeletsToBuild;
717
718 for( int start = 0, end = 1; end <= (int)mortonPrims.size(); ++end )
719 {
720 const uint32_t mask = 0b00111111111111000000000000000000;
721
722 if( ( end == (int) mortonPrims.size() )
723 || ( ( mortonPrims[start].mortonCode & mask )
724 != ( mortonPrims[end].mortonCode & mask ) ) )
725 {
726 // Add entry to _treeletsToBuild_ for this treelet
727 const int numPrimitives = end - start;
728 const int maxBVHNodes = 2 * numPrimitives;
729
730 // !TODO: implement a memory arena
731 BVHBuildNode *nodes = static_cast<BVHBuildNode *>( malloc( maxBVHNodes *
732 sizeof( BVHBuildNode ) ) );
733
734 m_nodesToFree.push_back( nodes );
735
736 for( int i = 0; i < maxBVHNodes; ++i )
737 {
738 nodes[i].bounds.Reset();
739 nodes[i].firstPrimOffset = 0;
740 nodes[i].nPrimitives = 0;
741 nodes[i].splitAxis = 0;
742 nodes[i].children[0] = nullptr;
743 nodes[i].children[1] = nullptr;
744 }
745
746 LBVHTreelet tmpTreelet;
747
748 tmpTreelet.startIndex = start;
749 tmpTreelet.numPrimitives = numPrimitives;
750 tmpTreelet.buildNodes = nodes;
751
752 treeletsToBuild.push_back( tmpTreelet );
753
754 start = end;
755 }
756 }
757
758 // Create LBVHs for treelets in parallel
759 int atomicTotal = 0;
760 int orderedPrimsOffset = 0;
761
762 orderedPrims.resize( m_primitives.size() );
763
764 for( int index = 0; index < (int)treeletsToBuild.size(); ++index )
765 {
766 // Generate _index_th LBVH treelet
767 int nodesCreated = 0;
768 const int firstBit = 29 - 12;
769
770 LBVHTreelet &tr = treeletsToBuild[index];
771
772 wxASSERT( tr.startIndex < (int)mortonPrims.size() );
773
774 tr.buildNodes = emitLBVH( tr.buildNodes, primitiveInfo, &mortonPrims[tr.startIndex],
775 tr.numPrimitives, &nodesCreated, orderedPrims,
776 &orderedPrimsOffset, firstBit );
777
778 atomicTotal += nodesCreated;
779 }
780
781 *totalNodes = atomicTotal;
782
783 // Initialize _finishedTreelets_ with treelet root node pointers
784 std::vector<BVHBuildNode *> finishedTreelets;
785 finishedTreelets.reserve( treeletsToBuild.size() );
786
787 for( int index = 0; index < (int)treeletsToBuild.size(); ++index )
788 finishedTreelets.push_back( treeletsToBuild[index].buildNodes );
789
790 // Create and return SAH BVH from LBVH treelets
791 return buildUpperSAH( finishedTreelets, 0, finishedTreelets.size(), totalNodes );
792}
793
794
796 const std::vector<BVHPrimitiveInfo>& primitiveInfo,
797 MortonPrimitive* mortonPrims, int nPrimitives, int* totalNodes,
798 CONST_VECTOR_OBJECT& orderedPrims, int *orderedPrimsOffset,
799 int bit )
800{
801 wxASSERT( nPrimitives > 0 );
802 wxASSERT( totalNodes != nullptr );
803 wxASSERT( orderedPrimsOffset != nullptr );
804 wxASSERT( nPrimitives > 0 );
805 wxASSERT( mortonPrims != nullptr );
806
807 if( ( bit == -1 ) || ( nPrimitives < m_maxPrimsInNode ) )
808 {
809 // Create and return leaf node of LBVH treelet
810 (*totalNodes)++;
811
812 BVHBuildNode *node = buildNodes++;
813 BBOX_3D bounds;
814 bounds.Reset();
815
816 int firstPrimOffset = *orderedPrimsOffset;
817 *orderedPrimsOffset += nPrimitives;
818
819 wxASSERT( ( firstPrimOffset + ( nPrimitives - 1 ) ) < (int) orderedPrims.size() );
820
821 for( int i = 0; i < nPrimitives; ++i )
822 {
823 const int primitiveIndex = mortonPrims[i].primitiveIndex;
824
825 wxASSERT( primitiveIndex < (int)m_primitives.size() );
826
827 orderedPrims[firstPrimOffset + i] = m_primitives[primitiveIndex];
828 bounds.Union( primitiveInfo[primitiveIndex].bounds );
829 }
830
831 node->InitLeaf( firstPrimOffset, nPrimitives, bounds );
832
833 return node;
834 }
835 else
836 {
837 int mask = 1 << bit;
838
839 // Advance to next subtree level if there's no LBVH split for this bit
840 if( ( mortonPrims[0].mortonCode & mask ) ==
841 ( mortonPrims[nPrimitives - 1].mortonCode & mask ) )
842 return emitLBVH( buildNodes, primitiveInfo, mortonPrims, nPrimitives, totalNodes,
843 orderedPrims, orderedPrimsOffset, bit - 1 );
844
845 // Find LBVH split point for this dimension
846 int searchStart = 0;
847 int searchEnd = nPrimitives - 1;
848
849 while( searchStart + 1 != searchEnd )
850 {
851 wxASSERT(searchStart != searchEnd);
852
853 const int mid = ( searchStart + searchEnd ) / 2;
854
855 if( ( mortonPrims[searchStart].mortonCode & mask ) ==
856 ( mortonPrims[mid].mortonCode & mask ) )
857 searchStart = mid;
858 else
859 {
860 wxASSERT( ( mortonPrims[mid].mortonCode & mask ) ==
861 ( mortonPrims[searchEnd].mortonCode & mask ) );
862 searchEnd = mid;
863 }
864 }
865
866 const int splitOffset = searchEnd;
867
868 wxASSERT( splitOffset <= ( nPrimitives - 1 ) );
869 wxASSERT( ( mortonPrims[splitOffset - 1].mortonCode & mask )
870 != ( mortonPrims[splitOffset].mortonCode & mask ) );
871
872 // Create and return interior LBVH node
873 (*totalNodes)++;
874
875 BVHBuildNode *node = buildNodes++;
876 BVHBuildNode *lbvh[2];
877
878 lbvh[0] = emitLBVH( buildNodes, primitiveInfo, mortonPrims, splitOffset,
879 totalNodes, orderedPrims, orderedPrimsOffset, bit - 1 );
880
881 lbvh[1] = emitLBVH( buildNodes, primitiveInfo, &mortonPrims[splitOffset],
882 nPrimitives - splitOffset, totalNodes, orderedPrims,
883 orderedPrimsOffset, bit - 1 );
884
885 const int axis = bit % 3;
886
887 node->InitInterior( axis, lbvh[0], lbvh[1] );
888
889 return node;
890 }
891}
892
893
894BVHBuildNode *BVH_PBRT::buildUpperSAH( std::vector<BVHBuildNode*>& treeletRoots, int start,
895 int end, int* totalNodes )
896{
897 wxASSERT( totalNodes != nullptr );
898 wxASSERT( start < end );
899 wxASSERT( end <= (int)treeletRoots.size() );
900
901 int nNodes = end - start;
902
903 if( nNodes == 1 )
904 return treeletRoots[start];
905
906 (*totalNodes)++;
907
908 BVHBuildNode* node = static_cast<BVHBuildNode*>( malloc( sizeof( BVHBuildNode ) ) );
909
910 m_nodesToFree.push_back( node );
911
912 node->bounds.Reset();
913 node->firstPrimOffset = 0;
914 node->nPrimitives = 0;
915 node->splitAxis = 0;
916 node->children[0] = nullptr;
917 node->children[1] = nullptr;
918
919 // Compute bounds of all nodes under this HLBVH node
920 BBOX_3D bounds;
921 bounds.Reset();
922
923 for( int i = start; i < end; ++i )
924 bounds.Union( treeletRoots[i]->bounds );
925
926 // Compute bound of HLBVH node centroids, choose split dimension _dim_
927 BBOX_3D centroidBounds;
928 centroidBounds.Reset();
929
930 for( int i = start; i < end; ++i )
931 {
932 SFVEC3F centroid = ( treeletRoots[i]->bounds.Min() + treeletRoots[i]->bounds.Max() ) * 0.5f;
933
934 centroidBounds.Union( centroid );
935 }
936
937 const int dim = centroidBounds.MaxDimension();
938
939 // FIXME: if this hits, what do we need to do?
940 // Make sure the SAH split below does something... ?
941 wxASSERT( centroidBounds.Max()[dim] != centroidBounds.Min()[dim] );
942
943 // Allocate _BucketInfo_ for SAH partition buckets
944 const int nBuckets = 12;
945
946 BucketInfo buckets[nBuckets];
947
948 for( int i = 0; i < nBuckets; ++i )
949 {
950 buckets[i].count = 0;
951 buckets[i].bounds.Reset();
952 }
953
954 // Initialize _BucketInfo_ for HLBVH SAH partition buckets
955 for( int i = start; i < end; ++i )
956 {
957 const float centroid = ( treeletRoots[i]->bounds.Min()[dim] +
958 treeletRoots[i]->bounds.Max()[dim] ) * 0.5f;
959 int b = nBuckets * ( ( centroid - centroidBounds.Min()[dim] )
960 / ( centroidBounds.Max()[dim] - centroidBounds.Min()[dim] ) );
961
962 if( b == nBuckets )
963 b = nBuckets - 1;
964
965 wxASSERT( ( b >= 0 ) && ( b < nBuckets ) );
966
967 buckets[b].count++;
968 buckets[b].bounds.Union( treeletRoots[i]->bounds );
969 }
970
971 // Compute costs for splitting after each bucket
972 float cost[nBuckets - 1];
973
974 for( int i = 0; i < nBuckets - 1; ++i )
975 {
976 BBOX_3D b0, b1;
977 b0.Reset();
978 b1.Reset();
979
980 int count0 = 0, count1 = 0;
981
982 for( int j = 0; j <= i; ++j )
983 {
984 if( buckets[j].count )
985 {
986 count0 += buckets[j].count;
987 b0.Union( buckets[j].bounds );
988 }
989 }
990
991 for( int j = i + 1; j < nBuckets; ++j )
992 {
993 if( buckets[j].count )
994 {
995 count1 += buckets[j].count;
996 b1.Union( buckets[j].bounds );
997 }
998 }
999
1000 cost[i] = .125f + ( count0 * b0.SurfaceArea() + count1 * b1.SurfaceArea() ) /
1001 bounds.SurfaceArea();
1002 }
1003
1004 // Find bucket to split at that minimizes SAH metric
1005 float minCost = cost[0];
1006 int minCostSplitBucket = 0;
1007
1008 for( int i = 1; i < nBuckets - 1; ++i )
1009 {
1010 if( cost[i] < minCost )
1011 {
1012 minCost = cost[i];
1013 minCostSplitBucket = i;
1014 }
1015 }
1016
1017 // Split nodes and create interior HLBVH SAH node
1018 BVHBuildNode **pmid = std::partition( &treeletRoots[start], &treeletRoots[end - 1] + 1,
1019 HLBVH_SAH_Evaluator( minCostSplitBucket, nBuckets,
1020 dim, centroidBounds ) );
1021
1022 const int mid = pmid - &treeletRoots[0];
1023
1024 wxASSERT( ( mid > start ) && ( mid < end ) );
1025
1026 node->InitInterior( dim,
1027 buildUpperSAH( treeletRoots, start, mid, totalNodes ),
1028 buildUpperSAH( treeletRoots, mid, end, totalNodes ) );
1029
1030 return node;
1031}
1032
1033
1034int BVH_PBRT::flattenBVHTree( BVHBuildNode* node, uint32_t* offset )
1035{
1036 LinearBVHNode *linearNode = &m_nodes[*offset];
1037
1038 linearNode->bounds = node->bounds;
1039
1040 int myOffset = (*offset)++;
1041
1042 if( node->nPrimitives > 0 )
1043 {
1044 wxASSERT( ( !node->children[0] ) && ( !node->children[1] ) );
1045 wxASSERT( node->nPrimitives < 65536 );
1046
1047 linearNode->primitivesOffset = node->firstPrimOffset;
1048 linearNode->nPrimitives = node->nPrimitives;
1049 }
1050 else
1051 {
1052 // Create interior flattened BVH node
1053 linearNode->axis = node->splitAxis;
1054 linearNode->nPrimitives = 0;
1055 flattenBVHTree( node->children[0], offset );
1056 linearNode->secondChildOffset = flattenBVHTree( node->children[1], offset );
1057 }
1058
1059 return myOffset;
1060}
1061
1062
1063#define MAX_TODOS 64
1064
1065
1066bool BVH_PBRT::Intersect( const RAY& aRay, HITINFO& aHitInfo ) const
1067{
1068 if( !m_nodes )
1069 return false;
1070
1071 bool hit = false;
1072
1073 // Follow ray through BVH nodes to find primitive intersections
1074 int todoOffset = 0, nodeNum = 0;
1075 int todo[MAX_TODOS];
1076
1077 while( true )
1078 {
1079 const LinearBVHNode *node = &m_nodes[nodeNum];
1080
1081 wxASSERT( todoOffset < MAX_TODOS );
1082
1083 // Check ray against BVH node
1084 float hitBox = 0.0f;
1085
1086 const bool hitted = node->bounds.Intersect( aRay, &hitBox );
1087
1088 if( hitted && ( hitBox < aHitInfo.m_tHit ) )
1089 {
1090 if( node->nPrimitives > 0 )
1091 {
1092 // Intersect ray with primitives in leaf BVH node
1093 for( int i = 0; i < node->nPrimitives; ++i )
1094 {
1095 if( m_primitives[node->primitivesOffset + i]->Intersect( aRay, aHitInfo ) )
1096 {
1097 aHitInfo.m_acc_node_info = nodeNum;
1098 hit = true;
1099 }
1100 }
1101 }
1102 else
1103 {
1104 // Put far BVH node on _todo_ stack, advance to near node
1105 if( aRay.m_dirIsNeg[node->axis] )
1106 {
1107 todo[todoOffset++] = nodeNum + 1;
1108 nodeNum = node->secondChildOffset;
1109 }
1110 else
1111 {
1112 todo[todoOffset++] = node->secondChildOffset;
1113 nodeNum = nodeNum + 1;
1114 }
1115
1116 continue;
1117 }
1118 }
1119
1120 if( todoOffset == 0 )
1121 break;
1122
1123 nodeNum = todo[--todoOffset];
1124 }
1125
1126 return hit;
1127}
1128
1129
1131bool BVH_PBRT::Intersect( const RAY& aRay, HITINFO& aHitInfo, unsigned int aAccNodeInfo ) const
1132{
1133 if( !m_nodes )
1134 return false;
1135
1136 bool hit = false;
1137
1138 // Follow ray through BVH nodes to find primitive intersections
1139 int todoOffset = 0, nodeNum = aAccNodeInfo;
1140 int todo[MAX_TODOS];
1141
1142 while( true )
1143 {
1144 const LinearBVHNode *node = &m_nodes[nodeNum];
1145
1146 wxASSERT( todoOffset < MAX_TODOS );
1147
1148 // Check ray against BVH node
1149 float hitBox = 0.0f;
1150
1151 const bool hitted = node->bounds.Intersect( aRay, &hitBox );
1152
1153 if( hitted && ( hitBox < aHitInfo.m_tHit ) )
1154 {
1155 if( node->nPrimitives > 0 )
1156 {
1157 // Intersect ray with primitives in leaf BVH node
1158 for( int i = 0; i < node->nPrimitives; ++i )
1159 {
1160 if( m_primitives[node->primitivesOffset + i]->Intersect( aRay, aHitInfo ) )
1161 {
1162 //aHitInfo.m_acc_node_info = nodeNum;
1163 hit = true;
1164 }
1165 }
1166 }
1167 else
1168 {
1169 // Put far BVH node on _todo_ stack, advance to near node
1170 if( aRay.m_dirIsNeg[node->axis] )
1171 {
1172 todo[todoOffset++] = nodeNum + 1;
1173 nodeNum = node->secondChildOffset;
1174 }
1175 else
1176 {
1177 todo[todoOffset++] = node->secondChildOffset;
1178 nodeNum = nodeNum + 1;
1179 }
1180
1181 continue;
1182 }
1183 }
1184
1185 if( todoOffset == 0 )
1186 break;
1187
1188 nodeNum = todo[--todoOffset];
1189 }
1190
1191 return hit;
1192}
1193
1194
1195bool BVH_PBRT::IntersectP( const RAY& aRay, float aMaxDistance ) const
1196{
1197 if( !m_nodes )
1198 return false;
1199
1200 // Follow ray through BVH nodes to find primitive intersections
1201 int todoOffset = 0, nodeNum = 0;
1202 int todo[MAX_TODOS];
1203
1204 while( true )
1205 {
1206 const LinearBVHNode* node = &m_nodes[nodeNum];
1207
1208 wxASSERT( todoOffset < MAX_TODOS );
1209
1210 // Check ray against BVH node
1211 float hitBox = 0.0f;
1212
1213 const bool hitted = node->bounds.Intersect( aRay, &hitBox );
1214
1215 if( hitted && ( hitBox < aMaxDistance ) )
1216 {
1217 if( node->nPrimitives > 0 )
1218 {
1219 // Intersect ray with primitives in leaf BVH node
1220 for( int i = 0; i < node->nPrimitives; ++i )
1221 {
1222 const OBJECT_3D* obj = m_primitives[node->primitivesOffset + i];
1223
1224 if( obj->GetMaterial()->GetCastShadows()
1225 && obj->IntersectP( aRay, aMaxDistance ) )
1226 return true;
1227 }
1228 }
1229 else
1230 {
1231 // Put far BVH node on _todo_ stack, advance to near node
1232 if( aRay.m_dirIsNeg[node->axis] )
1233 {
1234 todo[todoOffset++] = nodeNum + 1;
1235 nodeNum = node->secondChildOffset;
1236 }
1237 else
1238 {
1239 todo[todoOffset++] = node->secondChildOffset;
1240 nodeNum = nodeNum + 1;
1241 }
1242
1243 continue;
1244 }
1245 }
1246
1247 if( todoOffset == 0 )
1248 break;
1249
1250 nodeNum = todo[--todoOffset];
1251 }
1252
1253 return false;
1254}
uint32_t EncodeMorton3(const SFVEC3F &v)
Definition: bvh_pbrt.cpp:168
uint32_t LeftShift3(uint32_t x)
Definition: bvh_pbrt.cpp:148
#define MAX_TODOS
Definition: bvh_pbrt.cpp:1063
static void RadixSort(std::vector< MortonPrimitive > *v)
Definition: bvh_pbrt.cpp:178
This BVH implementation is based on the source code implementation from the book "Physically Based Re...
SPLITMETHOD
Definition: bvh_pbrt.h:102
std::list< void * > m_nodesToFree
Definition: bvh_pbrt.h:147
int flattenBVHTree(BVHBuildNode *node, uint32_t *offset)
Definition: bvh_pbrt.cpp:1034
LinearBVHNode * m_nodes
Definition: bvh_pbrt.h:145
unsigned int m_I[RAYPACKET_RAYS_PER_PACKET]
Definition: bvh_pbrt.h:150
BVHBuildNode * emitLBVH(BVHBuildNode *&buildNodes, const std::vector< BVHPrimitiveInfo > &primitiveInfo, MortonPrimitive *mortonPrims, int nPrimitives, int *totalNodes, CONST_VECTOR_OBJECT &orderedPrims, int *orderedPrimsOffset, int bit)
TODO: after implement memory arena, put const back to this functions.
Definition: bvh_pbrt.cpp:795
BVH_PBRT(const CONTAINER_3D_BASE &aObjectContainer, int aMaxPrimsInNode=4, SPLITMETHOD aSplitMethod=SPLITMETHOD::SAH)
Definition: bvh_pbrt.cpp:235
BVHBuildNode * buildUpperSAH(std::vector< BVHBuildNode * > &treeletRoots, int start, int end, int *totalNodes)
Definition: bvh_pbrt.cpp:894
const int m_maxPrimsInNode
Definition: bvh_pbrt.h:142
bool Intersect(const RAY &aRay, HITINFO &aHitInfo) const override
Definition: bvh_pbrt.cpp:1066
CONST_VECTOR_OBJECT m_primitives
Definition: bvh_pbrt.h:144
bool IntersectP(const RAY &aRay, float aMaxDistance) const override
Definition: bvh_pbrt.cpp:1195
BVHBuildNode * HLBVHBuild(const std::vector< BVHPrimitiveInfo > &primitiveInfo, int *totalNodes, CONST_VECTOR_OBJECT &orderedPrims)
Definition: bvh_pbrt.cpp:678
BVHBuildNode * recursiveBuild(std::vector< BVHPrimitiveInfo > &primitiveInfo, int start, int end, int *totalNodes, CONST_VECTOR_OBJECT &orderedPrims)
Definition: bvh_pbrt.cpp:426
SPLITMETHOD m_splitMethod
Definition: bvh_pbrt.h:143
void ConvertTo(CONST_VECTOR_OBJECT &aOutVector) const
const LIST_OBJECT & GetList() const
Definition: container_3d.h:59
bool GetCastShadows() const
Definition: material.h:309
virtual bool IntersectP(const RAY &aRay, float aMaxDistance) const =0
const MATERIAL * GetMaterial() const
Definition: object_3d.h:64
std::vector< const OBJECT_3D * > CONST_VECTOR_OBJECT
Definition: container_3d.h:38
This file contains miscellaneous commonly used macros and functions.
#define KI_FALLTHROUGH
The KI_FALLTHROUGH macro is to be used when switch statement cases should purposely fallthrough from ...
Definition: macros.h:83
Definition: bitmap.cpp:64
#define RAYPACKET_RAYS_PER_PACKET
Definition: raypacket.h:35
static std::vector< std::string > split(const std::string &aStr, const std::string &aDelim)
Split the input string into a vector of output strings.
Definition: string_utils.h:296
Manage a bounding box defined by two SFVEC3F min max points.
Definition: bbox_3d.h:42
bool Intersect(const RAY &aRay, float *t) const
Definition: bbox_3d_ray.cpp:46
void Union(const SFVEC3F &aPoint)
Recalculate the bounding box adding a point.
Definition: bbox_3d.cpp:102
SFVEC3F GetCenter() const
Return the center point of the bounding box.
Definition: bbox_3d.cpp:132
const SFVEC3F & Min() const
Return the minimum vertex pointer.
Definition: bbox_3d.h:183
SFVEC3F Offset(const SFVEC3F &p) const
Definition: bbox_3d.cpp:251
const SFVEC3F & Max() const
Return the maximum vertex pointer.
Definition: bbox_3d.h:190
void Set(const SFVEC3F &aPbMin, const SFVEC3F &aPbMax)
Set bounding box with new parameters.
Definition: bbox_3d.cpp:68
void Reset()
Reset the bounding box to zero and de-initialize it.
Definition: bbox_3d.cpp:95
unsigned int MaxDimension() const
Definition: bbox_3d.cpp:151
float SurfaceArea() const
Definition: bbox_3d.cpp:174
int firstPrimOffset
Definition: bvh_pbrt.cpp:129
void InitLeaf(int first, int n, const BBOX_3D &b)
Definition: bvh_pbrt.cpp:109
void InitInterior(int axis, BVHBuildNode *c0, BVHBuildNode *c1)
Definition: bvh_pbrt.cpp:117
BBOX_3D bounds
Definition: bvh_pbrt.cpp:127
BVHBuildNode * children[2]
Definition: bvh_pbrt.cpp:128
BVHPrimitiveInfo(int aPrimitiveNumber, const BBOX_3D &aBounds)
Definition: bvh_pbrt.cpp:95
SFVEC3F centroid
Definition: bvh_pbrt.cpp:102
BBOX_3D bounds
Definition: bvh_pbrt.cpp:422
bool operator()(const BVHPrimitiveInfo &a, const BVHPrimitiveInfo &b) const
Definition: bvh_pbrt.cpp:328
ComparePoints(int d)
Definition: bvh_pbrt.cpp:324
const BBOX_3D & centroidBounds
Definition: bvh_pbrt.cpp:363
CompareToBucket(int split, int num, int d, const BBOX_3D &b)
Definition: bvh_pbrt.cpp:351
bool operator()(const BVHPrimitiveInfo &p) const
Definition: bvh_pbrt.cpp:367
CompareToMid(int d, float m)
Definition: bvh_pbrt.cpp:337
bool operator()(const BVHPrimitiveInfo &a) const
Definition: bvh_pbrt.cpp:342
Stores the hit information of a ray with a point on the surface of a object.
Definition: hitinfo.h:36
unsigned int m_acc_node_info
( 4) The acc stores here the node that it hits
Definition: hitinfo.h:42
float m_tHit
( 4) distance
Definition: hitinfo.h:38
HLBVH_SAH_Evaluator(int split, int num, int d, const BBOX_3D &b)
Definition: bvh_pbrt.cpp:387
const BBOX_3D & centroidBounds
Definition: bvh_pbrt.cpp:397
bool operator()(const BVHBuildNode *node) const
Definition: bvh_pbrt.cpp:401
int numPrimitives
Definition: bvh_pbrt.cpp:142
int startIndex
Definition: bvh_pbrt.cpp:142
BVHBuildNode * buildNodes
Definition: bvh_pbrt.cpp:143
int secondChildOffset
interior
Definition: bvh_pbrt.h:91
uint8_t axis
interior node: xyz
Definition: bvh_pbrt.h:96
BBOX_3D bounds
Definition: bvh_pbrt.h:85
uint16_t nPrimitives
0 -> interior node
Definition: bvh_pbrt.h:95
int primitivesOffset
leaf
Definition: bvh_pbrt.h:90
uint32_t mortonCode
Definition: bvh_pbrt.cpp:136
Definition: ray.h:63
unsigned int m_dirIsNeg[3]
Definition: ray.h:75
glm::vec3 SFVEC3F
Definition: xv3d_types.h:44