dmech.bvh source code

1 /*
2 Copyright (c) 2011-2014 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.compound;
35 import dlib.math.utils;
36 import dlib.math.vector;
37 import dlib.geometry.aabb;
38 import dlib.geometry.sphere;
39 import dlib.geometry.ray;
40 
41 /*
42  * Bounding Volume Hierarchy implementation
43  */
44 
45 // Returns the axis that has the largest length
46 Axis boxGetMainAxis(AABB box)
47 {
48     float xl = box.size.x;
49     float yl = box.size.y;
50     float zl = box.size.z;
51          
52     if (xl < yl)
53     {
54         if (yl < zl)
55            return Axis.z;
56         return Axis.y;
57     }
58     else if (xl < zl)
59         return Axis.z;
60     return Axis.x;        
61 }
62 
63 struct SplitPlane
64 {
65     public:
66     float split;
67     Axis axis;
68     
69     this(float s, Axis ax)
70     {
71         split = s;
72         axis = ax;
73     }
74 }
75 
76 SplitPlane boxGetSplitPlaneForAxis(AABB box, Axis a)
77 {
78     return SplitPlane(box.center[a], a);
79 }
80 
81 Compound!(AABB, AABB) boxSplitWithPlane(AABB box, SplitPlane sp)
82 {
83     Vector3f minLP = box.pmin;
84     Vector3f maxLP = box.pmax;
85     maxLP[sp.axis] = sp.split;
86     
87     Vector3f minRP = box.pmin;
88     Vector3f maxRP = box.pmax;
89     minRP[sp.axis] = sp.split;
90 
91     AABB leftB = boxFromMinMaxPoints(minLP, maxLP);
92     AABB rightB = boxFromMinMaxPoints(minRP, maxRP);
93 
94     return compound(leftB, rightB);
95 }
96 
97 AABB enclosingAABB(T)(T[] objects)
98 {
99     Vector3f pmin = objects[0].boundingBox.pmin;
100     Vector3f pmax = pmin;
101     
102     void adjustMinPoint(Vector3f p)
103     {    
104         if (p.x < pmin.x) pmin.x = p.x;
105         if (p.y < pmin.y) pmin.y = p.y;
106         if (p.z < pmin.z) pmin.z = p.z;
107     }
108     
109     void adjustMaxPoint(Vector3f p)
110     {
111         if (p.x > pmax.x) pmax.x = p.x;
112         if (p.y > pmax.y) pmax.y = p.y;
113         if (p.z > pmax.z) pmax.z = p.z;
114     }
115 
116     foreach(ref obj; objects)
117     {
118         adjustMinPoint(obj.boundingBox.pmin);
119         adjustMaxPoint(obj.boundingBox.pmax);
120     }
121     
122     return boxFromMinMaxPoints(pmin, pmax);
123 }
124 
125 class BVHNode(T)
126 {
127     T[] objects;
128     AABB aabb;
129     BVHNode[2] child;
130     uint userData;
131 
132     this(T[] objs)
133     {
134         objects = objs;
135         aabb = enclosingAABB(objects);
136     }
137 }
138 
139 void traverseBySphere(T)(BVHNode!T node, ref Sphere sphere, void delegate(ref T) func)
140 {
141     Vector3f cn;
142     float pd;
143     if (node.aabb.intersectsSphere(sphere, cn, pd))
144     {
145         if (node.child[0] !is null)
146             node.child[0].traverseBySphere(sphere, func);
147         if (node.child[1] !is null)
148             node.child[1].traverseBySphere(sphere, func);
149 
150         foreach(ref obj; node.objects)
151             func(obj);
152     }
153 }
154 
155 void traverse(T)(BVHNode!T node, void delegate(BVHNode!T) func)
156 {
157     if (node.child[0] !is null)
158         node.child[0].traverse(func);
159     if (node.child[1] !is null)
160         node.child[1].traverse(func);
161 
162     func(node);
163 }
164 
165 void traverseByRay(T)(BVHNode!T node, Ray ray, void delegate(ref T) func)
166 {
167     float it = 0.0f;
168     if (node.aabb.intersectsSegment(ray.p0, ray.p1, it))
169     {
170         if (node.child[0] !is null)
171             node.child[0].traverseByRay(ray, func);
172         if (node.child[1] !is null)
173             node.child[1].traverseByRay(ray, func);
174 
175         foreach(ref obj; node.objects)
176             func(obj);
177     }
178 }
179 
180 // TODO:
181 // - support multithreading (2 children = 2 threads)
182 // - add ESC (Early Split Clipping)
183 enum Heuristic
184 {
185     HMA, // Half Main Axis
186     SAH, // Surface Area Heuristic
187     //ESC  // Early Split Clipping
188 }
189 
190 class BVHTree(T)
191 {
192     BVHNode!T root;
193 
194     this(T[] objects, 
195          uint maxObjectsPerNode = 8,
196          Heuristic splitHeuristic = Heuristic.SAH)
197     {
198         root = construct(objects, maxObjectsPerNode, splitHeuristic);
199     }
200 
201     BVHNode!T construct(
202          T[] objects, 
203          uint maxObjectsPerNode,
204          Heuristic splitHeuristic)
205     {
206         AABB box = enclosingAABB(objects);
207         
208         SplitPlane sp;
209         if (splitHeuristic == Heuristic.HMA)
210             sp = getHalfMainAxisSplitPlane(objects, box);
211         else if (splitHeuristic == Heuristic.SAH)
212             sp = getSAHSplitPlane(objects, box);
213         else
214             assert(0, "BVH: unsupported split heuristic");
215             
216         auto boxes = boxSplitWithPlane(box, sp);
217 
218         T[] leftObjects;
219         T[] rightObjects;
220     
221         foreach(obj; objects)
222         {
223             if (boxes[0].intersectsAABB(obj.boundingBox))
224                 leftObjects ~= obj;
225             else if (boxes[1].intersectsAABB(obj.boundingBox))
226                 rightObjects ~= obj;
227         }
228     
229         BVHNode!T node = new BVHNode!T(objects);
230 
231         if (objects.length <= maxObjectsPerNode)
232             return node;
233         
234         if (leftObjects.length > 0 || rightObjects.length > 0)
235             node.objects = [];
236 
237         if (leftObjects.length > 0)
238             node.child[0] = construct(leftObjects, maxObjectsPerNode, splitHeuristic);
239         else
240             node.child[0] = null;
241     
242         if (rightObjects.length > 0)
243             node.child[1] = construct(rightObjects, maxObjectsPerNode, splitHeuristic);
244         else
245             node.child[1] = null;
246 
247         return node;    
248     }
249 
250     SplitPlane getHalfMainAxisSplitPlane(ref T[] objects, ref AABB box)
251     {
252         Axis axis = boxGetMainAxis(box);
253         return boxGetSplitPlaneForAxis(box, axis);
254     }
255 
256     SplitPlane getSAHSplitPlane(ref T[] objects, ref AABB box)
257     {
258         Axis axis = boxGetMainAxis(box);
259         
260         float minAlongSplitPlane = box.pmin[axis];
261         float maxAlongSplitPlane = box.pmax[axis];
262         
263         float bestSAHCost = float.nan;
264         float bestSplitPoint = float.nan;
265 
266         int iterations = 12;
267         foreach (i; 0..iterations)
268         {
269             float valueOfSplit = minAlongSplitPlane + 
270                                ((maxAlongSplitPlane - minAlongSplitPlane) / (iterations + 1.0f) * (i + 1.0f));
271 
272             SplitPlane SAHSplitPlane = SplitPlane(valueOfSplit, axis);
273             auto boxes = boxSplitWithPlane(box, SAHSplitPlane);
274 
275             uint leftObjectsLength = 0;
276             uint rightObjectsLength = 0;
277 
278             foreach(obj; objects)
279             {
280                 if (boxes[0].intersectsAABB(obj.boundingBox))
281                     leftObjectsLength++;
282                 else if (boxes[1].intersectsAABB(obj.boundingBox))
283                     rightObjectsLength++;
284             }
285 
286             if (leftObjectsLength > 0 && rightObjectsLength > 0)
287             {
288                 float SAHCost = getSAHCost(boxes[0], leftObjectsLength, 
289                                            boxes[1], rightObjectsLength, box);
290 
291                 if (bestSAHCost.isNaN || SAHCost < bestSAHCost)
292                 {
293                     bestSAHCost = SAHCost;
294                     bestSplitPoint = valueOfSplit;
295                 }
296             }
297         }
298         
299         return SplitPlane(bestSplitPoint, axis);
300     }
301 
302     float getSAHCost(AABB leftBox, uint numLeftObjects, 
303                      AABB rightBox, uint numRightObjects,
304                      AABB parentBox)
305     {
306         return getSurfaceArea(leftBox) / getSurfaceArea(parentBox) * numLeftObjects
307              + getSurfaceArea(rightBox) / getSurfaceArea(parentBox) * numRightObjects;
308     }
309 
310     float getSurfaceArea(AABB bbox)
311     {
312         float width = bbox.pmax.x - bbox.pmin.x;
313         float height = bbox.pmax.y - bbox.pmin.y;
314         float depth = bbox.pmax.z - bbox.pmin.z;
315         return 2.0f * (width * height + width * depth + height * depth);
316     }
317 }
318