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 && theBit >= 0)
 185   {
 186     std::vector<BVH_EncodedLink>::iterator aPosition = std::lower_bound (
 187       theStart, theFinal, BVH_EncodedLink(), BVH::BitComparator (theBit));
 188
 189     if (aPosition == theStart || aPosition == theFinal)
 190     {
 191       return EmitHierachy (theBVH, theBit - 1, theShift, theStart, theFinal);
 192     }
 193
 194     // Build inner node
 195     const Standard_Integer aNode = theBVH->AddInnerNode (0, 0);
 196
 197     const Standard_Integer aRghNode = theShift + static_cast<Standard_Integer> (aPosition - theStart);
 198
 199     const Standard_Integer aLftChild = EmitHierachy (theBVH, theBit - 1, theShift, theStart, aPosition);
 200     const Standard_Integer aRghChild = EmitHierachy (theBVH, theBit - 1, aRghNode, aPosition, theFinal);
 201
 202     theBVH->NodeInfoBuffer()[aNode].y() = aLftChild;
 203     theBVH->NodeInfoBuffer()[aNode].z() = aRghChild;
 204
 205     return aNode;
 206   }
 207   else
 208   {
 209     // Build leaf node
 210     return theBVH->AddLeafNode (theShift, theShift + static_cast<Standard_Integer> (theFinal - theStart) - 1);
 211   }
 212 }
 213
 214 #ifdef HAVE_TBB
 215
 216 namespace BVH
 217 {
 218   //! TBB task for parallel radix sort.
 219   class RadixSortTask : public tbb::task
 220   {
 221     typedef std::vector<BVH_EncodedLink>::iterator LinkIterator;
 222
 223   private:
 224
 225     //! Start range element.
 226     LinkIterator myStart;
 227
 228     //! Final range element.
 229     LinkIterator myFinal;
 230
 231     //! Bit position for range partition.
 232     Standard_Integer myDigit;
 233
 234   public:
 235
 236     //! Creates new TBB radix sort task.
 237     RadixSortTask (LinkIterator theStart, LinkIterator theFinal, Standard_Integer theDigit)
 238     : myStart (theStart),
 239       myFinal (theFinal),
 240       myDigit (theDigit)
 241     {
 242       //
 243     }
 244
 245     //! Executes the task.
 246     tbb::task* execute()
 247     {
 248       if (myDigit < 28)
 249       {
 250         BVH::RadixSorter::Perform (myStart, myFinal, myDigit);
 251       }
 252       else
 253       {
 254         LinkIterator anOffset = std::partition (myStart, myFinal, BitPredicate (myDigit));
 255
 256         tbb::task_list aList;
 257
 258         aList.push_back (*new ( allocate_child() )
 259           RadixSortTask (myStart, anOffset, myDigit - 1));
 260
 261         aList.push_back (*new ( allocate_child() )
 262           RadixSortTask (anOffset, myFinal, myDigit - 1));
 263
 264         set_ref_count (3); // count + 1
 265         spawn_and_wait_for_all (aList);
 266       }
 267
 268       return NULL;
 269     }
 270   };
 271
 272   //! TBB task for parallel bounds updating.
 273   template<class T, int N>
 274   class UpdateBoundTask: public tbb::task
 275   {
 276     //! Set of geometric objects.
 277     BVH_Set<T, N>* mySet;
 278
 279     //! BVH tree built over the set.
 280     BVH_Tree<T, N>* myBVH;
 281
 282     //! BVH node to update bounding box.
 283     Standard_Integer myNode;
 284
 285     //! Level of the processed BVH node.
 286     Standard_Integer myLevel;
 287
 288     //! Height of the processed BVH node.
 289     Standard_Integer* myHeight;
 290
 291   public:
 292
 293     //! Creates new TBB parallel bound update task.
 294     UpdateBoundTask (BVH_Set<T, N>*    theSet,
 295                      BVH_Tree<T, N>*   theBVH,
 296                      Standard_Integer  theNode,
 297                      Standard_Integer  theLevel,
 298                      Standard_Integer* theHeight)
 299     : mySet    (theSet),
 300       myBVH    (theBVH),
 301       myNode   (theNode),
 302       myLevel  (theLevel),
 303       myHeight (theHeight)
 304     {
 305       //
 306     }
 307
 308     //! Executes the task.
 309     tbb::task* execute()
 310     {
 311       if (myBVH->IsOuter (myNode) || myLevel > 2)
 312       {
 313         *myHeight = BVH::UpdateBounds (mySet, myBVH, myNode);
 314       }
 315       else
 316       {
 317         Standard_Integer aLftHeight = 0;
 318         Standard_Integer aRghHeight = 0;
 319
 320         tbb::task_list aList;
 321
 322         const Standard_Integer aLftChild = myBVH->NodeInfoBuffer()[myNode].y();
 323         const Standard_Integer aRghChild = myBVH->NodeInfoBuffer()[myNode].z();
 324
 325         Standard_Integer aCount = 1;
 326
 327         if (!myBVH->IsOuter (aLftChild))
 328         {
 329           ++aCount;
 330           aList.push_back (*new ( allocate_child() )
 331             UpdateBoundTask (mySet, myBVH, aLftChild, myLevel + 1, &aLftHeight));
 332         }
 333         else
 334         {
 335           aLftHeight = BVH::UpdateBounds (mySet, myBVH, aLftChild);
 336         }
 337
 338         if (!myBVH->IsOuter (aRghChild))
 339         {
 340           ++aCount;
 341           aList.push_back (*new( allocate_child() )
 342             UpdateBoundTask (mySet, myBVH, aRghChild, myLevel + 1, &aRghHeight));
 343         }
 344         else
 345         {
 346           aRghHeight = BVH::UpdateBounds (mySet, myBVH, aRghChild);
 347         }
 348
 349         if (aCount > 1)
 350         {
 351           set_ref_count (aCount);
 352           spawn_and_wait_for_all (aList);
 353         }
 354
 355         typename BVH_Box<T, N>::BVH_VecNt aLftMinPoint = myBVH->MinPointBuffer()[aLftChild];
 356         typename BVH_Box<T, N>::BVH_VecNt aLftMaxPoint = myBVH->MaxPointBuffer()[aLftChild];
 357         typename BVH_Box<T, N>::BVH_VecNt aRghMinPoint = myBVH->MinPointBuffer()[aRghChild];
 358         typename BVH_Box<T, N>::BVH_VecNt aRghMaxPoint = myBVH->MaxPointBuffer()[aRghChild];
 359
 360         BVH::BoxMinMax<T, N>::CwiseMin (aLftMinPoint, aRghMinPoint);
 361         BVH::BoxMinMax<T, N>::CwiseMax (aLftMaxPoint, aRghMaxPoint);
 362
 363         myBVH->MinPointBuffer()[myNode] = aLftMinPoint;
 364         myBVH->MaxPointBuffer()[myNode] = aLftMaxPoint;
 365
 366         *myHeight = Max (aLftHeight, aRghHeight) + 1;
 367       }
 368
 369       return NULL;
 370     }
 371   };
 372 }
 373
 374 #endif
 375
 376 // =======================================================================
 377 // function : Build
 378 // purpose  :
 379 // =======================================================================
 380 template<class T, int N>
 381 void BVH_LinearBuilder<T, N>::Build (BVH_Set<T, N>*       theSet,
 382                                      BVH_Tree<T, N>*      theBVH,
 383                                      const BVH_Box<T, N>& theBox)
 384 {
 385   Standard_STATIC_ASSERT (N == 3 || N == 4);
 386
 387   if (theBVH == NULL || theSet->Size() == 0)
 388   {
 389     return;
 390   }
 391
 392   theBVH->Clear();
 393
 394   const Standard_Integer aDimensionX = 1024;
 395   const Standard_Integer aDimensionY = 1024;
 396   const Standard_Integer aDimensionZ = 1024;
 397
 398   const BVH_VecNt aSceneMin = theBox.CornerMin();
 399   const BVH_VecNt aSceneMax = theBox.CornerMax();
 400
 401   const T aMinSize = static_cast<T> (BVH::THE_NODE_MIN_SIZE);
 402
 403   const T aReverseSizeX = static_cast<T> (aDimensionX) / Max (aMinSize, aSceneMax.x() - aSceneMin.x());
 404   const T aReverseSizeY = static_cast<T> (aDimensionY) / Max (aMinSize, aSceneMax.y() - aSceneMin.y());
 405   const T aReverseSizeZ = static_cast<T> (aDimensionZ) / Max (aMinSize, aSceneMax.z() - aSceneMin.z());
 406
 407   std::vector<BVH_EncodedLink> anEncodedLinks (theSet->Size(), BVH_EncodedLink());
 408
 409   // Step 1 -- Assign Morton code to each primitive
 410   for (Standard_Integer aPrimIdx = 0; aPrimIdx < theSet->Size(); ++aPrimIdx)
 411   {
 412     const BVH_VecNt aCenter = theSet->Box (aPrimIdx).Center();
 413
 414     Standard_Integer aVoxelX = BVH::IntFloor ((aCenter.x() - aSceneMin.x()) * aReverseSizeX);
 415     Standard_Integer aVoxelY = BVH::IntFloor ((aCenter.y() - aSceneMin.y()) * aReverseSizeY);
 416     Standard_Integer aVoxelZ = BVH::IntFloor ((aCenter.z() - aSceneMin.z()) * aReverseSizeZ);
 417
 418     aVoxelX = Max (0, Min (aVoxelX, aDimensionX - 1));
 419     aVoxelY = Max (0, Min (aVoxelY, aDimensionY - 1));
 420     aVoxelZ = Max (0, Min (aVoxelZ, aDimensionZ - 1));
 421
 422     aVoxelX = (aVoxelX | (aVoxelX << 16)) & 0x030000FF;
 423     aVoxelX = (aVoxelX | (aVoxelX <<  8)) & 0x0300F00F;
 424     aVoxelX = (aVoxelX | (aVoxelX <<  4)) & 0x030C30C3;
 425     aVoxelX = (aVoxelX | (aVoxelX <<  2)) & 0x09249249;
 426
 427     aVoxelY = (aVoxelY | (aVoxelY << 16)) & 0x030000FF;
 428     aVoxelY = (aVoxelY | (aVoxelY <<  8)) & 0x0300F00F;
 429     aVoxelY = (aVoxelY | (aVoxelY <<  4)) & 0x030C30C3;
 430     aVoxelY = (aVoxelY | (aVoxelY <<  2)) & 0x09249249;
 431
 432     aVoxelZ = (aVoxelZ | (aVoxelZ << 16)) & 0x030000FF;
 433     aVoxelZ = (aVoxelZ | (aVoxelZ <<  8)) & 0x0300F00F;
 434     aVoxelZ = (aVoxelZ | (aVoxelZ <<  4)) & 0x030C30C3;
 435     aVoxelZ = (aVoxelZ | (aVoxelZ <<  2)) & 0x09249249;
 436
 437     anEncodedLinks[aPrimIdx] = BVH_EncodedLink (
 438       aVoxelX | (aVoxelY << 1) | (aVoxelZ << 2), aPrimIdx);
 439   }
 440
 441   // Step 2 -- Sort primitives by their Morton codes using radix sort
 442 #ifdef HAVE_TBB
 443
 444   BVH::RadixSortTask& aSortTask = *new ( tbb::task::allocate_root() )
 445     BVH::RadixSortTask (anEncodedLinks.begin(), anEncodedLinks.end(), 29);
 446
 447   tbb::task::spawn_root_and_wait (aSortTask);
 448
 449 #else
 450
 451   BVH::RadixSorter::Perform (anEncodedLinks.begin(), anEncodedLinks.end());
 452
 453 #endif
 454
 455   // Step 3 -- Emitting BVH hierarchy from sorted Morton codes
 456   EmitHierachy (theBVH, 29, 0, anEncodedLinks.begin(), anEncodedLinks.end());
 457
 458   NCollection_Array1<Standard_Integer> aLinkMap (0, theSet->Size() - 1);
 459
 460   for (Standard_Integer aLinkIdx = 0; aLinkIdx < theSet->Size(); ++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 < theSet->Size())
 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 }