1 /*
2 Copyright (c) 2011-2017 Timur Gafarov 
3 
4 Boost Software License - Version 1.0 - August 17th, 2003
5 
6 Permission is hereby granted, free of charge, to any person or organization
7 obtaining a copy of the software and accompanying documentation covered by
8 this license (the "Software") to use, reproduce, display, distribute,
9 execute, and transmit the Software, and to prepare derivative works of the
10 Software, and to permit third-parties to whom the Software is furnished to
11 do so, all subject to the following:
12 
13 The copyright notices in the Software and this entire statement, including
14 the above license grant, this restriction and the following disclaimer,
15 must be included in all copies of the Software, in whole or in part, and
16 all derivative works of the Software, unless such copies or derivative
17 works are solely in the form of machine-executable object code generated by
18 a source language processor.
19 
20 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
21 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
23 SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
24 FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
25 ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 DEALINGS IN THE SOFTWARE.
27 */
28 
29 module dmech.bvh;
30 
31 import std.array;
32 import std.math;
33 
34 import dlib.core.memory;
35 import dlib.core.compound;
36 import dlib.container.array;
37 import dlib.math.utils;
38 import dlib.math.vector;
39 import dlib.geometry.aabb;
40 import dlib.geometry.sphere;
41 import dlib.geometry.ray;
42 
43 /*
44  * Bounding Volume Hierarchy implementation
45  */
46 
47 // Returns the axis that has the largest length
48 Axis boxGetMainAxis(AABB box)
49 {
50     float xl = box.size.x;
51     float yl = box.size.y;
52     float zl = box.size.z;
53          
54     if (xl < yl)
55     {
56         if (yl < zl)
57            return Axis.z;
58         return Axis.y;
59     }
60     else if (xl < zl)
61         return Axis.z;
62     return Axis.x;        
63 }
64 
65 struct SplitPlane
66 {
67     public:
68     float split;
69     Axis axis;
70     
71     this(float s, Axis ax)
72     {
73         split = s;
74         axis = ax;
75     }
76 }
77 
78 SplitPlane boxGetSplitPlaneForAxis(AABB box, Axis a)
79 {
80     return SplitPlane(box.center[a], a);
81 }
82 
83 Compound!(AABB, AABB) boxSplitWithPlane(AABB box, SplitPlane sp)
84 {
85     Vector3f minLP = box.pmin;
86     Vector3f maxLP = box.pmax;
87     maxLP[sp.axis] = sp.split;
88     
89     Vector3f minRP = box.pmin;
90     Vector3f maxRP = box.pmax;
91     minRP[sp.axis] = sp.split;
92 
93     AABB leftB = boxFromMinMaxPoints(minLP, maxLP);
94     AABB rightB = boxFromMinMaxPoints(minRP, maxRP);
95 
96     return compound(leftB, rightB);
97 }
98 
99 AABB enclosingAABB(T)(T[] objects)
100 {
101     Vector3f pmin = objects[0].boundingBox.pmin;
102     Vector3f pmax = pmin;
103     
104     void adjustMinPoint(Vector3f p)
105     {    
106         if (p.x < pmin.x) pmin.x = p.x;
107         if (p.y < pmin.y) pmin.y = p.y;
108         if (p.z < pmin.z) pmin.z = p.z;
109     }
110     
111     void adjustMaxPoint(Vector3f p)
112     {
113         if (p.x > pmax.x) pmax.x = p.x;
114         if (p.y > pmax.y) pmax.y = p.y;
115         if (p.z > pmax.z) pmax.z = p.z;
116     }
117 
118     foreach(ref obj; objects)
119     {
120         adjustMinPoint(obj.boundingBox.pmin);
121         adjustMaxPoint(obj.boundingBox.pmax);
122     }
123     
124     return boxFromMinMaxPoints(pmin, pmax);
125 }
126 
127 class BVHNode(T)
128 {
129     DynamicArray!T objects;
130     AABB aabb;
131     BVHNode[2] child;
132     uint userData;
133 
134     this(DynamicArray!T objs)
135     {
136         objects = objs;
137         aabb = enclosingAABB(objects.data);
138     }
139     
140     void free()
141     {
142         objects.free();
143         if (child[0] !is null) child[0].free();
144         if (child[1] !is null) child[1].free();
145         Delete(this);
146     }
147     
148     SphereTraverseAggregate!T traverseBySphere(Sphere* sphere)
149     {
150         return SphereTraverseAggregate!(T)(this, sphere);
151     }
152     
153     RayTraverseAggregate!T traverseByRay(Ray* ray)
154     {
155         return RayTraverseAggregate!(T)(this, ray);
156     }
157 }
158 
159 struct SphereTraverseAggregate(T)
160 {
161     BVHNode!T node;
162     Sphere* sphere;
163     
164     int opApply(int delegate(ref T) dg)
165     {
166         int result = 0;
167         
168         Vector3f cn;
169         float pd;
170         if (node.aabb.intersectsSphere(*sphere, cn, pd))
171         {        
172             if (node.child[0] !is null)
173             {
174                 result = node.child[0].traverseBySphere(sphere).opApply(dg);
175                 if (result)
176                     return result;
177             }
178             
179             if (node.child[1] !is null)
180             {
181                 result = node.child[1].traverseBySphere(sphere).opApply(dg);
182                 if (result)
183                     return result;
184             }
185             
186             foreach(ref obj; node.objects.data)
187                 dg(obj);
188         }
189         else
190             return result;
191             
192         return result;
193     }
194 }
195 
196 struct RayTraverseAggregate(T)
197 {
198     BVHNode!T node;
199     Ray* ray;
200     
201     int opApply(int delegate(ref T) dg) // TODO: nearest intersection point
202     {
203         int result = 0;
204         
205         float it = 0.0f;
206         if (node.aabb.intersectsSegment(ray.p0, ray.p1, it))    
207         { 
208             if (node.child[0] !is null)
209             {
210                 result = node.child[0].traverseByRay(ray).opApply(dg);
211                 if (result)
212                     return result;
213             }
214             
215             if (node.child[1] !is null)
216             {
217                 result = node.child[1].traverseByRay(ray).opApply(dg);
218                 if (result)
219                     return result;
220             }
221             
222             foreach(ref obj; node.objects.data)
223                 dg(obj);
224         }
225         else
226             return result;
227             
228         return result;
229     }
230 }
231 
232 /+
233 void traverseBySphere(T)(BVHNode!T node, ref Sphere sphere /*, void delegate(ref T) func*/)
234 {
235     Vector3f cn;
236     float pd;
237     if (node.aabb.intersectsSphere(sphere, cn, pd))
238     {
239         //if (node.child[0] !is null)
240         //    node.child[0].traverseBySphere(sphere, func);
241         //if (node.child[1] !is null)
242         //    node.child[1].traverseBySphere(sphere, func);
243 
244         //foreach(ref obj; node.objects.data)
245         //    func(obj);
246     }
247 }
248 +/
249 /*
250 void traverse(T)(BVHNode!T node, void delegate(BVHNode!T) func)
251 {
252     if (node.child[0] !is null)
253         node.child[0].traverse(func);
254     if (node.child[1] !is null)
255         node.child[1].traverse(func);
256 
257     func(node);
258 }
259 */
260 /*
261 void traverseByRay(T)(BVHNode!T node, Ray ray, void delegate(ref T) func)
262 {
263     float it = 0.0f;
264     if (node.aabb.intersectsSegment(ray.p0, ray.p1, it))
265     {
266         if (node.child[0] !is null)
267             node.child[0].traverseByRay(ray, func);
268         if (node.child[1] !is null)
269             node.child[1].traverseByRay(ray, func);
270 
271         foreach(ref obj; node.objects.data)
272             func(obj);
273     }
274 }
275 */
276 
277 // TODO:
278 // - support multithreading (2 children = 2 threads)
279 // - add ESC (Early Split Clipping)
280 enum Heuristic
281 {
282     HMA, // Half Main Axis
283     SAH, // Surface Area Heuristic
284     //ESC  // Early Split Clipping
285 }
286 
287 DynamicArray!T duplicate(T)(DynamicArray!T arr)
288 {
289     DynamicArray!T res;
290     foreach(v; arr.data)
291         res.append(v);
292     return res;
293 }
294 
295 class BVHTree(T)
296 {
297     BVHNode!T root;
298 
299     this(DynamicArray!T objects, 
300          uint maxObjectsPerNode = 8,
301          uint maxRecursionDepth = 10,
302          Heuristic splitHeuristic = Heuristic.SAH)
303     {
304         root = construct(objects, 0, maxObjectsPerNode, maxRecursionDepth, splitHeuristic);
305     }
306     
307     void free()
308     {
309         root.free();
310         Delete(this);
311     }
312     
313     import std.stdio;
314 
315     BVHNode!T construct(
316          DynamicArray!T objects, 
317          uint rec,
318          uint maxObjectsPerNode,
319          uint maxRecursionDepth,
320          Heuristic splitHeuristic)
321     {
322         BVHNode!T node = New!(BVHNode!T)(duplicate(objects));
323 
324         if (node.objects.data.length <= maxObjectsPerNode)
325         {
326             return node;
327         }
328         
329         if (rec == maxRecursionDepth)
330         {
331             return node;
332         }
333         
334         AABB box = enclosingAABB(node.objects.data);
335 
336         SplitPlane sp;
337         if (splitHeuristic == Heuristic.HMA)
338             sp = getHalfMainAxisSplitPlane(node.objects.data, box);
339         else if (splitHeuristic == Heuristic.SAH)
340             sp = getSAHSplitPlane(node.objects.data, box);
341         else
342             assert(0, "BVH: unsupported split heuristic");
343 
344         auto boxes = boxSplitWithPlane(box, sp);
345 
346         DynamicArray!T leftObjects;
347         DynamicArray!T rightObjects;
348         
349         foreach(obj; node.objects.data)
350         {
351             if (boxes[0].intersectsAABB(obj.boundingBox))
352                 leftObjects.append(obj);
353             else if (boxes[1].intersectsAABB(obj.boundingBox))
354                 rightObjects.append(obj);
355         }
356 
357         if (leftObjects.data.length > 0 || rightObjects.data.length > 0)
358             node.objects.free();
359 
360         if (leftObjects.data.length > 0)
361             node.child[0] = construct(leftObjects, rec + 1, maxObjectsPerNode, maxRecursionDepth, splitHeuristic);
362         else
363             node.child[0] = null;
364 
365         if (rightObjects.data.length > 0)
366             node.child[1] = construct(rightObjects, rec + 1, maxObjectsPerNode, maxRecursionDepth, splitHeuristic);
367         else
368             node.child[1] = null;
369 
370         leftObjects.free();
371         rightObjects.free();
372 
373         return node;    
374     }
375 
376     SplitPlane getHalfMainAxisSplitPlane(T[] objects, ref AABB box)
377     {
378         Axis axis = boxGetMainAxis(box);
379         return boxGetSplitPlaneForAxis(box, axis);
380     }
381 
382     SplitPlane getSAHSplitPlane(T[] objects, ref AABB box)
383     {
384         Axis axis = boxGetMainAxis(box);
385         
386         float minAlongSplitPlane = box.pmin[axis];
387         float maxAlongSplitPlane = box.pmax[axis];
388         
389         float bestSAHCost = float.nan;
390         float bestSplitPoint = float.nan;
391 
392         int iterations = 12;
393         foreach (i; 0..iterations)
394         {
395             float valueOfSplit = minAlongSplitPlane + 
396                                ((maxAlongSplitPlane - minAlongSplitPlane) / (iterations + 1.0f) * (i + 1.0f));
397 
398             SplitPlane SAHSplitPlane = SplitPlane(valueOfSplit, axis);
399             auto boxes = boxSplitWithPlane(box, SAHSplitPlane);
400 
401             uint leftObjectsLength = 0;
402             uint rightObjectsLength = 0;
403 
404             foreach(obj; objects)
405             {
406                 if (boxes[0].intersectsAABB(obj.boundingBox))
407                     leftObjectsLength++;
408                 else if (boxes[1].intersectsAABB(obj.boundingBox))
409                     rightObjectsLength++;
410             }
411 
412             if (leftObjectsLength > 0 && rightObjectsLength > 0)
413             {
414                 float SAHCost = getSAHCost(boxes[0], leftObjectsLength, 
415                                            boxes[1], rightObjectsLength, box);
416 
417                 if (bestSAHCost.isNaN || SAHCost < bestSAHCost)
418                 {
419                     bestSAHCost = SAHCost;
420                     bestSplitPoint = valueOfSplit;
421                 }
422             }
423         }
424         
425         return SplitPlane(bestSplitPoint, axis);
426     }
427 
428     float getSAHCost(AABB leftBox, uint numLeftObjects, 
429                      AABB rightBox, uint numRightObjects,
430                      AABB parentBox)
431     {
432         return getSurfaceArea(leftBox) / getSurfaceArea(parentBox) * numLeftObjects
433              + getSurfaceArea(rightBox) / getSurfaceArea(parentBox) * numRightObjects;
434     }
435 
436     float getSurfaceArea(AABB bbox)
437     {
438         float width = bbox.pmax.x - bbox.pmin.x;
439         float height = bbox.pmax.y - bbox.pmin.y;
440         float depth = bbox.pmax.z - bbox.pmin.z;
441         return 2.0f * (width * height + width * depth + height * depth);
442     }
443 }
444