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