0026364: Foundation Classes, TKMath - Optimize BVH binned algorithm
[occt.git] / src / BVH / BVH_BinnedBuilder.lxx
1 // Created on: 2013-12-20
2 // Created by: Denis BOGOLEPOV
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 // =======================================================================
17 // function : BVH_BinnedBuilder
18 // purpose  :
19 // =======================================================================
20 template<class T, int N, int Bins>
21 BVH_BinnedBuilder<T, N, Bins>::BVH_BinnedBuilder (const Standard_Integer theLeafNodeSize,
22                                                   const Standard_Integer theMaxTreeDepth,
23                                                   const Standard_Boolean theDoMainSplits,
24                                                   const Standard_Integer theNumOfThreads)
25 : BVH_QueueBuilder<T, N> (theLeafNodeSize,
26                           theMaxTreeDepth,
27                           theNumOfThreads),
28   myUseMainAxis (theDoMainSplits)
29 {
30   //
31 }
32
33 // =======================================================================
34 // function : ~BVH_BinnedBuilder
35 // purpose  :
36 // =======================================================================
37 template<class T, int N, int Bins>
38 BVH_BinnedBuilder<T, N, Bins>::~BVH_BinnedBuilder()
39 {
40   //
41 }
42
43 // =======================================================================
44 // function : GetSubVolumes
45 // purpose  :
46 // =======================================================================
47 template<class T, int N, int Bins>
48 void BVH_BinnedBuilder<T, N, Bins>::GetSubVolumes (BVH_Set<T, N>*         theSet,
49                                                    BVH_Tree<T, N>*        theBVH,
50                                                    const Standard_Integer theNode,
51                                                    BVH_BinVector&         theBins,
52                                                    const Standard_Integer theAxis)
53 {
54   const T aMin = BVH::VecComp<T, N>::Get (theBVH->MinPoint (theNode), theAxis);
55   const T aMax = BVH::VecComp<T, N>::Get (theBVH->MaxPoint (theNode), theAxis);
56
57   const T anInverseStep = static_cast<T> (Bins) / (aMax - aMin);
58
59   for (Standard_Integer anIdx = theBVH->BegPrimitive (theNode); anIdx <= theBVH->EndPrimitive (theNode); ++anIdx)
60   {
61     typename BVH_Set<T, N>::BVH_BoxNt aBox = theSet->Box (anIdx);
62
63     Standard_Integer aBinIndex = BVH::IntFloor<T> (
64       (theSet->Center (anIdx, theAxis) - aMin) * anInverseStep);
65
66     if (aBinIndex < 0)
67     {
68       aBinIndex = 0;
69     }
70     else if (aBinIndex >= Bins)
71     {
72       aBinIndex = Bins - 1;
73     }
74
75     theBins[aBinIndex].Count++;
76     theBins[aBinIndex].Box.Combine (aBox);
77   }
78 }
79
80 namespace BVH
81 {
82   // =======================================================================
83   // function : SplitPrimitives
84   // purpose  :
85   // =======================================================================
86   template<class T, int N>
87   Standard_Integer SplitPrimitives (BVH_Set<T, N>*         theSet,
88                                     const BVH_Box<T, N>&   theBox,
89                                     const Standard_Integer theBeg,
90                                     const Standard_Integer theEnd,
91                                     const Standard_Integer theBin,
92                                     const Standard_Integer theAxis,
93                                     const Standard_Integer theBins)
94   {
95     const T aMin = BVH::VecComp<T, N>::Get (theBox.CornerMin(), theAxis);
96     const T aMax = BVH::VecComp<T, N>::Get (theBox.CornerMax(), theAxis);
97
98     const T anInverseStep = static_cast<T> (theBins) / (aMax - aMin);
99
100     Standard_Integer aLftIdx (theBeg);
101     Standard_Integer aRghIdx (theEnd);
102
103     do
104     {
105       while (BVH::IntFloor<T> ((theSet->Center (aLftIdx, theAxis) - aMin) * anInverseStep) <= theBin && aLftIdx < theEnd)
106       {
107         ++aLftIdx;
108       }
109       while (BVH::IntFloor<T> ((theSet->Center (aRghIdx, theAxis) - aMin) * anInverseStep) > theBin && aRghIdx > theBeg)
110       {
111         --aRghIdx;
112       }
113
114       if (aLftIdx <= aRghIdx)
115       {
116         if (aLftIdx != aRghIdx)
117         {
118           theSet->Swap (aLftIdx, aRghIdx);
119         }
120
121         ++aLftIdx;
122         --aRghIdx;
123       }
124     } while (aLftIdx <= aRghIdx);
125
126     return aLftIdx;
127   }
128 }
129
130 #if defined (_WIN32) && defined (max)
131   #undef max
132 #endif
133
134 #include <limits>
135
136 namespace BVH
137 {
138   template<class T, int N>
139   struct BVH_AxisSelector
140   {
141     typedef typename BVH::VectorType<T, N>::Type BVH_VecNt;
142
143     // =======================================================================
144     // function : MainAxis
145     // purpose  :
146     // =======================================================================
147     static Standard_Integer MainAxis (const BVH_VecNt& theSize)
148     {
149       if (theSize.y() > theSize.x())
150       {
151         return theSize.y() > theSize.z() ? 1 : 2;
152       }
153       else
154       {
155         return theSize.z() > theSize.x() ? 2 : 0;
156       }
157     }
158   };
159
160   template<class T>
161   struct BVH_AxisSelector<T, 2>
162   {
163     typedef typename BVH::VectorType<T, 2>::Type BVH_VecNt;
164
165     // =======================================================================
166     // function : MainAxis
167     // purpose  :
168     // =======================================================================
169     static Standard_Integer MainAxis (const BVH_VecNt& theSize)
170     {
171       return theSize.x() > theSize.y() ? 0 : 1;
172     }
173   };
174 }
175
176 // =======================================================================
177 // function : BuildNode
178 // purpose  :
179 // =======================================================================
180 template<class T, int N, int Bins>
181 typename BVH_QueueBuilder<T, N>::BVH_ChildNodes BVH_BinnedBuilder<T, N, Bins>::BuildNode (BVH_Set<T, N>*         theSet,
182                                                                                           BVH_Tree<T, N>*        theBVH,
183                                                                                           const Standard_Integer theNode)
184 {
185   const Standard_Integer aNodeBegPrimitive = theBVH->BegPrimitive (theNode);
186   const Standard_Integer aNodeEndPrimitive = theBVH->EndPrimitive (theNode);
187
188   if (aNodeEndPrimitive - aNodeBegPrimitive < BVH_Builder<T, N>::myLeafNodeSize)
189   {
190     return typename BVH_QueueBuilder<T, N>::BVH_ChildNodes(); // node does not require partitioning
191   }
192
193   const BVH_Box<T, N> anAABB (theBVH->MinPoint (theNode),
194                               theBVH->MaxPoint (theNode));
195
196   const typename BVH_Box<T, N>::BVH_VecNt aSize = anAABB.Size();
197
198   // Parameters for storing best split
199   Standard_Integer aMinSplitAxis   = -1;
200   Standard_Integer aMinSplitIndex  =  0;
201   Standard_Integer aMinSplitNumLft =  0;
202   Standard_Integer aMinSplitNumRgh =  0;
203
204   BVH_Box<T, N> aMinSplitBoxLft;
205   BVH_Box<T, N> aMinSplitBoxRgh;
206
207   Standard_Real aMinSplitCost = std::numeric_limits<Standard_Real>::max();
208
209   const Standard_Integer aMainAxis = BVH::BVH_AxisSelector<T, N>::MainAxis (aSize);
210
211   // Find best split
212   for (Standard_Integer anAxis = myUseMainAxis ? aMainAxis : 0; anAxis <= (myUseMainAxis ? aMainAxis : Min (N - 1, 2)); ++anAxis)
213   {
214     if (BVH::VecComp<T, N>::Get (aSize, anAxis) <= BVH::THE_NODE_MIN_SIZE)
215     {
216       continue;
217     }
218
219     BVH_BinVector aBinVector;
220     GetSubVolumes (theSet, theBVH, theNode, aBinVector, anAxis);
221
222     BVH_SplitPlanes aSplitPlanes;
223     for (Standard_Integer aLftSplit = 1, aRghSplit = Bins - 1; aLftSplit < Bins; ++aLftSplit, --aRghSplit)
224     {
225       aSplitPlanes[aLftSplit].LftVoxel.Count = aSplitPlanes[aLftSplit - 1].LftVoxel.Count + aBinVector[aLftSplit - 1].Count;
226       aSplitPlanes[aRghSplit].RghVoxel.Count = aSplitPlanes[aRghSplit + 1].RghVoxel.Count + aBinVector[aRghSplit + 0].Count;
227
228       aSplitPlanes[aLftSplit].LftVoxel.Box = aSplitPlanes[aLftSplit - 1].LftVoxel.Box;
229       aSplitPlanes[aRghSplit].RghVoxel.Box = aSplitPlanes[aRghSplit + 1].RghVoxel.Box;
230
231       aSplitPlanes[aLftSplit].LftVoxel.Box.Combine (aBinVector[aLftSplit - 1].Box);
232       aSplitPlanes[aRghSplit].RghVoxel.Box.Combine (aBinVector[aRghSplit + 0].Box);
233     }
234
235     // Choose the best split (with minimum SAH cost)
236     for (Standard_Integer aSplit = 1; aSplit < Bins; ++aSplit)
237     {
238       // Simple SAH evaluation
239       Standard_Real aCost =
240         (static_cast<Standard_Real> (aSplitPlanes[aSplit].LftVoxel.Box.Area()) /* / S(N) */) * aSplitPlanes[aSplit].LftVoxel.Count
241       + (static_cast<Standard_Real> (aSplitPlanes[aSplit].RghVoxel.Box.Area()) /* / S(N) */) * aSplitPlanes[aSplit].RghVoxel.Count;
242
243       if (aCost <= aMinSplitCost)
244       {
245         aMinSplitCost   = aCost;
246         aMinSplitAxis   = anAxis;
247         aMinSplitIndex  = aSplit;
248         aMinSplitBoxLft = aSplitPlanes[aSplit].LftVoxel.Box;
249         aMinSplitBoxRgh = aSplitPlanes[aSplit].RghVoxel.Box;
250         aMinSplitNumLft = aSplitPlanes[aSplit].LftVoxel.Count;
251         aMinSplitNumRgh = aSplitPlanes[aSplit].RghVoxel.Count;
252       }
253     }
254   }
255
256   theBVH->SetInner (theNode);
257
258   Standard_Integer aMiddle = -1;
259
260   if (aMinSplitNumLft == 0 || aMinSplitNumRgh == 0 || aMinSplitAxis == -1) // case of objects with the same center
261   {
262     aMinSplitBoxLft.Clear();
263     aMinSplitBoxRgh.Clear();
264
265     aMiddle = std::max (aNodeBegPrimitive + 1,
266       static_cast<Standard_Integer> ((aNodeBegPrimitive + aNodeEndPrimitive) / 2.f));
267
268     aMinSplitNumLft = aMiddle - aNodeBegPrimitive;
269
270     for (Standard_Integer anIndex = aNodeBegPrimitive; anIndex < aMiddle; ++anIndex)
271     {
272       aMinSplitBoxLft.Combine (theSet->Box (anIndex));
273     }
274
275     aMinSplitNumRgh = aNodeEndPrimitive - aMiddle + 1;
276
277     for (Standard_Integer anIndex = aNodeEndPrimitive; anIndex >= aMiddle; --anIndex)
278     {
279       aMinSplitBoxRgh.Combine (theSet->Box (anIndex));
280     }
281   }
282   else
283   {
284     aMiddle = BVH::SplitPrimitives<T, N> (theSet,
285                                           anAABB,
286                                           aNodeBegPrimitive,
287                                           aNodeEndPrimitive,
288                                           aMinSplitIndex - 1,
289                                           aMinSplitAxis,
290                                           Bins);
291   }
292
293   typedef typename BVH_QueueBuilder<T, N>::BVH_PrimitiveRange Range;
294
295   return typename BVH_QueueBuilder<T, N>::BVH_ChildNodes (aMinSplitBoxLft,
296                                                           aMinSplitBoxRgh,
297                                                           Range (aNodeBegPrimitive, aMiddle - 1),
298                                                           Range (aMiddle,     aNodeEndPrimitive));
299 }