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