0030476: Visualization, Path Tracing - Adaptive Screen Sampling leads to unstable...
[occt.git] / src / BVH / BVH_RadixSorter.hxx
1 // Created on: 2016-04-13
2 // Created by: Denis BOGOLEPOV
3 // Copyright (c) 2013-2016 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 #ifndef _BVH_RadixSorter_Header
17 #define _BVH_RadixSorter_Header
18
19 #include <BVH_Sorter.hxx>
20 #include <BVH_Builder.hxx>
21 #include <NCollection_Array1.hxx>
22 #include <NCollection_Shared.hxx>
23 #include <OSD_Parallel.hxx>
24
25 #include <algorithm>
26
27 //! Pair of Morton code and primitive ID.
28 typedef std::pair<Standard_Integer, Standard_Integer> BVH_EncodedLink;
29
30 //! Performs radix sort of a BVH primitive set using
31 //! 10-bit Morton codes (or 1024 x 1024 x 1024 grid).
32 template<class T, int N>
33 class BVH_RadixSorter : public BVH_Sorter<T, N>
34 {
35 public:
36
37   typedef typename BVH::VectorType<T, N>::Type BVH_VecNt;
38
39 public:
40
41   //! Creates new BVH radix sorter for the given AABB.
42   BVH_RadixSorter (const BVH_Box<T, N>& theBox) : myBox (theBox) { }
43
44   //! Sorts the set.
45   virtual void Perform (BVH_Set<T, N>* theSet) Standard_OVERRIDE { Perform (theSet, 0, theSet->Size() - 1); }
46
47   //! Sorts the given (inclusive) range in the set.
48   virtual void Perform (BVH_Set<T, N>* theSet, const Standard_Integer theStart, const Standard_Integer theFinal) Standard_OVERRIDE;
49
50   //! Returns Morton codes assigned to BVH primitives.
51   const NCollection_Array1<BVH_EncodedLink>& EncodedLinks() const { return *myEncodedLinks; }
52
53 protected:
54
55   //! Axis-aligned bounding box (AABB) to perform sorting.
56   BVH_Box<T, N> myBox;
57
58   //! Morton codes assigned to BVH primitives.
59   Handle(NCollection_Shared<NCollection_Array1<BVH_EncodedLink> >) myEncodedLinks;
60
61 };
62
63 namespace BVH
64 {
65   // Radix sort STL predicate for 32-bit integer.
66   struct BitPredicate
67   {
68     Standard_Integer myBit;
69
70     //! Creates new radix sort predicate.
71     BitPredicate (const Standard_Integer theBit) : myBit (theBit) {}
72
73     //! Returns predicate value.
74     bool operator() (const BVH_EncodedLink theLink) const
75     {
76       const Standard_Integer aMask = 1 << myBit;
77       return !(theLink.first & aMask); // 0-bit to the left side
78     }
79   };
80
81   //! STL compare tool used in binary search algorithm.
82   struct BitComparator
83   {
84     Standard_Integer myBit;
85
86     //! Creates new STL comparator.
87     BitComparator (const Standard_Integer theBit) : myBit (theBit) {}
88
89     //! Checks left value for the given bit.
90     bool operator() (BVH_EncodedLink theLink1, BVH_EncodedLink /*theLink2*/)
91     {
92       return !(theLink1.first & (1 << myBit));
93     }
94   };
95
96   //! Tool object for sorting link array using radix sort algorithm.
97   class RadixSorter
98   {
99   public:
100
101     typedef NCollection_Array1<BVH_EncodedLink>::iterator LinkIterator;
102
103   private:
104
105     //! Structure defining sorting range.
106     struct SortRange
107     {
108       LinkIterator     myStart; //!< Start element of exclusive sorting range
109       LinkIterator     myFinal; //!< Final element of exclusive sorting range
110       Standard_Integer myDigit; //!< Bit number used for partition operation
111     };
112
113     //! Functor class to run sorting in parallel.
114     class Functor
115     {
116     public:
117       Functor(const SortRange (&aSplits)[2], const Standard_Boolean isParallel)
118       : mySplits     (aSplits),
119         myIsParallel (isParallel)
120       {
121       }
122       
123       //! Runs sorting function for the given range.
124       void operator()(const Standard_Integer theIndex) const
125       {
126         RadixSorter::Sort (mySplits[theIndex].myStart, mySplits[theIndex].myFinal,
127                            mySplits[theIndex].myDigit, myIsParallel);
128       }
129
130     private:
131       void operator=(const Functor&);
132       
133     private:
134       const SortRange (&mySplits)[2];
135       Standard_Boolean myIsParallel;
136     };
137
138   public:
139
140     static void Sort (LinkIterator theStart, LinkIterator theFinal, Standard_Integer theDigit, const Standard_Boolean isParallel)
141     {
142       if (theDigit < 24)
143       {
144         BVH::RadixSorter::perform (theStart, theFinal, theDigit);
145       }
146       else
147       {
148         LinkIterator anOffset = std::partition (theStart, theFinal, BitPredicate (theDigit));
149         SortRange aSplits[2] = {
150           {theStart, anOffset, theDigit - 1},
151           {anOffset, theFinal, theDigit - 1}
152         };
153
154         OSD_Parallel::For (0, 2, Functor (aSplits, isParallel), !isParallel);
155       }
156     }
157
158   protected:
159
160     // Performs MSD (most significant digit) radix sort.
161     static void perform (LinkIterator theStart, LinkIterator theFinal, Standard_Integer theBit = 29)
162     {
163       while (theStart != theFinal && theBit >= 0)
164       {
165         LinkIterator anOffset = std::partition (theStart, theFinal, BitPredicate (theBit--));
166         perform (theStart, anOffset, theBit);
167         theStart = anOffset;
168       }
169     }
170   };
171 }
172
173 // =======================================================================
174 // function : Perform
175 // purpose  :
176 // =======================================================================
177 template<class T, int N>
178 void BVH_RadixSorter<T, N>::Perform (BVH_Set<T, N>* theSet, const Standard_Integer theStart, const Standard_Integer theFinal)
179 {
180   Standard_STATIC_ASSERT (N == 2 || N == 3 || N == 4);
181
182   const Standard_Integer aDimension = 1024;
183   const Standard_Integer aNbEffComp = N == 2 ? 2 : 3; // 4th component is ignored
184
185   const BVH_VecNt aSceneMin = myBox.CornerMin();
186   const BVH_VecNt aSceneMax = myBox.CornerMax();
187
188   BVH_VecNt aNodeMinSizeVecT (static_cast<T>(BVH::THE_NODE_MIN_SIZE));
189   BVH::BoxMinMax<T, N>::CwiseMax (aNodeMinSizeVecT, aSceneMax - aSceneMin);
190
191   const BVH_VecNt aReverseSize = BVH_VecNt (static_cast<T>(aDimension)) / aNodeMinSizeVecT;
192
193   myEncodedLinks = new NCollection_Shared<NCollection_Array1<BVH_EncodedLink> >(theStart, theFinal);
194
195   // Step 1 -- Assign Morton code to each primitive
196   for (Standard_Integer aPrimIdx = theStart; aPrimIdx <= theFinal; ++aPrimIdx)
197   {
198     const BVH_VecNt aCenter = theSet->Box (aPrimIdx).Center();
199     const BVH_VecNt aVoxelF = (aCenter - aSceneMin) * aReverseSize;
200
201     Standard_Integer aMortonCode = 0;
202     for (Standard_Integer aCompIter = 0; aCompIter < aNbEffComp; ++aCompIter)
203     {
204       Standard_Integer aVoxel = BVH::IntFloor (BVH::VecComp<T, N>::Get (aVoxelF, aCompIter));
205
206       aVoxel = Max (0, Min (aVoxel, aDimension - 1));
207
208       aVoxel = (aVoxel | (aVoxel << 16)) & 0x030000FF;
209       aVoxel = (aVoxel | (aVoxel <<  8)) & 0x0300F00F;
210       aVoxel = (aVoxel | (aVoxel <<  4)) & 0x030C30C3;
211       aVoxel = (aVoxel | (aVoxel <<  2)) & 0x09249249;
212
213       aMortonCode |= (aVoxel << aCompIter);
214     }
215
216     myEncodedLinks->ChangeValue (aPrimIdx) = BVH_EncodedLink (aMortonCode, aPrimIdx);
217   }
218
219   // Step 2 -- Sort primitives by their Morton codes using radix sort
220   BVH::RadixSorter::Sort (myEncodedLinks->begin(), myEncodedLinks->end(), 29, this->IsParallel());
221
222   NCollection_Array1<Standard_Integer> aLinkMap (theStart, theFinal);
223   for (Standard_Integer aLinkIdx = theStart; aLinkIdx <= theFinal; ++aLinkIdx)
224   {
225     aLinkMap (myEncodedLinks->Value (aLinkIdx).second) = aLinkIdx;
226   }
227
228   // Step 3 -- Rearranging primitive list according to Morton codes (in place)
229   Standard_Integer aPrimIdx = theStart;
230   while (aPrimIdx <= theFinal)
231   {
232     const Standard_Integer aSortIdx = aLinkMap (aPrimIdx);
233     if (aPrimIdx != aSortIdx)
234     {
235       theSet->Swap (aPrimIdx, aSortIdx);
236       std::swap (aLinkMap (aPrimIdx),
237                  aLinkMap (aSortIdx));
238     }
239     else
240     {
241       ++aPrimIdx;
242     }
243   }
244 }
245
246 #endif // _BVH_RadixSorter_Header