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  {
90  primitiveNumber = 0;
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
148 inline 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 
168 inline 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 
178 static 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 
235 BVH_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();
293  m_nodes[i].primitivesOffset = 0;
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  {
354  splitBucket = split;
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] ) /
374  ( centroidBounds.Max()[dim] - 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] ) /
408  ( centroidBounds.Max()[dim] - 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 
426 BVHBuildNode *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  {
514  case SPLITMETHOD::MIDDLE:
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 
678 BVHBuildNode *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 
894 BVHBuildNode *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 
1034 int 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 
1066 bool 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 
1131 bool 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 
1195 bool 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 }
bool operator()(const BVHPrimitiveInfo &a, const BVHPrimitiveInfo &b) const
Definition: bvh_pbrt.cpp:328
unsigned int MaxDimension() const
Definition: bbox_3d.cpp:151
BBOX_3D bounds
Definition: bvh_pbrt.cpp:422
int primitivesOffset
leaf
Definition: bvh_pbrt.h:90
const SFVEC3F & Max() const
Return the maximum vertex pointer.
Definition: bbox_3d.h:198
CONST_VECTOR_OBJECT m_primitives
Definition: bvh_pbrt.h:144
void InitLeaf(int first, int n, const BBOX_3D &b)
Definition: bvh_pbrt.cpp:109
SPLITMETHOD m_splitMethod
Definition: bvh_pbrt.h:143
uint16_t nPrimitives
0 -> interior node
Definition: bvh_pbrt.h:95
BVHBuildNode * recursiveBuild(std::vector< BVHPrimitiveInfo > &primitiveInfo, int start, int end, int *totalNodes, CONST_VECTOR_OBJECT &orderedPrims)
Definition: bvh_pbrt.cpp:426
Manage a bounding box defined by two SFVEC3F min max points.
Definition: bbox_3d.h:41
bool Intersect(const RAY &aRay, float *t) const
Definition: bbox_3d_ray.cpp:46
uint8_t axis
interior node: xyz
Definition: bvh_pbrt.h:96
BVHBuildNode * HLBVHBuild(const std::vector< BVHPrimitiveInfo > &primitiveInfo, int *totalNodes, CONST_VECTOR_OBJECT &orderedPrims)
Definition: bvh_pbrt.cpp:678
const LIST_OBJECT & GetList() const
Definition: container_3d.h:59
Definition: bitmap.cpp:64
std::vector< const OBJECT_3D * > CONST_VECTOR_OBJECT
Definition: container_3d.h:38
#define KI_FALLTHROUGH
The KI_FALLTHROUGH macro is to be used when switch statement cases should purposely fallthrough from ...
Definition: macros.h:83
int firstPrimOffset
Definition: bvh_pbrt.cpp:129
Definition: ray.h:62
This BVH implementation is based on the source code implementation from the book "Physically Based Re...
SPLITMETHOD
Definition: bvh_pbrt.h:101
SFVEC3F Offset(const SFVEC3F &p) const
Definition: bbox_3d.cpp:251
bool operator()(const BVHPrimitiveInfo &p) const
Definition: bvh_pbrt.cpp:367
int startIndex
Definition: bvh_pbrt.cpp:142
float m_tHit
( 4) distance
Definition: hitinfo.h:38
void InitInterior(int axis, BVHBuildNode *c0, BVHBuildNode *c1)
Definition: bvh_pbrt.cpp:117
void Union(const SFVEC3F &aPoint)
Recalculate the bounding box adding a point.
Definition: bbox_3d.cpp:102
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
static void RadixSort(std::vector< MortonPrimitive > *v)
Definition: bvh_pbrt.cpp:178
float SurfaceArea() const
Definition: bbox_3d.cpp:174
const SFVEC3F & Min() const
Return the minimum vertex pointer.
Definition: bbox_3d.h:191
This file contains miscellaneous commonly used macros and functions.
SFVEC3F GetCenter() const
Return the center point of the bounding box.
Definition: bbox_3d.cpp:132
ComparePoints(int d)
Definition: bvh_pbrt.cpp:324
BBOX_3D bounds
Definition: bvh_pbrt.cpp:127
BVHBuildNode * buildUpperSAH(std::vector< BVHBuildNode * > &treeletRoots, int start, int end, int *totalNodes)
Definition: bvh_pbrt.cpp:894
int numPrimitives
Definition: bvh_pbrt.cpp:142
const int m_maxPrimsInNode
Definition: bvh_pbrt.h:142
unsigned int m_dirIsNeg[3]
Definition: ray.h:75
BVHBuildNode * children[2]
Definition: bvh_pbrt.cpp:128
bool IntersectP(const RAY &aRay, float aMaxDistance) const override
Definition: bvh_pbrt.cpp:1195
void ConvertTo(CONST_VECTOR_OBJECT &aOutVector) const
BVHBuildNode * buildNodes
Definition: bvh_pbrt.cpp:143
uint32_t EncodeMorton3(const SFVEC3F &v)
Definition: bvh_pbrt.cpp:168
HLBVH_SAH_Evaluator(int split, int num, int d, const BBOX_3D &b)
Definition: bvh_pbrt.cpp:387
virtual bool IntersectP(const RAY &aRay, float aMaxDistance) const =0
void Set(const SFVEC3F &aPbMin, const SFVEC3F &aPbMax)
Set bounding box with new parameters.
Definition: bbox_3d.cpp:68
uint32_t LeftShift3(uint32_t x)
Definition: bvh_pbrt.cpp:148
uint32_t mortonCode
Definition: bvh_pbrt.cpp:136
const BBOX_3D & centroidBounds
Definition: bvh_pbrt.cpp:363
std::list< void * > m_nodesToFree
Definition: bvh_pbrt.h:147
BVHPrimitiveInfo(int aPrimitiveNumber, const BBOX_3D &aBounds)
Definition: bvh_pbrt.cpp:95
Stores the hit information of a ray with a point on the surface of a object.
Definition: hitinfo.h:35
#define RAYPACKET_RAYS_PER_PACKET
Definition: raypacket.h:35
CompareToMid(int d, float m)
Definition: bvh_pbrt.cpp:337
unsigned int m_acc_node_info
( 4) The acc stores here the node that it hits
Definition: hitinfo.h:42
glm::vec3 SFVEC3F
Definition: xv3d_types.h:44
bool Intersect(const RAY &aRay, HITINFO &aHitInfo) const override
Definition: bvh_pbrt.cpp:1066
SFVEC3F centroid
Definition: bvh_pbrt.cpp:102
int flattenBVHTree(BVHBuildNode *node, uint32_t *offset)
Definition: bvh_pbrt.cpp:1034
int secondChildOffset
interior
Definition: bvh_pbrt.h:91
BVH_PBRT(const CONTAINER_3D_BASE &aObjectContainer, int aMaxPrimsInNode=4, SPLITMETHOD aSplitMethod=SPLITMETHOD::SAH)
Definition: bvh_pbrt.cpp:235
#define MAX_TODOS
Definition: bvh_pbrt.cpp:1063
bool operator()(const BVHPrimitiveInfo &a) const
Definition: bvh_pbrt.cpp:342
CompareToBucket(int split, int num, int d, const BBOX_3D &b)
Definition: bvh_pbrt.cpp:351
void Reset()
Reset the bounding box to zero and de-initialize it.
Definition: bbox_3d.cpp:95
const BBOX_3D & centroidBounds
Definition: bvh_pbrt.cpp:397
unsigned int m_I[RAYPACKET_RAYS_PER_PACKET]
Definition: bvh_pbrt.h:150
bool GetCastShadows() const
Definition: material.h:309
LinearBVHNode * m_nodes
Definition: bvh_pbrt.h:145
bool operator()(const BVHBuildNode *node) const
Definition: bvh_pbrt.cpp:401
const MATERIAL * GetMaterial() const
Definition: object_3d.h:64
BBOX_3D bounds
Definition: bvh_pbrt.h:85
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:293