src/BVH/BVH_LinearBuilder.lxx

   1 // Created on: 2014-09-11
   2 // Created by: Danila ULYANOV
   3 // Copyright (c) 2013-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   7 // This library is free software; you can redistribute it and/or modify it under
   8 // the terms of the GNU Lesser General Public License version 2.1 as published
   9 // by the Free Software Foundation, with special exception defined in the file
  10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  11 // distribution for complete text of the license and disclaimer of any warranty.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16 #include <algorithm>
  17
  18 #include <Standard_Assert.hxx>
  19
  20 #include <NCollection_Array1.hxx>
  21
  22 #ifdef HAVE_TBB
  23   // On Windows, function TryEnterCriticalSection has appeared in Windows NT
  24   // and is surrounded by #ifdef in MS VC++ 7.1 headers.
  25   // Thus to use it we need to define appropriate macro saying that we will
  26   // run on Windows NT 4.0 at least
  27   #if defined(_WIN32) && !defined(_WIN32_WINNT)
  28     #define _WIN32_WINNT 0x0501
  29   #endif
  30
  31   #include <tbb/task.h>
  32 #endif
  33
  34 // =======================================================================
  35 // function : BVH_LinearBuilder
  36 // purpose  :
  37 // =======================================================================
  38 template<class T, int N>
  39 BVH_LinearBuilder<T, N>::BVH_LinearBuilder (const Standard_Integer theLeafNodeSize,
  40                                             const Standard_Integer theMaxTreeDepth)
  41 : BVH_Builder<T, N> (theLeafNodeSize,
  42                      theMaxTreeDepth)
  43 {
  44   //
  45 }
  46
  47 // =======================================================================
  48 // function : ~BVH_LinearBuilder
  49 // purpose  :
  50 // =======================================================================
  51 template<class T, int N>
  52 BVH_LinearBuilder<T, N>::~BVH_LinearBuilder()
  53 {
  54   //
  55 }
  56
  57 namespace BVH
  58 {
  59   // Radix sort STL predicate for 32-bit integer.
  60   class BitPredicate
  61   {
  62     Standard_Integer myBit;
  63
  64   public:
  65
  66     //! Creates new radix sort predicate.
  67     BitPredicate (const Standard_Integer theBit) : myBit (theBit)
  68     {
  69       //
  70     }
  71
  72     //! Returns predicate value.
  73     bool operator() (const BVH_EncodedLink theLink) const
  74     {
  75       const Standard_Integer aMask = 1 << myBit;
  76
  77       return !(theLink.first & aMask); // 0-bit to the left side
  78     }
  79   };
  80
  81   //! STL compare tool used in binary search algorithm.
  82   class BitComparator
  83   {
  84     Standard_Integer myBit;
  85
  86   public:
  87
  88     //! Creates new STL comparator.
  89     BitComparator (const Standard_Integer theBit) : myBit (theBit)
  90     {
  91       //
  92     }
  93
  94     //! Checks left value for the given bit.
  95     bool operator() (BVH_EncodedLink theLink1, BVH_EncodedLink /*theLink2*/)
  96     {
  97       return !(theLink1.first & (1 << myBit));
  98     }
  99   };
 100
 101   //! Tool object for sorting link array using radix sort algorithm.
 102   struct RadixSorter
 103   {
 104     typedef std::vector<BVH_EncodedLink>::iterator LinkIterator;
 105
 106     // Performs MSD (most significant digit) radix sort.
 107     static void Perform (LinkIterator theStart, LinkIterator theFinal, Standard_Integer theBit = 29)
 108     {
 109       while (theStart != theFinal && theBit >= 0)
 110       {
 111         LinkIterator anOffset = std::partition (theStart, theFinal, BitPredicate (theBit--));
 112
 113         Perform (theStart, anOffset, theBit);
 114
 115         theStart = anOffset;
 116       }
 117     }
 118   };
 119
 120   //! Calculates bounding boxes (AABBs) for the given BVH tree.
 121   template<class T, int N>
 122   Standard_Integer UpdateBounds (BVH_Set<T, N>* theSet, BVH_Tree<T, N>* theTree, const Standard_Integer theNode = 0)
 123   {
 124     const BVH_Vec4i aData = theTree->NodeInfoBuffer()[theNode];
 125
 126     if (aData.x() == 0)
 127     {
 128       const Standard_Integer aLftChild = theTree->NodeInfoBuffer()[theNode].y();
 129       const Standard_Integer aRghChild = theTree->NodeInfoBuffer()[theNode].z();
 130
 131       const Standard_Integer aLftDepth = UpdateBounds (theSet, theTree, aLftChild);
 132       const Standard_Integer aRghDepth = UpdateBounds (theSet, theTree, aRghChild);
 133
 134       typename BVH_Box<T, N>::BVH_VecNt aLftMinPoint = theTree->MinPointBuffer()[aLftChild];
 135       typename BVH_Box<T, N>::BVH_VecNt aLftMaxPoint = theTree->MaxPointBuffer()[aLftChild];
 136       typename BVH_Box<T, N>::BVH_VecNt aRghMinPoint = theTree->MinPointBuffer()[aRghChild];
 137       typename BVH_Box<T, N>::BVH_VecNt aRghMaxPoint = theTree->MaxPointBuffer()[aRghChild];
 138
 139       BVH::BoxMinMax<T, N>::CwiseMin (aLftMinPoint, aRghMinPoint);
 140       BVH::BoxMinMax<T, N>::CwiseMax (aLftMaxPoint, aRghMaxPoint);
 141
 142       theTree->MinPointBuffer()[theNode] = aLftMinPoint;
 143       theTree->MaxPointBuffer()[theNode] = aLftMaxPoint;
 144
 145       return Max (aLftDepth, aRghDepth) + 1;
 146     }
 147     else
 148     {
 149       typename BVH_Box<T, N>::BVH_VecNt& aMinPoint = theTree->MinPointBuffer()[theNode];
 150       typename BVH_Box<T, N>::BVH_VecNt& aMaxPoint = theTree->MaxPointBuffer()[theNode];
 151
 152       for (Standard_Integer aPrimIdx = aData.y(); aPrimIdx <= aData.z(); ++aPrimIdx)
 153       {
 154         const BVH_Box<T, N> aBox = theSet->Box (aPrimIdx);
 155
 156         if (aPrimIdx == aData.y())
 157         {
 158           aMinPoint = aBox.CornerMin();
 159           aMaxPoint = aBox.CornerMax();
 160         }
 161         else
 162         {
 163           BVH::BoxMinMax<T, N>::CwiseMin (aMinPoint, aBox.CornerMin());
 164           BVH::BoxMinMax<T, N>::CwiseMax (aMaxPoint, aBox.CornerMax());
 165         }
 166       }
 167     }
 168
 169     return 0;
 170   }
 171 }
 172
 173 // =======================================================================
 174 // function : EmitHierachy
 175 // purpose  : Emits hierarchy from sorted Morton codes
 176 // =======================================================================
 177 template<class T, int N>
 178 Standard_Integer BVH_LinearBuilder<T, N>::EmitHierachy (BVH_Tree<T, N>*                        theBVH,
 179                                                         const Standard_Integer                 theBit,
 180                                                         const Standard_Integer                 theShift,
 181                                                         std::vector<BVH_EncodedLink>::iterator theStart,
 182                                                         std::vector<BVH_EncodedLink>::iterator theFinal)
 183 {
 184   if (theFinal - theStart > BVH_Builder<T, N>::myLeafNodeSize)
 185   {
 186     std::vector<BVH_EncodedLink>::iterator aPosition =
 187       (theBit >= 0) ? std::lower_bound (theStart, theFinal, BVH_EncodedLink(), BVH::BitComparator (theBit))
 188                     :                   theStart + ((theFinal - theStart) / 2);
 189
 190     if (aPosition == theStart || aPosition == theFinal)
 191     {
 192       return EmitHierachy (theBVH, theBit - 1, theShift, theStart, theFinal);
 193     }
 194
 195     // Build inner node
 196     const Standard_Integer aNode = theBVH->AddInnerNode (0, 0);
 197
 198     const Standard_Integer aRghNode = theShift + static_cast<Standard_Integer> (aPosition - theStart);
 199
 200     const Standard_Integer aLftChild = EmitHierachy (theBVH, theBit - 1, theShift, theStart, aPosition);
 201     const Standard_Integer aRghChild = EmitHierachy (theBVH, theBit - 1, aRghNode, aPosition, theFinal);
 202
 203     theBVH->NodeInfoBuffer()[aNode].y() = aLftChild;
 204     theBVH->NodeInfoBuffer()[aNode].z() = aRghChild;
 205
 206     return aNode;
 207   }
 208   else
 209   {
 210     // Build leaf node
 211     return theBVH->AddLeafNode (theShift, theShift + static_cast<Standard_Integer> (theFinal - theStart) - 1);
 212   }
 213 }
 214
 215 #ifdef HAVE_TBB
 216
 217 namespace BVH
 218 {
 219   //! TBB task for parallel radix sort.
 220   class RadixSortTask : public tbb::task
 221   {
 222     typedef std::vector<BVH_EncodedLink>::iterator LinkIterator;
 223
 224   private:
 225
 226     //! Start range element.
 227     LinkIterator myStart;
 228
 229     //! Final range element.
 230     LinkIterator myFinal;
 231
 232     //! Bit position for range partition.
 233     Standard_Integer myDigit;
 234
 235   public:
 236
 237     //! Creates new TBB radix sort task.
 238     RadixSortTask (LinkIterator theStart, LinkIterator theFinal, Standard_Integer theDigit)
 239     : myStart (theStart),
 240       myFinal (theFinal),
 241       myDigit (theDigit)
 242     {
 243       //
 244     }
 245
 246     //! Executes the task.
 247     tbb::task* execute()
 248     {
 249       if (myDigit < 28)
 250       {
 251         BVH::RadixSorter::Perform (myStart, myFinal, myDigit);
 252       }
 253       else
 254       {
 255         LinkIterator anOffset = std::partition (myStart, myFinal, BitPredicate (myDigit));
 256
 257         tbb::task_list aList;
 258
 259         aList.push_back (*new ( allocate_child() )
 260           RadixSortTask (myStart, anOffset, myDigit - 1));
 261
 262         aList.push_back (*new ( allocate_child() )
 263           RadixSortTask (anOffset, myFinal, myDigit - 1));
 264
 265         set_ref_count (3); // count + 1
 266         spawn_and_wait_for_all (aList);
 267       }
 268
 269       return NULL;
 270     }
 271   };
 272
 273   //! TBB task for parallel bounds updating.
 274   template<class T, int N>
 275   class UpdateBoundTask: public tbb::task
 276   {
 277     //! Set of geometric objects.
 278     BVH_Set<T, N>* mySet;
 279
 280     //! BVH tree built over the set.
 281     BVH_Tree<T, N>* myBVH;
 282
 283     //! BVH node to update bounding box.
 284     Standard_Integer myNode;
 285
 286     //! Level of the processed BVH node.
 287     Standard_Integer myLevel;
 288
 289     //! Height of the processed BVH node.
 290     Standard_Integer* myHeight;
 291
 292   public:
 293
 294     //! Creates new TBB parallel bound update task.
 295     UpdateBoundTask (BVH_Set<T, N>*    theSet,
 296                      BVH_Tree<T, N>*   theBVH,
 297                      Standard_Integer  theNode,
 298                      Standard_Integer  theLevel,
 299                      Standard_Integer* theHeight)
 300     : mySet    (theSet),
 301       myBVH    (theBVH),
 302       myNode   (theNode),
 303       myLevel  (theLevel),
 304       myHeight (theHeight)
 305     {
 306       //
 307     }
 308
 309     //! Executes the task.
 310     tbb::task* execute()
 311     {
 312       if (myBVH->IsOuter (myNode) || myLevel > 2)
 313       {
 314         *myHeight = BVH::UpdateBounds (mySet, myBVH, myNode);
 315       }
 316       else
 317       {
 318         Standard_Integer aLftHeight = 0;
 319         Standard_Integer aRghHeight = 0;
 320
 321         tbb::task_list aList;
 322
 323         const Standard_Integer aLftChild = myBVH->NodeInfoBuffer()[myNode].y();
 324         const Standard_Integer aRghChild = myBVH->NodeInfoBuffer()[myNode].z();
 325
 326         Standard_Integer aCount = 1;
 327
 328         if (!myBVH->IsOuter (aLftChild))
 329         {
 330           ++aCount;
 331           aList.push_back (*new ( allocate_child() )
 332             UpdateBoundTask (mySet, myBVH, aLftChild, myLevel + 1, &aLftHeight));
 333         }
 334         else
 335         {
 336           aLftHeight = BVH::UpdateBounds (mySet, myBVH, aLftChild);
 337         }
 338
 339         if (!myBVH->IsOuter (aRghChild))
 340         {
 341           ++aCount;
 342           aList.push_back (*new( allocate_child() )
 343             UpdateBoundTask (mySet, myBVH, aRghChild, myLevel + 1, &aRghHeight));
 344         }
 345         else
 346         {
 347           aRghHeight = BVH::UpdateBounds (mySet, myBVH, aRghChild);
 348         }
 349
 350         if (aCount > 1)
 351         {
 352           set_ref_count (aCount);
 353           spawn_and_wait_for_all (aList);
 354         }
 355
 356         typename BVH_Box<T, N>::BVH_VecNt aLftMinPoint = myBVH->MinPointBuffer()[aLftChild];
 357         typename BVH_Box<T, N>::BVH_VecNt aLftMaxPoint = myBVH->MaxPointBuffer()[aLftChild];
 358         typename BVH_Box<T, N>::BVH_VecNt aRghMinPoint = myBVH->MinPointBuffer()[aRghChild];
 359         typename BVH_Box<T, N>::BVH_VecNt aRghMaxPoint = myBVH->MaxPointBuffer()[aRghChild];
 360
 361         BVH::BoxMinMax<T, N>::CwiseMin (aLftMinPoint, aRghMinPoint);
 362         BVH::BoxMinMax<T, N>::CwiseMax (aLftMaxPoint, aRghMaxPoint);
 363
 364         myBVH->MinPointBuffer()[myNode] = aLftMinPoint;
 365         myBVH->MaxPointBuffer()[myNode] = aLftMaxPoint;
 366
 367         *myHeight = Max (aLftHeight, aRghHeight) + 1;
 368       }
 369
 370       return NULL;
 371     }
 372   };
 373 }
 374
 375 #endif
 376
 377 // =======================================================================
 378 // function : Build
 379 // purpose  :
 380 // =======================================================================
 381 template<class T, int N>
 382 void BVH_LinearBuilder<T, N>::Build (BVH_Set<T, N>*       theSet,
 383                                      BVH_Tree<T, N>*      theBVH,
 384                                      const BVH_Box<T, N>& theBox)
 385 {
 386   Standard_STATIC_ASSERT (N == 3 || N == 4);
 387   const Standard_Integer aSetSize = theSet->Size();
 388   if (theBVH == NULL || aSetSize == 0)
 389   {
 390     return;
 391   }
 392
 393   theBVH->Clear();
 394
 395   const Standard_Integer aDimensionX = 1024;
 396   const Standard_Integer aDimensionY = 1024;
 397   const Standard_Integer aDimensionZ = 1024;
 398
 399   const BVH_VecNt aSceneMin = theBox.CornerMin();
 400   const BVH_VecNt aSceneMax = theBox.CornerMax();
 401
 402   const T aMinSize = static_cast<T> (BVH::THE_NODE_MIN_SIZE);
 403
 404   const T aReverseSizeX = static_cast<T> (aDimensionX) / Max (aMinSize, aSceneMax.x() - aSceneMin.x());
 405   const T aReverseSizeY = static_cast<T> (aDimensionY) / Max (aMinSize, aSceneMax.y() - aSceneMin.y());
 406   const T aReverseSizeZ = static_cast<T> (aDimensionZ) / Max (aMinSize, aSceneMax.z() - aSceneMin.z());
 407
 408   std::vector<BVH_EncodedLink> anEncodedLinks (aSetSize, BVH_EncodedLink());
 409
 410   // Step 1 -- Assign Morton code to each primitive
 411   for (Standard_Integer aPrimIdx = 0; aPrimIdx < aSetSize; ++aPrimIdx)
 412   {
 413     const BVH_VecNt aCenter = theSet->Box (aPrimIdx).Center();
 414
 415     Standard_Integer aVoxelX = BVH::IntFloor ((aCenter.x() - aSceneMin.x()) * aReverseSizeX);
 416     Standard_Integer aVoxelY = BVH::IntFloor ((aCenter.y() - aSceneMin.y()) * aReverseSizeY);
 417     Standard_Integer aVoxelZ = BVH::IntFloor ((aCenter.z() - aSceneMin.z()) * aReverseSizeZ);
 418
 419     aVoxelX = Max (0, Min (aVoxelX, aDimensionX - 1));
 420     aVoxelY = Max (0, Min (aVoxelY, aDimensionY - 1));
 421     aVoxelZ = Max (0, Min (aVoxelZ, aDimensionZ - 1));
 422
 423     aVoxelX = (aVoxelX | (aVoxelX << 16)) & 0x030000FF;
 424     aVoxelX = (aVoxelX | (aVoxelX <<  8)) & 0x0300F00F;
 425     aVoxelX = (aVoxelX | (aVoxelX <<  4)) & 0x030C30C3;
 426     aVoxelX = (aVoxelX | (aVoxelX <<  2)) & 0x09249249;
 427
 428     aVoxelY = (aVoxelY | (aVoxelY << 16)) & 0x030000FF;
 429     aVoxelY = (aVoxelY | (aVoxelY <<  8)) & 0x0300F00F;
 430     aVoxelY = (aVoxelY | (aVoxelY <<  4)) & 0x030C30C3;
 431     aVoxelY = (aVoxelY | (aVoxelY <<  2)) & 0x09249249;
 432
 433     aVoxelZ = (aVoxelZ | (aVoxelZ << 16)) & 0x030000FF;
 434     aVoxelZ = (aVoxelZ | (aVoxelZ <<  8)) & 0x0300F00F;
 435     aVoxelZ = (aVoxelZ | (aVoxelZ <<  4)) & 0x030C30C3;
 436     aVoxelZ = (aVoxelZ | (aVoxelZ <<  2)) & 0x09249249;
 437
 438     anEncodedLinks[aPrimIdx] = BVH_EncodedLink (
 439       aVoxelX | (aVoxelY << 1) | (aVoxelZ << 2), aPrimIdx);
 440   }
 441
 442   // Step 2 -- Sort primitives by their Morton codes using radix sort
 443 #ifdef HAVE_TBB
 444
 445   BVH::RadixSortTask& aSortTask = *new ( tbb::task::allocate_root() )
 446     BVH::RadixSortTask (anEncodedLinks.begin(), anEncodedLinks.end(), 29);
 447
 448   tbb::task::spawn_root_and_wait (aSortTask);
 449
 450 #else
 451
 452   BVH::RadixSorter::Perform (anEncodedLinks.begin(), anEncodedLinks.end());
 453
 454 #endif
 455
 456   // Step 3 -- Emitting BVH hierarchy from sorted Morton codes
 457   EmitHierachy (theBVH, 29, 0, anEncodedLinks.begin(), anEncodedLinks.end());
 458
 459   NCollection_Array1<Standard_Integer> aLinkMap (0, aSetSize - 1);
 460   for (Standard_Integer aLinkIdx = 0; aLinkIdx < aSetSize; ++aLinkIdx)
 461   {
 462     aLinkMap (anEncodedLinks[aLinkIdx].second) = aLinkIdx;
 463   }
 464
 465   // Step 4 -- Rearranging primitive list according to Morton codes (in place)
 466   Standard_Integer aPrimIdx = 0;
 467
 468   while (aPrimIdx < aSetSize)
 469   {
 470     const Standard_Integer aSortIdx = aLinkMap (aPrimIdx);
 471
 472     if (aPrimIdx != aSortIdx)
 473     {
 474       theSet->Swap (aPrimIdx, aSortIdx);
 475
 476       std::swap (aLinkMap (aPrimIdx),
 477                  aLinkMap (aSortIdx));
 478     }
 479     else
 480     {
 481       ++aPrimIdx;
 482     }
 483   }
 484
 485   // Step 5 -- Compute bounding boxes of BVH nodes
 486   theBVH->MinPointBuffer().resize (theBVH->NodeInfoBuffer().size());
 487   theBVH->MaxPointBuffer().resize (theBVH->NodeInfoBuffer().size());
 488
 489   Standard_Integer aDepth = 0;
 490
 491 #ifdef HAVE_TBB
 492
 493   // Note: Although TBB tasks are allocated using placement
 494   // new, we do not need to delete them explicitly
 495   BVH::UpdateBoundTask<T, N>& aRootTask = *new ( tbb::task::allocate_root() )
 496     BVH::UpdateBoundTask<T, N> (theSet, theBVH, 0, 0, &aDepth);
 497
 498   tbb::task::spawn_root_and_wait (aRootTask);
 499
 500 #else
 501
 502   aDepth = BVH::UpdateBounds (theSet, theBVH, 0);
 503
 504 #endif
 505
 506   BVH_Builder<T, N>::UpdateDepth (theBVH, aDepth);
 507 }