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