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 <vector>
72 #include <boost/range/algorithm/partition.hpp>
73 #include <boost/range/algorithm/nth_element.hpp>
74 #include <stdlib.h>
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 
278  if( m_splitMethod == SPLIT_HLBVH )
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 SPLIT_MIDDLE: printf( "using SPLIT_MIDDLE\n" ); break;
318  case SPLIT_EQUALCOUNTS: printf( "using SPLIT_EQUALCOUNTS\n" ); break;
319  case SPLIT_SAH: printf( "using SPLIT_SAH\n" ); break;
320  case SPLIT_HLBVH: printf( "using SPLIT_HLBVH\n" ); break;
321  }
322 
323  printf( " BVH created with %d nodes (%.2f MB)\n",
324  totalNodes, float(treeBytes) / (1024.f * 1024.f) );
325  printf( "////////////////////////////////////////////////////////////////////////////////\n\n" );
326 #endif
327 }
328 
329 
331 {
332  if( !m_addresses_pointer_to_mm_free.empty() )
333  {
334  for( std::list<void *>::iterator ii = m_addresses_pointer_to_mm_free.begin();
335  ii != m_addresses_pointer_to_mm_free.end();
336  ++ii )
337  {
338  free( *ii );
339  *ii = NULL;
340  }
341  }
342 
344 }
345 
346 
348 {
349  explicit ComparePoints(int d) { dim = d; }
350 
351  int dim;
352 
353  bool operator()( const BVHPrimitiveInfo &a,
354  const BVHPrimitiveInfo &b ) const
355  {
356  return a.centroid[dim] < b.centroid[dim];
357  }
358 };
359 
360 
362 {
363  explicit CompareToMid( int d, float m ) { dim = d; mid = m; }
364 
365  int dim;
366  float mid;
367 
368  bool operator()( const BVHPrimitiveInfo &a ) const
369  {
370  return a.centroid[dim] < mid;
371  }
372 };
373 
374 
376 {
377  CompareToBucket( int split, int num, int d, const CBBOX &b )
378  : centroidBounds(b)
379  { splitBucket = split; nBuckets = num; dim = d; }
380 
381  bool operator()(const BVHPrimitiveInfo &p) const;
382 
384 
386 };
387 
388 
390 {
391  const float centroid = p.centroid[dim];
392 
393  int b = nBuckets *
394  // Computes the offset (0.0 - 1.0) for one axis
395  ( ( centroid - centroidBounds.Min()[dim] ) /
396  ( centroidBounds.Max()[dim] - centroidBounds.Min()[dim] ) );
397 
398  if( b == nBuckets )
399  b = nBuckets - 1;
400 
401  wxASSERT( (b >= 0) && (b < nBuckets) );
402 
403  return b <= splitBucket;
404 }
405 
406 
408 {
409  HLBVH_SAH_Evaluator( int split, int num, int d, const CBBOX &b )
410  : centroidBounds(b)
411  { minCostSplitBucket = split; nBuckets = num; dim = d; }
412 
413  bool operator()(const BVHBuildNode *node) const;
414 
417 };
418 
419 
421 {
422  const float centroid = node->bounds.GetCenter( dim );
423 
424  int b = nBuckets *
425  // Computes the offset (0.0 - 1.0) for one axis
426  ( ( centroid - centroidBounds.Min()[dim] ) /
427  ( centroidBounds.Max()[dim] - centroidBounds.Min()[dim] ) );
428 
429  if( b == nBuckets )
430  b = nBuckets - 1;
431 
432  wxASSERT( b >= 0 && b < nBuckets );
433 
434  return b <= minCostSplitBucket;
435 }
436 
437 
439 {
440  int count;
442 };
443 
444 
445 BVHBuildNode *CBVH_PBRT::recursiveBuild ( std::vector<BVHPrimitiveInfo> &primitiveInfo,
446  int start,
447  int end,
448  int *totalNodes,
449  CONST_VECTOR_OBJECT &orderedPrims )
450 {
451  wxASSERT( totalNodes != NULL );
452  wxASSERT( start >= 0 );
453  wxASSERT( end >= 0 );
454  wxASSERT( start != end );
455  wxASSERT( start < end );
456  wxASSERT( start <= (int)primitiveInfo.size() );
457  wxASSERT( end <= (int)primitiveInfo.size() );
458 
459  (*totalNodes)++;
460 
461  // !TODO: implement an memory Arena
462  BVHBuildNode *node = static_cast<BVHBuildNode *>( malloc( sizeof( BVHBuildNode ) ) );
463  m_addresses_pointer_to_mm_free.push_back( node );
464 
465  node->bounds.Reset();
466  node->firstPrimOffset = 0;
467  node->nPrimitives = 0;
468  node->splitAxis = 0;
469  node->children[0] = NULL;
470  node->children[1] = NULL;
471 
472  // Compute bounds of all primitives in BVH node
473  CBBOX bounds;
474  bounds.Reset();
475 
476  for( int i = start; i < end; ++i )
477  bounds.Union( primitiveInfo[i].bounds );
478 
479  int nPrimitives = end - start;
480 
481  if( nPrimitives == 1 )
482  {
483  // Create leaf _BVHBuildNode_
484  int firstPrimOffset = orderedPrims.size();
485 
486  for( int i = start; i < end; ++i )
487  {
488  int primitiveNr = primitiveInfo[i].primitiveNumber;
489  wxASSERT( primitiveNr < (int)m_primitives.size() );
490  orderedPrims.push_back( m_primitives[ primitiveNr ] );
491  }
492 
493  node->InitLeaf( firstPrimOffset, nPrimitives, bounds );
494  }
495  else
496  {
497  // Compute bound of primitive centroids, choose split dimension _dim_
498  CBBOX centroidBounds;
499  centroidBounds.Reset();
500 
501  for( int i = start; i < end; ++i )
502  centroidBounds.Union( primitiveInfo[i].centroid );
503 
504  const int dim = centroidBounds.MaxDimension();
505 
506  // Partition primitives into two sets and build children
507  int mid = (start + end) / 2;
508 
509  if( fabs( centroidBounds.Max()[dim] -
510  centroidBounds.Min()[dim] ) < (FLT_EPSILON + FLT_EPSILON) )
511  {
512  // Create leaf _BVHBuildNode_
513  const int firstPrimOffset = orderedPrims.size();
514 
515  for( int i = start; i < end; ++i )
516  {
517  int primitiveNr = primitiveInfo[i].primitiveNumber;
518 
519  wxASSERT( (primitiveNr >= 0) &&
520  (primitiveNr < (int)m_primitives.size()) );
521 
522  const COBJECT *obj = static_cast<const COBJECT *>( m_primitives[ primitiveNr ] );
523 
524  wxASSERT( obj != NULL );
525 
526  orderedPrims.push_back( obj );
527  }
528 
529  node->InitLeaf( firstPrimOffset, nPrimitives, bounds );
530  }
531  else
532  {
533  // Partition primitives based on _splitMethod_
534  switch( m_splitMethod )
535  {
536  case SPLIT_MIDDLE:
537  {
538  // Partition primitives through node's midpoint
539  float pmid = centroidBounds.GetCenter( dim );
540 
541  BVHPrimitiveInfo *midPtr = std::partition( &primitiveInfo[start],
542  &primitiveInfo[end - 1] + 1,
543  CompareToMid( dim, pmid ) );
544  mid = midPtr - &primitiveInfo[0];
545 
546  wxASSERT( (mid >= start) &&
547  (mid <= end) );
548 
549  if( (mid != start) && (mid != end) )
550  // for lots of prims with large overlapping bounding boxes, this
551  // may fail to partition; in that case don't break and fall through
552  // to SPLIT_EQUAL_COUNTS
553  break;
554  }
555 
556  case SPLIT_EQUALCOUNTS:
557  {
558  // Partition primitives into equally-sized subsets
559  mid = (start + end) / 2;
560 
561  std::nth_element( &primitiveInfo[start],
562  &primitiveInfo[mid],
563  &primitiveInfo[end - 1] + 1,
564  ComparePoints( dim ) );
565 
566  break;
567  }
568 
569  case SPLIT_SAH:
570  default:
571  {
572  // Partition primitives using approximate SAH
573  if( nPrimitives <= 2 )
574  {
575  // Partition primitives into equally-sized subsets
576  mid = (start + end) / 2;
577 
578  std::nth_element( &primitiveInfo[start],
579  &primitiveInfo[mid],
580  &primitiveInfo[end - 1] + 1,
581  ComparePoints( dim ) );
582  }
583  else
584  {
585  // Allocate _BucketInfo_ for SAH partition buckets
586  const int nBuckets = 12;
587 
588  BucketInfo buckets[nBuckets];
589 
590  for( int i = 0; i < nBuckets; ++i )
591  {
592  buckets[i].count = 0;
593  buckets[i].bounds.Reset();
594  }
595 
596  // Initialize _BucketInfo_ for SAH partition buckets
597  for( int i = start; i < end; ++i )
598  {
599  int b = nBuckets *
600  centroidBounds.Offset( primitiveInfo[i].centroid )[dim];
601 
602  if( b == nBuckets )
603  b = nBuckets - 1;
604 
605  wxASSERT( b >= 0 && b < nBuckets );
606 
607  buckets[b].count++;
608  buckets[b].bounds.Union( primitiveInfo[i].bounds );
609  }
610 
611  // Compute costs for splitting after each bucket
612  float cost[nBuckets - 1];
613 
614  for( int i = 0; i < (nBuckets - 1); ++i )
615  {
616  CBBOX b0, b1;
617 
618  b0.Reset();
619  b1.Reset();
620 
621  int count0 = 0;
622  int count1 = 0;
623 
624  for( int j = 0; j <= i; ++j )
625  {
626  if( buckets[j].count )
627  {
628  count0 += buckets[j].count;
629  b0.Union( buckets[j].bounds );
630  }
631  }
632 
633  for( int j = i + 1; j < nBuckets; ++j )
634  {
635  if( buckets[j].count )
636  {
637  count1 += buckets[j].count;
638  b1.Union( buckets[j].bounds );
639  }
640  }
641 
642  cost[i] = 1.0f +
643  ( count0 * b0.SurfaceArea() +
644  count1 * b1.SurfaceArea() ) /
645  bounds.SurfaceArea();
646  }
647 
648  // Find bucket to split at that minimizes SAH metric
649  float minCost = cost[0];
650  int minCostSplitBucket = 0;
651 
652  for( int i = 1; i < (nBuckets - 1); ++i )
653  {
654  if( cost[i] < minCost )
655  {
656  minCost = cost[i];
657  minCostSplitBucket = i;
658  }
659  }
660 
661  // Either create leaf or split primitives at selected SAH
662  // bucket
663  if( (nPrimitives > m_maxPrimsInNode) ||
664  (minCost < (float)nPrimitives) )
665  {
666  BVHPrimitiveInfo *pmid =
667  std::partition( &primitiveInfo[start],
668  &primitiveInfo[end - 1] + 1,
669  CompareToBucket( minCostSplitBucket,
670  nBuckets,
671  dim,
672  centroidBounds ) );
673  mid = pmid - &primitiveInfo[0];
674 
675  wxASSERT( (mid >= start) &&
676  (mid <= end) );
677  }
678  else
679  {
680  // Create leaf _BVHBuildNode_
681  const int firstPrimOffset = orderedPrims.size();
682 
683  for( int i = start; i < end; ++i )
684  {
685  const int primitiveNr = primitiveInfo[i].primitiveNumber;
686 
687  wxASSERT( primitiveNr < (int)m_primitives.size() );
688 
689  orderedPrims.push_back( m_primitives[ primitiveNr ] );
690  }
691 
692  node->InitLeaf( firstPrimOffset, nPrimitives, bounds );
693 
694  return node;
695  }
696  }
697  break;
698  }
699  }
700 
701  node->InitInterior( dim,
702  recursiveBuild( primitiveInfo,
703  start,
704  mid,
705  totalNodes,
706  orderedPrims ),
707  recursiveBuild( primitiveInfo,
708  mid,
709  end,
710  totalNodes,
711  orderedPrims) );
712  }
713  }
714 
715  return node;
716 }
717 
718 
719 BVHBuildNode *CBVH_PBRT::HLBVHBuild( const std::vector<BVHPrimitiveInfo> &primitiveInfo,
720  int *totalNodes,
721  CONST_VECTOR_OBJECT &orderedPrims )
722 {
723  // Compute bounding box of all primitive centroids
724  CBBOX bounds;
725  bounds.Reset();
726 
727  for( unsigned int i = 0; i < primitiveInfo.size(); ++i )
728  bounds.Union( primitiveInfo[i].centroid );
729 
730  // Compute Morton indices of primitives
731  std::vector<MortonPrimitive> mortonPrims( primitiveInfo.size() );
732 
733  for( int i = 0; i < (int)primitiveInfo.size(); ++i )
734  {
735  // Initialize _mortonPrims[i]_ for _i_th primitive
736  const int mortonBits = 10;
737  const int mortonScale = 1 << mortonBits;
738 
739  wxASSERT( primitiveInfo[i].primitiveNumber < (int)primitiveInfo.size() );
740 
741  mortonPrims[i].primitiveIndex = primitiveInfo[i].primitiveNumber;
742 
743  const SFVEC3F centroidOffset = bounds.Offset( primitiveInfo[i].centroid );
744 
745  wxASSERT( (centroidOffset.x >= 0.0f) && (centroidOffset.x <= 1.0f) );
746  wxASSERT( (centroidOffset.y >= 0.0f) && (centroidOffset.y <= 1.0f) );
747  wxASSERT( (centroidOffset.z >= 0.0f) && (centroidOffset.z <= 1.0f) );
748 
749  mortonPrims[i].mortonCode = EncodeMorton3( centroidOffset *
750  SFVEC3F( (float)mortonScale ) );
751  }
752 
753  // Radix sort primitive Morton indices
754  RadixSort( &mortonPrims );
755 
756  // Create LBVH treelets at bottom of BVH
757 
758  // Find intervals of primitives for each treelet
759  std::vector<LBVHTreelet> treeletsToBuild;
760 
761  for( int start = 0, end = 1; end <= (int)mortonPrims.size(); ++end )
762  {
763  const uint32_t mask = 0b00111111111111000000000000000000;
764 
765  if( (end == (int)mortonPrims.size()) ||
766  ( (mortonPrims[start].mortonCode & mask) !=
767  (mortonPrims[end].mortonCode & mask) ) )
768  {
769  // Add entry to _treeletsToBuild_ for this treelet
770  const int numPrimitives = end - start;
771  const int maxBVHNodes = 2 * numPrimitives;
772 
773  // !TODO: implement a memory arena
774  BVHBuildNode *nodes = static_cast<BVHBuildNode *>( malloc( maxBVHNodes *
775  sizeof( BVHBuildNode ) ) );
776 
777  m_addresses_pointer_to_mm_free.push_back( nodes );
778 
779  for( int i = 0; i < maxBVHNodes; ++i )
780  {
781  nodes[i].bounds.Reset();
782  nodes[i].firstPrimOffset = 0;
783  nodes[i].nPrimitives = 0;
784  nodes[i].splitAxis = 0;
785  nodes[i].children[0] = NULL;
786  nodes[i].children[1] = NULL;
787  }
788 
789  LBVHTreelet tmpTreelet;
790 
791  tmpTreelet.startIndex = start;
792  tmpTreelet.numPrimitives = numPrimitives;
793  tmpTreelet.buildNodes = nodes;
794 
795  treeletsToBuild.push_back( tmpTreelet );
796 
797  start = end;
798  }
799  }
800 
801  // Create LBVHs for treelets in parallel
802  int atomicTotal = 0;
803  int orderedPrimsOffset = 0;
804 
805  orderedPrims.resize( m_primitives.size() );
806 
807  for( int index = 0; index < (int)treeletsToBuild.size(); ++index )
808  {
809  // Generate _index_th LBVH treelet
810  int nodesCreated = 0;
811  const int firstBit = 29 - 12;
812 
813  LBVHTreelet &tr = treeletsToBuild[index];
814 
815  wxASSERT( tr.startIndex < (int)mortonPrims.size() );
816 
817  tr.buildNodes = emitLBVH( tr.buildNodes,
818  primitiveInfo,
819  &mortonPrims[tr.startIndex],
820  tr.numPrimitives,
821  &nodesCreated,
822  orderedPrims,
823  &orderedPrimsOffset,
824  firstBit );
825 
826  atomicTotal += nodesCreated;
827  }
828 
829  *totalNodes = atomicTotal;
830 
831  // Initialize _finishedTreelets_ with treelet root node pointers
832  std::vector<BVHBuildNode *> finishedTreelets;
833  finishedTreelets.reserve( treeletsToBuild.size() );
834 
835  for( int index = 0; index < (int)treeletsToBuild.size(); ++index )
836  finishedTreelets.push_back( treeletsToBuild[index].buildNodes );
837 
838  // Create and return SAH BVH from LBVH treelets
839  return buildUpperSAH( finishedTreelets,
840  0,
841  finishedTreelets.size(),
842  totalNodes );
843 }
844 
845 
847  BVHBuildNode *&buildNodes,
848  const std::vector<BVHPrimitiveInfo> &primitiveInfo,
849  MortonPrimitive *mortonPrims, int nPrimitives, int *totalNodes,
850  CONST_VECTOR_OBJECT &orderedPrims,
851  int *orderedPrimsOffset, int bit)
852 {
853  wxASSERT( nPrimitives > 0 );
854  wxASSERT( totalNodes != NULL );
855  wxASSERT( orderedPrimsOffset != NULL );
856  wxASSERT( nPrimitives > 0 );
857  wxASSERT( mortonPrims != NULL );
858 
859  if( (bit == -1) || (nPrimitives < m_maxPrimsInNode) )
860  {
861  // Create and return leaf node of LBVH treelet
862  (*totalNodes)++;
863 
864  BVHBuildNode *node = buildNodes++;
865  CBBOX bounds;
866  bounds.Reset();
867 
868  int firstPrimOffset = *orderedPrimsOffset;
869  *orderedPrimsOffset += nPrimitives;
870 
871  wxASSERT( (firstPrimOffset + (nPrimitives - 1)) < (int)orderedPrims.size() );
872 
873  for( int i = 0; i < nPrimitives; ++i )
874  {
875  const int primitiveIndex = mortonPrims[i].primitiveIndex;
876 
877  wxASSERT( primitiveIndex < (int)m_primitives.size() );
878 
879  orderedPrims[firstPrimOffset + i] = m_primitives[primitiveIndex];
880  bounds.Union( primitiveInfo[primitiveIndex].bounds );
881  }
882 
883  node->InitLeaf( firstPrimOffset, nPrimitives, bounds );
884 
885  return node;
886  }
887  else
888  {
889  int mask = 1 << bit;
890 
891  // Advance to next subtree level if there's no LBVH split for this bit
892  if( (mortonPrims[0].mortonCode & mask) ==
893  (mortonPrims[nPrimitives - 1].mortonCode & mask) )
894  return emitLBVH( buildNodes, primitiveInfo, mortonPrims, nPrimitives,
895  totalNodes, orderedPrims, orderedPrimsOffset,
896  bit - 1 );
897 
898  // Find LBVH split point for this dimension
899  int searchStart = 0;
900  int searchEnd = nPrimitives - 1;
901 
902  while( searchStart + 1 != searchEnd )
903  {
904  wxASSERT(searchStart != searchEnd);
905 
906  const int mid = (searchStart + searchEnd) / 2;
907 
908  if( (mortonPrims[searchStart].mortonCode & mask) ==
909  (mortonPrims[mid].mortonCode & mask) )
910  searchStart = mid;
911  else
912  {
913  wxASSERT( (mortonPrims[mid].mortonCode & mask) ==
914  (mortonPrims[searchEnd].mortonCode & mask) );
915  searchEnd = mid;
916  }
917  }
918 
919  const int splitOffset = searchEnd;
920 
921  wxASSERT( splitOffset <= (nPrimitives - 1) );
922  wxASSERT( (mortonPrims[splitOffset - 1].mortonCode & mask) !=
923  (mortonPrims[splitOffset].mortonCode & mask) );
924 
925  // Create and return interior LBVH node
926  (*totalNodes)++;
927 
928  BVHBuildNode *node = buildNodes++;
929  BVHBuildNode *lbvh[2];
930 
931  lbvh[0] = emitLBVH( buildNodes, primitiveInfo, mortonPrims, splitOffset,
932  totalNodes, orderedPrims, orderedPrimsOffset, bit - 1 );
933 
934  lbvh[1] = emitLBVH( buildNodes, primitiveInfo, &mortonPrims[splitOffset],
935  nPrimitives - splitOffset, totalNodes, orderedPrims,
936  orderedPrimsOffset, bit - 1 );
937 
938  const int axis = bit % 3;
939 
940  node->InitInterior( axis, lbvh[0], lbvh[1] );
941 
942  return node;
943  }
944 }
945 
946 
948  std::vector<BVHBuildNode *> &treeletRoots,
949  int start, int end,
950  int *totalNodes )
951 {
952  wxASSERT( totalNodes != NULL );
953  wxASSERT( start < end );
954  wxASSERT( end <= (int)treeletRoots.size() );
955 
956  int nNodes = end - start;
957 
958  if( nNodes == 1 )
959  return treeletRoots[start];
960 
961 
962  (*totalNodes)++;
963 
964  BVHBuildNode *node = static_cast<BVHBuildNode *>( malloc( sizeof( BVHBuildNode ) ) );
965 
966  m_addresses_pointer_to_mm_free.push_back( node );
967 
968  node->bounds.Reset();
969  node->firstPrimOffset = 0;
970  node->nPrimitives = 0;
971  node->splitAxis = 0;
972  node->children[0] = NULL;
973  node->children[1] = NULL;
974 
975  // Compute bounds of all nodes under this HLBVH node
976  CBBOX bounds;
977  bounds.Reset();
978 
979  for( int i = start; i < end; ++i )
980  bounds.Union( treeletRoots[i]->bounds );
981 
982  // Compute bound of HLBVH node centroids, choose split dimension _dim_
983  CBBOX centroidBounds;
984  centroidBounds.Reset();
985 
986  for( int i = start; i < end; ++i )
987  {
988  SFVEC3F centroid =
989  (treeletRoots[i]->bounds.Min() + treeletRoots[i]->bounds.Max()) *
990  0.5f;
991 
992  centroidBounds.Union(centroid);
993  }
994 
995  const int dim = centroidBounds.MaxDimension();
996 
997  // FIXME: if this hits, what do we need to do?
998  // Make sure the SAH split below does something... ?
999  wxASSERT( centroidBounds.Max()[dim] != centroidBounds.Min()[dim] );
1000 
1001  // Allocate _BucketInfo_ for SAH partition buckets
1002  const int nBuckets = 12;
1003 
1004  BucketInfo buckets[nBuckets];
1005 
1006  for( int i = 0; i < nBuckets; ++i )
1007  {
1008  buckets[i].count = 0;
1009  buckets[i].bounds.Reset();
1010  }
1011 
1012  // Initialize _BucketInfo_ for HLBVH SAH partition buckets
1013  for( int i = start; i < end; ++i )
1014  {
1015  const float centroid = ( treeletRoots[i]->bounds.Min()[dim] +
1016  treeletRoots[i]->bounds.Max()[dim] ) *
1017  0.5f;
1018  int b =
1019  nBuckets * ( (centroid - centroidBounds.Min()[dim] ) /
1020  (centroidBounds.Max()[dim] - centroidBounds.Min()[dim] ) );
1021 
1022  if( b == nBuckets )
1023  b = nBuckets - 1;
1024 
1025  wxASSERT( (b >= 0) && (b < nBuckets) );
1026 
1027  buckets[b].count++;
1028  buckets[b].bounds.Union( treeletRoots[i]->bounds );
1029  }
1030 
1031  // Compute costs for splitting after each bucket
1032  float cost[nBuckets - 1];
1033 
1034  for( int i = 0; i < nBuckets - 1; ++i )
1035  {
1036  CBBOX b0, b1;
1037  b0.Reset();
1038  b1.Reset();
1039 
1040  int count0 = 0, count1 = 0;
1041 
1042  for( int j = 0; j <= i; ++j )
1043  {
1044  if( buckets[j].count )
1045  {
1046  count0 += buckets[j].count;
1047  b0.Union( buckets[j].bounds );
1048  }
1049  }
1050 
1051  for( int j = i + 1; j < nBuckets; ++j )
1052  {
1053  if( buckets[j].count )
1054  {
1055  count1 += buckets[j].count;
1056  b1.Union( buckets[j].bounds );
1057  }
1058  }
1059 
1060  cost[i] = .125f +
1061  ( count0 * b0.SurfaceArea() + count1 * b1.SurfaceArea() ) /
1062  bounds.SurfaceArea();
1063  }
1064 
1065  // Find bucket to split at that minimizes SAH metric
1066  float minCost = cost[0];
1067  int minCostSplitBucket = 0;
1068 
1069  for( int i = 1; i < nBuckets - 1; ++i )
1070  {
1071  if( cost[i] < minCost )
1072  {
1073  minCost = cost[i];
1074  minCostSplitBucket = i;
1075  }
1076  }
1077 
1078  // Split nodes and create interior HLBVH SAH node
1079  BVHBuildNode **pmid = std::partition( &treeletRoots[start],
1080  &treeletRoots[end - 1] + 1,
1081  HLBVH_SAH_Evaluator( minCostSplitBucket,
1082  nBuckets,
1083  dim,
1084  centroidBounds ) );
1085 
1086  const int mid = pmid - &treeletRoots[0];
1087 
1088  wxASSERT( (mid > start) && (mid < end) );
1089 
1090  node->InitInterior( dim,
1091  buildUpperSAH( treeletRoots, start, mid, totalNodes ),
1092  buildUpperSAH( treeletRoots, mid, end, totalNodes ) );
1093 
1094  return node;
1095 }
1096 
1097 
1098 int CBVH_PBRT::flattenBVHTree( BVHBuildNode *node, uint32_t *offset )
1099 {
1100  LinearBVHNode *linearNode = &m_nodes[*offset];
1101 
1102  linearNode->bounds = node->bounds;
1103 
1104  int myOffset = (*offset)++;
1105 
1106  if( node->nPrimitives > 0 )
1107  {
1108  wxASSERT( (!node->children[0]) && (!node->children[1]) );
1109  wxASSERT( node->nPrimitives < 65536 );
1110 
1111  linearNode->primitivesOffset = node->firstPrimOffset;
1112  linearNode->nPrimitives = node->nPrimitives;
1113  }
1114  else
1115  {
1116  // Creater interior flattened BVH node
1117  linearNode->axis = node->splitAxis;
1118  linearNode->nPrimitives = 0;
1119  flattenBVHTree( node->children[0], offset );
1120  linearNode->secondChildOffset = flattenBVHTree( node->children[1], offset );
1121  }
1122 
1123  return myOffset;
1124 }
1125 
1126 
1127 #define MAX_TODOS 64
1128 
1129 bool CBVH_PBRT::Intersect( const RAY &aRay, HITINFO &aHitInfo ) const
1130 {
1131  if( !m_nodes )
1132  return false;
1133 
1134  bool hit = false;
1135 
1136  // Follow ray through BVH nodes to find primitive intersections
1137  int todoOffset = 0, nodeNum = 0;
1138  int todo[MAX_TODOS];
1139 
1140  while( true )
1141  {
1142  const LinearBVHNode *node = &m_nodes[nodeNum];
1143 
1144  wxASSERT( todoOffset < MAX_TODOS );
1145 
1146  // Check ray against BVH node
1147  float hitBox = 0.0f;
1148 
1149  const bool hitted = node->bounds.Intersect( aRay, &hitBox );
1150 
1151  if( hitted && (hitBox < aHitInfo.m_tHit) )
1152  {
1153  if( node->nPrimitives > 0 )
1154  {
1155  // Intersect ray with primitives in leaf BVH node
1156  for( int i = 0; i < node->nPrimitives; ++i )
1157  {
1158  if( m_primitives[node->primitivesOffset + i]->Intersect( aRay,
1159  aHitInfo ) )
1160  {
1161  aHitInfo.m_acc_node_info = nodeNum;
1162  hit = true;
1163  }
1164  }
1165  }
1166  else
1167  {
1168  // Put far BVH node on _todo_ stack, advance to near node
1169  if( aRay.m_dirIsNeg[node->axis] )
1170  {
1171  todo[todoOffset++] = nodeNum + 1;
1172  nodeNum = node->secondChildOffset;
1173  }
1174  else
1175  {
1176  todo[todoOffset++] = node->secondChildOffset;
1177  nodeNum = nodeNum + 1;
1178  }
1179 
1180  continue;
1181  }
1182  }
1183 
1184  if( todoOffset == 0 )
1185  break;
1186 
1187  nodeNum = todo[--todoOffset];
1188  }
1189 
1190  return hit;
1191 }
1192 
1193 // !TODO: this may be optimized
1194 bool CBVH_PBRT::Intersect( const RAY &aRay,
1195  HITINFO &aHitInfo,
1196  unsigned int aAccNodeInfo ) const
1197 {
1198  if( !m_nodes )
1199  return false;
1200 
1201  bool hit = false;
1202 
1203  // Follow ray through BVH nodes to find primitive intersections
1204  int todoOffset = 0, nodeNum = aAccNodeInfo;
1205  int todo[MAX_TODOS];
1206 
1207  while( true )
1208  {
1209  const LinearBVHNode *node = &m_nodes[nodeNum];
1210 
1211  wxASSERT( todoOffset < MAX_TODOS );
1212 
1213  // Check ray against BVH node
1214  float hitBox = 0.0f;
1215 
1216  const bool hitted = node->bounds.Intersect( aRay, &hitBox );
1217 
1218  if( hitted && (hitBox < aHitInfo.m_tHit) )
1219  {
1220  if( node->nPrimitives > 0 )
1221  {
1222  // Intersect ray with primitives in leaf BVH node
1223  for( int i = 0; i < node->nPrimitives; ++i )
1224  {
1225  if( m_primitives[node->primitivesOffset + i]->Intersect( aRay,
1226  aHitInfo ) )
1227  {
1228  //aHitInfo.m_acc_node_info = nodeNum;
1229  hit = true;
1230  }
1231  }
1232  }
1233  else
1234  {
1235  // Put far BVH node on _todo_ stack, advance to near node
1236  if( aRay.m_dirIsNeg[node->axis] )
1237  {
1238  todo[todoOffset++] = nodeNum + 1;
1239  nodeNum = node->secondChildOffset;
1240  }
1241  else
1242  {
1243  todo[todoOffset++] = node->secondChildOffset;
1244  nodeNum = nodeNum + 1;
1245  }
1246 
1247  continue;
1248  }
1249  }
1250 
1251  if( todoOffset == 0 )
1252  break;
1253 
1254  nodeNum = todo[--todoOffset];
1255  }
1256 
1257  return hit;
1258 }
1259 
1260 
1261 bool CBVH_PBRT::IntersectP( const RAY &aRay, float aMaxDistance ) const
1262 {
1263  if( !m_nodes )
1264  return false;
1265 
1266  // Follow ray through BVH nodes to find primitive intersections
1267  int todoOffset = 0, nodeNum = 0;
1268  int todo[MAX_TODOS];
1269 
1270  while( true )
1271  {
1272  const LinearBVHNode *node = &m_nodes[nodeNum];
1273 
1274  wxASSERT( todoOffset < MAX_TODOS );
1275 
1276  // Check ray against BVH node
1277  float hitBox = 0.0f;
1278 
1279  const bool hitted = node->bounds.Intersect( aRay, &hitBox );
1280 
1281  if( hitted && (hitBox < aMaxDistance) )
1282  {
1283  if( node->nPrimitives > 0 )
1284  {
1285  // Intersect ray with primitives in leaf BVH node
1286  for( int i = 0; i < node->nPrimitives; ++i )
1287  {
1288  const COBJECT *obj = m_primitives[node->primitivesOffset + i];
1289 
1290  if( obj->GetMaterial()->GetCastShadows() )
1291  if( obj->IntersectP( aRay, aMaxDistance ) )
1292  return true;
1293  }
1294  }
1295  else
1296  {
1297  // Put far BVH node on _todo_ stack, advance to near node
1298  if( aRay.m_dirIsNeg[node->axis] )
1299  {
1300  todo[todoOffset++] = nodeNum + 1;
1301  nodeNum = node->secondChildOffset;
1302  }
1303  else
1304  {
1305  todo[todoOffset++] = node->secondChildOffset;
1306  nodeNum = nodeNum + 1;
1307  }
1308 
1309  continue;
1310  }
1311  }
1312 
1313  if( todoOffset == 0 )
1314  break;
1315 
1316  nodeNum = todo[--todoOffset];
1317  }
1318 
1319  return false;
1320 }
uint32_t EncodeMorton3(const SFVEC3F &v)
Definition: cbvh_pbrt.cpp:166
CONST_VECTOR_OBJECT m_primitives
Definition: cbvh_pbrt.h:159
float SurfaceArea() const
Function SurfaceArea.
Definition: cbbox.cpp:184
int primitivesOffset
leaf
Definition: cbvh_pbrt.h:90
uint16_t nPrimitives
0 -> interior node
Definition: cbvh_pbrt.h:95
bool operator()(const BVHPrimitiveInfo &p) const
Definition: cbvh_pbrt.cpp:389
std::vector< const COBJECT * > CONST_VECTOR_OBJECT
Definition: ccontainer.h:39
uint8_t axis
interior node: xyz
Definition: cbvh_pbrt.h:96
const SFVEC3F & Min() const
Function Min return the minimun vertex pointer.
Definition: cbbox.h:205
CBVH_PBRT(const CGENERICCONTAINER &aObjectContainer, int aMaxPrimsInNode=4, SPLITMETHOD aSplitMethod=SPLIT_SAH)
Definition: cbvh_pbrt.cpp:233
const CMATERIAL * GetMaterial() const
Definition: cobject.h:63
SFVEC3F Offset(const SFVEC3F &p) const
Function Offset.
Definition: cbbox.cpp:263
int firstPrimOffset
Definition: cbvh_pbrt.cpp:127
Definition: ray.h:43
bool operator()(const BVHPrimitiveInfo &a) const
Definition: cbvh_pbrt.cpp:368
HLBVH_SAH_Evaluator(int split, int num, int d, const CBBOX &b)
Definition: cbvh_pbrt.cpp:409
const int m_maxPrimsInNode
Definition: cbvh_pbrt.h:157
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
BVHPrimitiveInfo(int aPrimitiveNumber, const CBBOX &aBounds)
Definition: cbvh_pbrt.cpp:93
CBBOX bounds
Definition: cbvh_pbrt.cpp:441
ComparePoints(int d)
Definition: cbvh_pbrt.cpp:349
void ConvertTo(CONST_VECTOR_OBJECT &aOutVector) const
Definition: ccontainer.cpp:64
bool IntersectP(const RAY &aRay, float aMaxDistance) const override
Definition: cbvh_pbrt.cpp:1261
void InitLeaf(int first, int n, const CBBOX &b)
Definition: cbvh_pbrt.cpp:107
int numPrimitives
Definition: cbvh_pbrt.cpp:140
bool GetCastShadows() const
Definition: cmaterial.h:199
unsigned int m_dirIsNeg[3]
Definition: ray.h:56
BVHBuildNode * buildUpperSAH(std::vector< BVHBuildNode * > &treeletRoots, int start, int end, int *totalNodes)
Definition: cbvh_pbrt.cpp:947
CompareToBucket(int split, int num, int d, const CBBOX &b)
Definition: cbvh_pbrt.cpp:377
BVHBuildNode * children[2]
Definition: cbvh_pbrt.cpp:126
bool Intersect(const RAY &aRay, HITINFO &aHitInfo) const override
Definition: cbvh_pbrt.cpp:1129
virtual bool IntersectP(const RAY &aRay, float aMaxDistance) const =0
Functions Intersect for shadow test.
const SFVEC3F & Max() const
Function Max return the maximum vertex pointer.
Definition: cbbox.h:212
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:846
bool Intersect(const RAY &aRay, float *t) const
Function Intersect.
Definition: cbbox_ray.cpp:46
BVHBuildNode * recursiveBuild(std::vector< BVHPrimitiveInfo > &primitiveInfo, int start, int end, int *totalNodes, CONST_VECTOR_OBJECT &orderedPrims)
Definition: cbvh_pbrt.cpp:445
BVHBuildNode * buildNodes
Definition: cbvh_pbrt.cpp:141
#define MAX_TODOS
Definition: cbvh_pbrt.cpp:1127
SFVEC3F GetCenter() const
Function GetCenter return the center point of the bounding box.
Definition: cbbox.cpp:135
bool operator()(const BVHBuildNode *node) const
Definition: cbvh_pbrt.cpp:420
LinearBVHNode * m_nodes
Definition: cbvh_pbrt.h:160
uint32_t LeftShift3(uint32_t x)
Definition: cbvh_pbrt.cpp:146
unsigned int MaxDimension() const
Function MaxDimension.
Definition: cbbox.cpp:154
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:363
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:165
static long bit[8][9]
Definition: solve.cpp:826
int flattenBVHTree(BVHBuildNode *node, uint32_t *offset)
Definition: cbvh_pbrt.cpp:1098
bool operator()(const BVHPrimitiveInfo &a, const BVHPrimitiveInfo &b) const
Definition: cbvh_pbrt.cpp:353
BVHBuildNode * HLBVHBuild(const std::vector< BVHPrimitiveInfo > &primitiveInfo, int *totalNodes, CONST_VECTOR_OBJECT &orderedPrims)
Definition: cbvh_pbrt.cpp:719
int secondChildOffset
interior
Definition: cbvh_pbrt.h:91
SPLITMETHOD
Definition: cbvh_pbrt.h:101
CBBOX bounds
Definition: cbvh_pbrt.h:85
Class 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...
const LIST_OBJECT & GetList() const
Definition: ccontainer.h:63
const CBBOX & centroidBounds
Definition: cbvh_pbrt.cpp:385
const CBBOX & centroidBounds
Definition: cbvh_pbrt.cpp:416
#define min(a, b)
Definition: auxiliary.h:85
std::list< void * > m_addresses_pointer_to_mm_free
Definition: cbvh_pbrt.h:162
SPLITMETHOD m_splitMethod
Definition: cbvh_pbrt.h:158