Cleanup: Remove more #if 0 blocks
[blender.git] / source / blender / blenlib / intern / BLI_kdopbvh.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2006 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): Daniel Genrich, Andre Pinto
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/blenlib/intern/BLI_kdopbvh.c
29  *  \ingroup bli
30  *  \brief BVH-tree implementation.
31  *
32  * k-DOP BVH (Discrete Oriented Polytope, Bounding Volume Hierarchy).
33  * A k-DOP is represented as k/2 pairs of min , max values for k/2 directions (intervals, "slabs").
34  *
35  * See: http://www.gris.uni-tuebingen.de/people/staff/jmezger/papers/bvh.pdf
36  *
37  * implements a bvh-tree structure with support for:
38  *
39  * - Ray-cast:
40  *   #BLI_bvhtree_ray_cast, #BVHRayCastData
41  * - Nearest point on surface:
42  *   #BLI_bvhtree_find_nearest, #BVHNearestData
43  * - Overlapping 2 trees:
44  *   #BLI_bvhtree_overlap, #BVHOverlapData_Shared, #BVHOverlapData_Thread
45  * - Range Query:
46  *   #BLI_bvhtree_range_query
47  */
48
49 #include <assert.h>
50
51 #include "MEM_guardedalloc.h"
52
53 #include "BLI_utildefines.h"
54 #include "BLI_alloca.h"
55 #include "BLI_stack.h"
56 #include "BLI_kdopbvh.h"
57 #include "BLI_math.h"
58 #include "BLI_task.h"
59
60 #include "BLI_strict_flags.h"
61
62 /* used for iterative_raycast */
63 // #define USE_SKIP_LINKS
64
65 /* Use to print balanced output. */
66 // #define USE_PRINT_TREE
67
68 /* Check tree is valid. */
69 // #define USE_VERIFY_TREE
70
71
72 #define MAX_TREETYPE 32
73
74 /* Setting zero so we can catch bugs in BLI_task/KDOPBVH.
75  * TODO(sergey): Deduplicate the limits with PBVH from BKE.
76  */
77 #ifdef DEBUG
78 #  define KDOPBVH_THREAD_LEAF_THRESHOLD 0
79 #else
80 #  define KDOPBVH_THREAD_LEAF_THRESHOLD 1024
81 #endif
82
83
84 /* -------------------------------------------------------------------- */
85
86 /** \name Struct Definitions
87  * \{ */
88
89 typedef unsigned char axis_t;
90
91 typedef struct BVHNode {
92         struct BVHNode **children;
93         struct BVHNode *parent; /* some user defined traversed need that */
94 #ifdef USE_SKIP_LINKS
95         struct BVHNode *skip[2];
96 #endif
97         float *bv;      /* Bounding volume of all nodes, max 13 axis */
98         int index;      /* face, edge, vertex index */
99         char totnode;   /* how many nodes are used, used for speedup */
100         char main_axis; /* Axis used to split this node */
101 } BVHNode;
102
103 /* keep under 26 bytes for speed purposes */
104 struct BVHTree {
105         BVHNode **nodes;
106         BVHNode *nodearray;     /* pre-alloc branch nodes */
107         BVHNode **nodechild;    /* pre-alloc childs for nodes */
108         float   *nodebv;        /* pre-alloc bounding-volumes for nodes */
109         float epsilon;          /* epslion is used for inflation of the k-dop      */
110         int totleaf;            /* leafs */
111         int totbranch;
112         axis_t start_axis, stop_axis;  /* bvhtree_kdop_axes array indices according to axis */
113         axis_t axis;                   /* kdop type (6 => OBB, 7 => AABB, ...) */
114         char tree_type;                /* type of tree (4 => quadtree) */
115 };
116
117 /* optimization, ensure we stay small */
118 BLI_STATIC_ASSERT((sizeof(void *) == 8 && sizeof(BVHTree) <= 48) ||
119                   (sizeof(void *) == 4 && sizeof(BVHTree) <= 32),
120                   "over sized")
121
122 /* avoid duplicating vars in BVHOverlapData_Thread */
123 typedef struct BVHOverlapData_Shared {
124         const BVHTree *tree1, *tree2;
125         axis_t start_axis, stop_axis;
126
127         /* use for callbacks */
128         BVHTree_OverlapCallback callback;
129         void *userdata;
130 } BVHOverlapData_Shared;
131
132 typedef struct BVHOverlapData_Thread {
133         BVHOverlapData_Shared *shared;
134         struct BLI_Stack *overlap;  /* store BVHTreeOverlap */
135         /* use for callbacks */
136         int thread;
137 } BVHOverlapData_Thread;
138
139 typedef struct BVHNearestData {
140         const BVHTree *tree;
141         const float *co;
142         BVHTree_NearestPointCallback callback;
143         void    *userdata;
144         float proj[13];         /* coordinates projection over axis */
145         BVHTreeNearest nearest;
146
147 } BVHNearestData;
148
149 typedef struct BVHRayCastData {
150         const BVHTree *tree;
151
152         BVHTree_RayCastCallback callback;
153         void    *userdata;
154
155
156         BVHTreeRay ray;
157
158 #ifdef USE_KDOPBVH_WATERTIGHT
159         struct IsectRayPrecalc isect_precalc;
160 #endif
161
162         /* initialized by bvhtree_ray_cast_data_precalc */
163         float ray_dot_axis[13];
164         float idot_axis[13];
165         int index[6];
166
167         BVHTreeRayHit hit;
168 } BVHRayCastData;
169
170 typedef struct BVHNearestProjectedData {
171         const BVHTree *tree;
172         struct DistProjectedAABBPrecalc precalc;
173         bool closest_axis[3];
174         float clip_plane[6][4];
175         int clip_plane_len;
176         BVHTree_NearestProjectedCallback callback;
177         void *userdata;
178         BVHTreeNearest nearest;
179
180 } BVHNearestProjectedData;
181
182 /** \} */
183
184
185 /**
186  * Bounding Volume Hierarchy Definition
187  *
188  * Notes: From OBB until 26-DOP --> all bounding volumes possible, just choose type below
189  * Notes: You have to choose the type at compile time ITM
190  * Notes: You can choose the tree type --> binary, quad, octree, choose below
191  */
192
193 const float bvhtree_kdop_axes[13][3] = {
194         {1.0, 0, 0}, {0, 1.0, 0}, {0, 0, 1.0},
195         {1.0, 1.0, 1.0}, {1.0, -1.0, 1.0}, {1.0, 1.0, -1.0}, {1.0, -1.0, -1.0},
196         {1.0, 1.0, 0}, {1.0, 0, 1.0}, {0, 1.0, 1.0}, {1.0, -1.0, 0}, {1.0, 0, -1.0}, {0, 1.0, -1.0}
197 };
198
199
200 /* -------------------------------------------------------------------- */
201
202 /** \name Utility Functions
203  * \{ */
204
205 MINLINE axis_t min_axis(axis_t a, axis_t b)
206 {
207         return (a < b) ? a : b;
208 }
209 #if 0
210 MINLINE axis_t max_axis(axis_t a, axis_t b)
211 {
212         return (b < a) ? a : b;
213 }
214 #endif
215
216 /**
217  * Introsort
218  * with permission deriven from the following Java code:
219  * http://ralphunden.net/content/tutorials/a-guide-to-introsort/
220  * and he derived it from the SUN STL
221  */
222
223
224
225 static void node_minmax_init(const BVHTree *tree, BVHNode *node)
226 {
227         axis_t axis_iter;
228         float (*bv)[2] = (float (*)[2])node->bv;
229
230         for (axis_iter = tree->start_axis; axis_iter != tree->stop_axis; axis_iter++) {
231                 bv[axis_iter][0] =  FLT_MAX;
232                 bv[axis_iter][1] = -FLT_MAX;
233         }
234 }
235
236 /** \} */
237
238
239 /* -------------------------------------------------------------------- */
240
241 /** \name Balance Utility Functions
242  * \{ */
243
244 /**
245  * Insertion sort algorithm
246  */
247 static void bvh_insertionsort(BVHNode **a, int lo, int hi, int axis)
248 {
249         int i, j;
250         BVHNode *t;
251         for (i = lo; i < hi; i++) {
252                 j = i;
253                 t = a[i];
254                 while ((j != lo) && (t->bv[axis] < (a[j - 1])->bv[axis])) {
255                         a[j] = a[j - 1];
256                         j--;
257                 }
258                 a[j] = t;
259         }
260 }
261
262 static int bvh_partition(BVHNode **a, int lo, int hi, BVHNode *x, int axis)
263 {
264         int i = lo, j = hi;
265         while (1) {
266                 while (a[i]->bv[axis] < x->bv[axis]) {
267                         i++;
268                 }
269                 j--;
270                 while (x->bv[axis] < a[j]->bv[axis]) {
271                         j--;
272                 }
273                 if (!(i < j)) {
274                         return i;
275                 }
276                 SWAP(BVHNode *, a[i], a[j]);
277                 i++;
278         }
279 }
280
281 static BVHNode *bvh_medianof3(BVHNode **a, int lo, int mid, int hi, int axis)  /* returns Sortable */
282 {
283         if ((a[mid])->bv[axis] < (a[lo])->bv[axis]) {
284                 if ((a[hi])->bv[axis] < (a[mid])->bv[axis])
285                         return a[mid];
286                 else {
287                         if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
288                                 return a[hi];
289                         else
290                                 return a[lo];
291                 }
292         }
293         else {
294                 if ((a[hi])->bv[axis] < (a[mid])->bv[axis]) {
295                         if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
296                                 return a[lo];
297                         else
298                                 return a[hi];
299                 }
300                 else
301                         return a[mid];
302         }
303 }
304
305 /**
306  * \note after a call to this function you can expect one of:
307  * - every node to left of a[n] are smaller or equal to it
308  * - every node to the right of a[n] are greater or equal to it */
309 static void partition_nth_element(BVHNode **a, int begin, int end, const int n, const int axis)
310 {
311         while (end - begin > 3) {
312                 const int cut = bvh_partition(a, begin, end, bvh_medianof3(a, begin, (begin + end) / 2, end - 1, axis), axis);
313                 if (cut <= n) {
314                         begin = cut;
315                 }
316                 else {
317                         end = cut;
318                 }
319         }
320         bvh_insertionsort(a, begin, end, axis);
321 }
322
323 #ifdef USE_SKIP_LINKS
324 static void build_skip_links(BVHTree *tree, BVHNode *node, BVHNode *left, BVHNode *right)
325 {
326         int i;
327
328         node->skip[0] = left;
329         node->skip[1] = right;
330
331         for (i = 0; i < node->totnode; i++) {
332                 if (i + 1 < node->totnode)
333                         build_skip_links(tree, node->children[i], left, node->children[i + 1]);
334                 else
335                         build_skip_links(tree, node->children[i], left, right);
336
337                 left = node->children[i];
338         }
339 }
340 #endif
341
342 /*
343  * BVHTree bounding volumes functions
344  */
345 static void create_kdop_hull(const BVHTree *tree, BVHNode *node, const float *co, int numpoints, int moving)
346 {
347         float newminmax;
348         float *bv = node->bv;
349         int k;
350         axis_t axis_iter;
351
352         /* don't init boudings for the moving case */
353         if (!moving) {
354                 node_minmax_init(tree, node);
355         }
356
357         for (k = 0; k < numpoints; k++) {
358                 /* for all Axes. */
359                 for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
360                         newminmax = dot_v3v3(&co[k * 3], bvhtree_kdop_axes[axis_iter]);
361                         if (newminmax < bv[2 * axis_iter])
362                                 bv[2 * axis_iter] = newminmax;
363                         if (newminmax > bv[(2 * axis_iter) + 1])
364                                 bv[(2 * axis_iter) + 1] = newminmax;
365                 }
366         }
367 }
368
369 /**
370  * \note depends on the fact that the BVH's for each face is already built
371  */
372 static void refit_kdop_hull(const BVHTree *tree, BVHNode *node, int start, int end)
373 {
374         float newmin, newmax;
375         float *__restrict bv = node->bv;
376         int j;
377         axis_t axis_iter;
378
379         node_minmax_init(tree, node);
380
381         for (j = start; j < end; j++) {
382                 float *__restrict node_bv = tree->nodes[j]->bv;
383
384                 /* for all Axes. */
385                 for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
386                         newmin = node_bv[(2 * axis_iter)];
387                         if ((newmin < bv[(2 * axis_iter)]))
388                                 bv[(2 * axis_iter)] = newmin;
389
390                         newmax = node_bv[(2 * axis_iter) + 1];
391                         if ((newmax > bv[(2 * axis_iter) + 1]))
392                                 bv[(2 * axis_iter) + 1] = newmax;
393                 }
394         }
395
396 }
397
398 /**
399  * only supports x,y,z axis in the moment
400  * but we should use a plain and simple function here for speed sake */
401 static char get_largest_axis(const float *bv)
402 {
403         float middle_point[3];
404
405         middle_point[0] = (bv[1]) - (bv[0]); /* x axis */
406         middle_point[1] = (bv[3]) - (bv[2]); /* y axis */
407         middle_point[2] = (bv[5]) - (bv[4]); /* z axis */
408         if (middle_point[0] > middle_point[1]) {
409                 if (middle_point[0] > middle_point[2])
410                         return 1;  /* max x axis */
411                 else
412                         return 5;  /* max z axis */
413         }
414         else {
415                 if (middle_point[1] > middle_point[2])
416                         return 3;  /* max y axis */
417                 else
418                         return 5;  /* max z axis */
419         }
420 }
421
422 /**
423  * bottom-up update of bvh node BV
424  * join the children on the parent BV */
425 static void node_join(BVHTree *tree, BVHNode *node)
426 {
427         int i;
428         axis_t axis_iter;
429
430         node_minmax_init(tree, node);
431
432         for (i = 0; i < tree->tree_type; i++) {
433                 if (node->children[i]) {
434                         for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
435                                 /* update minimum */
436                                 if (node->children[i]->bv[(2 * axis_iter)] < node->bv[(2 * axis_iter)])
437                                         node->bv[(2 * axis_iter)] = node->children[i]->bv[(2 * axis_iter)];
438
439                                 /* update maximum */
440                                 if (node->children[i]->bv[(2 * axis_iter) + 1] > node->bv[(2 * axis_iter) + 1])
441                                         node->bv[(2 * axis_iter) + 1] = node->children[i]->bv[(2 * axis_iter) + 1];
442                         }
443                 }
444                 else
445                         break;
446         }
447 }
448
449 #ifdef USE_PRINT_TREE
450
451 /**
452  * Debug and information functions
453  */
454
455 static void bvhtree_print_tree(BVHTree *tree, BVHNode *node, int depth)
456 {
457         int i;
458         axis_t axis_iter;
459
460         for (i = 0; i < depth; i++) printf(" ");
461         printf(" - %d (%ld): ", node->index, (long int)(node - tree->nodearray));
462         for (axis_iter = (axis_t)(2 * tree->start_axis);
463              axis_iter < (axis_t)(2 * tree->stop_axis);
464              axis_iter++)
465         {
466                 printf("%.3f ", node->bv[axis_iter]);
467         }
468         printf("\n");
469
470         for (i = 0; i < tree->tree_type; i++)
471                 if (node->children[i])
472                         bvhtree_print_tree(tree, node->children[i], depth + 1);
473 }
474
475 static void bvhtree_info(BVHTree *tree)
476 {
477         printf("BVHTree Info: tree_type = %d, axis = %d, epsilon = %f\n",
478                tree->tree_type, tree->axis, tree->epsilon);
479         printf("nodes = %d, branches = %d, leafs = %d\n",
480                tree->totbranch + tree->totleaf,  tree->totbranch, tree->totleaf);
481         printf("Memory per node = %ubytes\n",
482                (uint)(sizeof(BVHNode) + sizeof(BVHNode *) * tree->tree_type + sizeof(float) * tree->axis));
483         printf("BV memory = %ubytes\n",
484                (uint)MEM_allocN_len(tree->nodebv));
485
486         printf("Total memory = %ubytes\n",
487                (uint)(sizeof(BVHTree) +
488                       MEM_allocN_len(tree->nodes) +
489                       MEM_allocN_len(tree->nodearray) +
490                       MEM_allocN_len(tree->nodechild) +
491                       MEM_allocN_len(tree->nodebv)));
492
493         bvhtree_print_tree(tree, tree->nodes[tree->totleaf], 0);
494 }
495 #endif  /* USE_PRINT_TREE */
496
497 #ifdef USE_VERIFY_TREE
498
499 static void bvhtree_verify(BVHTree *tree)
500 {
501         int i, j, check = 0;
502
503         /* check the pointer list */
504         for (i = 0; i < tree->totleaf; i++) {
505                 if (tree->nodes[i]->parent == NULL) {
506                         printf("Leaf has no parent: %d\n", i);
507                 }
508                 else {
509                         for (j = 0; j < tree->tree_type; j++) {
510                                 if (tree->nodes[i]->parent->children[j] == tree->nodes[i])
511                                         check = 1;
512                         }
513                         if (!check) {
514                                 printf("Parent child relationship doesn't match: %d\n", i);
515                         }
516                         check = 0;
517                 }
518         }
519
520         /* check the leaf list */
521         for (i = 0; i < tree->totleaf; i++) {
522                 if (tree->nodearray[i].parent == NULL) {
523                         printf("Leaf has no parent: %d\n", i);
524                 }
525                 else {
526                         for (j = 0; j < tree->tree_type; j++) {
527                                 if (tree->nodearray[i].parent->children[j] == &tree->nodearray[i])
528                                         check = 1;
529                         }
530                         if (!check) {
531                                 printf("Parent child relationship doesn't match: %d\n", i);
532                         }
533                         check = 0;
534                 }
535         }
536
537         printf("branches: %d, leafs: %d, total: %d\n",
538                tree->totbranch, tree->totleaf, tree->totbranch + tree->totleaf);
539 }
540 #endif  /* USE_VERIFY_TREE */
541
542 /* Helper data and structures to build a min-leaf generalized implicit tree
543  * This code can be easily reduced
544  * (basicly this is only method to calculate pow(k, n) in O(1).. and stuff like that) */
545 typedef struct BVHBuildHelper {
546         int tree_type;              /* */
547         int totleafs;               /* */
548
549         int leafs_per_child[32];    /* Min number of leafs that are archievable from a node at depth N */
550         int branches_on_level[32];  /* Number of nodes at depth N (tree_type^N) */
551
552         int remain_leafs;           /* Number of leafs that are placed on the level that is not 100% filled */
553
554 } BVHBuildHelper;
555
556 static void build_implicit_tree_helper(const BVHTree *tree, BVHBuildHelper *data)
557 {
558         int depth = 0;
559         int remain;
560         int nnodes;
561
562         data->totleafs = tree->totleaf;
563         data->tree_type = tree->tree_type;
564
565         /* Calculate the smallest tree_type^n such that tree_type^n >= num_leafs */
566         for (data->leafs_per_child[0] = 1;
567              data->leafs_per_child[0] <  data->totleafs;
568              data->leafs_per_child[0] *= data->tree_type)
569         {
570                 /* pass */
571         }
572
573         data->branches_on_level[0] = 1;
574
575         for (depth = 1; (depth < 32) && data->leafs_per_child[depth - 1]; depth++) {
576                 data->branches_on_level[depth] = data->branches_on_level[depth - 1] * data->tree_type;
577                 data->leafs_per_child[depth] = data->leafs_per_child[depth - 1] / data->tree_type;
578         }
579
580         remain = data->totleafs - data->leafs_per_child[1];
581         nnodes = (remain + data->tree_type - 2) / (data->tree_type - 1);
582         data->remain_leafs = remain + nnodes;
583 }
584
585 // return the min index of all the leafs archivable with the given branch
586 static int implicit_leafs_index(const BVHBuildHelper *data, const int depth, const int child_index)
587 {
588         int min_leaf_index = child_index * data->leafs_per_child[depth - 1];
589         if (min_leaf_index <= data->remain_leafs)
590                 return min_leaf_index;
591         else if (data->leafs_per_child[depth])
592                 return data->totleafs - (data->branches_on_level[depth - 1] - child_index) * data->leafs_per_child[depth];
593         else
594                 return data->remain_leafs;
595 }
596
597 /**
598  * Generalized implicit tree build
599  *
600  * An implicit tree is a tree where its structure is implied, thus there is no need to store child pointers or indexs.
601  * Its possible to find the position of the child or the parent with simple maths (multiplication and adittion).
602  * This type of tree is for example used on heaps.. where node N has its childs at indexs N*2 and N*2+1.
603  *
604  * Although in this case the tree type is general.. and not know until runtime.
605  * tree_type stands for the maximum number of childs that a tree node can have.
606  * All tree types >= 2 are supported.
607  *
608  * Advantages of the used trees include:
609  *  - No need to store child/parent relations (they are implicit);
610  *  - Any node child always has an index greater than the parent;
611  *  - Brother nodes are sequential in memory;
612  *
613  *
614  * Some math relations derived for general implicit trees:
615  *
616  *   K = tree_type, ( 2 <= K )
617  *   ROOT = 1
618  *   N child of node A = A * K + (2 - K) + N, (0 <= N < K)
619  *
620  * Util methods:
621  *   TODO...
622  *    (looping elements, knowing if its a leaf or not.. etc...)
623  */
624
625 /* This functions returns the number of branches needed to have the requested number of leafs. */
626 static int implicit_needed_branches(int tree_type, int leafs)
627 {
628         return max_ii(1, (leafs + tree_type - 3) / (tree_type - 1));
629 }
630
631 /**
632  * This function handles the problem of "sorting" the leafs (along the split_axis).
633  *
634  * It arranges the elements in the given partitions such that:
635  *  - any element in partition N is less or equal to any element in partition N+1.
636  *  - if all elements are different all partition will get the same subset of elements
637  *    as if the array was sorted.
638  *
639  * partition P is described as the elements in the range ( nth[P], nth[P+1] ]
640  *
641  * TODO: This can be optimized a bit by doing a specialized nth_element instead of K nth_elements
642  */
643 static void split_leafs(BVHNode **leafs_array, const int nth[], const int partitions, const int split_axis)
644 {
645         int i;
646         for (i = 0; i < partitions - 1; i++) {
647                 if (nth[i] >= nth[partitions])
648                         break;
649
650                 partition_nth_element(leafs_array, nth[i], nth[partitions], nth[i + 1], split_axis);
651         }
652 }
653
654 typedef struct BVHDivNodesData {
655         const BVHTree *tree;
656         BVHNode *branches_array;
657         BVHNode **leafs_array;
658
659         int tree_type;
660         int tree_offset;
661
662         const BVHBuildHelper *data;
663
664         int depth;
665         int i;
666         int first_of_next_level;
667 } BVHDivNodesData;
668
669 static void non_recursive_bvh_div_nodes_task_cb(
670         void *__restrict userdata,
671         const int j,
672         const ParallelRangeTLS *__restrict UNUSED(tls))
673 {
674         BVHDivNodesData *data = userdata;
675
676         int k;
677         const int parent_level_index = j - data->i;
678         BVHNode *parent = &data->branches_array[j];
679         int nth_positions[MAX_TREETYPE + 1];
680         char split_axis;
681
682         int parent_leafs_begin = implicit_leafs_index(data->data, data->depth, parent_level_index);
683         int parent_leafs_end   = implicit_leafs_index(data->data, data->depth, parent_level_index + 1);
684
685         /* This calculates the bounding box of this branch
686          * and chooses the largest axis as the axis to divide leafs */
687         refit_kdop_hull(data->tree, parent, parent_leafs_begin, parent_leafs_end);
688         split_axis = get_largest_axis(parent->bv);
689
690         /* Save split axis (this can be used on raytracing to speedup the query time) */
691         parent->main_axis = split_axis / 2;
692
693         /* Split the childs along the split_axis, note: its not needed to sort the whole leafs array
694          * Only to assure that the elements are partitioned on a way that each child takes the elements
695          * it would take in case the whole array was sorted.
696          * Split_leafs takes care of that "sort" problem. */
697         nth_positions[0] = parent_leafs_begin;
698         nth_positions[data->tree_type] = parent_leafs_end;
699         for (k = 1; k < data->tree_type; k++) {
700                 const int child_index = j * data->tree_type + data->tree_offset + k;
701                 const int child_level_index = child_index - data->first_of_next_level; /* child level index */
702                 nth_positions[k] = implicit_leafs_index(data->data, data->depth + 1, child_level_index);
703         }
704
705         split_leafs(data->leafs_array, nth_positions, data->tree_type, split_axis);
706
707         /* Setup children and totnode counters
708          * Not really needed but currently most of BVH code relies on having an explicit children structure */
709         for (k = 0; k < data->tree_type; k++) {
710                 const int child_index = j * data->tree_type + data->tree_offset + k;
711                 const int child_level_index = child_index - data->first_of_next_level; /* child level index */
712
713                 const int child_leafs_begin = implicit_leafs_index(data->data, data->depth + 1, child_level_index);
714                 const int child_leafs_end   = implicit_leafs_index(data->data, data->depth + 1, child_level_index + 1);
715
716                 if (child_leafs_end - child_leafs_begin > 1) {
717                         parent->children[k] = &data->branches_array[child_index];
718                         parent->children[k]->parent = parent;
719                 }
720                 else if (child_leafs_end - child_leafs_begin == 1) {
721                         parent->children[k] = data->leafs_array[child_leafs_begin];
722                         parent->children[k]->parent = parent;
723                 }
724                 else {
725                         break;
726                 }
727         }
728         parent->totnode = (char)k;
729 }
730
731 /**
732  * This functions builds an optimal implicit tree from the given leafs.
733  * Where optimal stands for:
734  *  - The resulting tree will have the smallest number of branches;
735  *  - At most only one branch will have NULL childs;
736  *  - All leafs will be stored at level N or N+1.
737  *
738  * This function creates an implicit tree on branches_array, the leafs are given on the leafs_array.
739  *
740  * The tree is built per depth levels. First branches at depth 1.. then branches at depth 2.. etc..
741  * The reason is that we can build level N+1 from level N without any data dependencies.. thus it allows
742  * to use multithread building.
743  *
744  * To archive this is necessary to find how much leafs are accessible from a certain branch, BVHBuildHelper
745  * #implicit_needed_branches and #implicit_leafs_index are auxiliary functions to solve that "optimal-split".
746  */
747 static void non_recursive_bvh_div_nodes(
748         const BVHTree *tree, BVHNode *branches_array, BVHNode **leafs_array, int num_leafs)
749 {
750         int i;
751
752         const int tree_type   = tree->tree_type;
753         const int tree_offset = 2 - tree->tree_type; /* this value is 0 (on binary trees) and negative on the others */
754         const int num_branches = implicit_needed_branches(tree_type, num_leafs);
755
756         BVHBuildHelper data;
757         int depth;
758
759         {
760                 /* set parent from root node to NULL */
761                 BVHNode *root = &branches_array[1];
762                 root->parent = NULL;
763
764                 /* Most of bvhtree code relies on 1-leaf trees having at least one branch
765                  * We handle that special case here */
766                 if (num_leafs == 1) {
767                         refit_kdop_hull(tree, root, 0, num_leafs);
768                         root->main_axis = get_largest_axis(root->bv) / 2;
769                         root->totnode = 1;
770                         root->children[0] = leafs_array[0];
771                         root->children[0]->parent = root;
772                         return;
773                 }
774         }
775
776         build_implicit_tree_helper(tree, &data);
777
778         BVHDivNodesData cb_data = {
779                 .tree = tree, .branches_array = branches_array, .leafs_array = leafs_array,
780                 .tree_type = tree_type, .tree_offset = tree_offset, .data = &data,
781                 .first_of_next_level = 0, .depth = 0, .i = 0,
782         };
783
784         /* Loop tree levels (log N) loops */
785         for (i = 1, depth = 1; i <= num_branches; i = i * tree_type + tree_offset, depth++) {
786                 const int first_of_next_level = i * tree_type + tree_offset;
787                 const int i_stop = min_ii(first_of_next_level, num_branches + 1);  /* index of last branch on this level */
788
789                 /* Loop all branches on this level */
790                 cb_data.first_of_next_level = first_of_next_level;
791                 cb_data.i = i;
792                 cb_data.depth = depth;
793
794                 if (true) {
795                         ParallelRangeSettings settings;
796                         BLI_parallel_range_settings_defaults(&settings);
797                         settings.use_threading = (num_leafs > KDOPBVH_THREAD_LEAF_THRESHOLD);
798                         BLI_task_parallel_range(
799                                 i, i_stop,
800                                 &cb_data,
801                                 non_recursive_bvh_div_nodes_task_cb,
802                                 &settings);
803                 }
804                 else {
805                         /* Less hassle for debugging. */
806                         ParallelRangeTLS tls = {0};
807                         for (int i_task = i; i_task < i_stop; i_task++) {
808                                 non_recursive_bvh_div_nodes_task_cb(&cb_data, i_task, &tls);
809                         }
810                 }
811         }
812 }
813
814 /** \} */
815
816
817 /* -------------------------------------------------------------------- */
818
819 /** \name BLI_bvhtree API
820  * \{ */
821
822 /**
823  * \note many callers don't check for ``NULL`` return.
824  */
825 BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis)
826 {
827         BVHTree *tree;
828         int numnodes, i;
829
830         BLI_assert(tree_type >= 2 && tree_type <= MAX_TREETYPE);
831
832         tree = MEM_callocN(sizeof(BVHTree), "BVHTree");
833
834         /* tree epsilon must be >= FLT_EPSILON
835          * so that tangent rays can still hit a bounding volume..
836          * this bug would show up when casting a ray aligned with a kdop-axis and with an edge of 2 faces */
837         epsilon = max_ff(FLT_EPSILON, epsilon);
838
839         if (tree) {
840                 tree->epsilon = epsilon;
841                 tree->tree_type = tree_type;
842                 tree->axis = axis;
843
844                 if (axis == 26) {
845                         tree->start_axis = 0;
846                         tree->stop_axis = 13;
847                 }
848                 else if (axis == 18) {
849                         tree->start_axis = 7;
850                         tree->stop_axis = 13;
851                 }
852                 else if (axis == 14) {
853                         tree->start_axis = 0;
854                         tree->stop_axis = 7;
855                 }
856                 else if (axis == 8) { /* AABB */
857                         tree->start_axis = 0;
858                         tree->stop_axis = 4;
859                 }
860                 else if (axis == 6) { /* OBB */
861                         tree->start_axis = 0;
862                         tree->stop_axis = 3;
863                 }
864                 else {
865                         /* should never happen! */
866                         BLI_assert(0);
867
868                         goto fail;
869                 }
870
871
872                 /* Allocate arrays */
873                 numnodes = maxsize + implicit_needed_branches(tree_type, maxsize) + tree_type;
874
875                 tree->nodes = MEM_callocN(sizeof(BVHNode *) * (size_t)numnodes, "BVHNodes");
876                 tree->nodebv = MEM_callocN(sizeof(float) * (size_t)(axis * numnodes), "BVHNodeBV");
877                 tree->nodechild = MEM_callocN(sizeof(BVHNode *) * (size_t)(tree_type * numnodes), "BVHNodeBV");
878                 tree->nodearray = MEM_callocN(sizeof(BVHNode) * (size_t)numnodes, "BVHNodeArray");
879
880                 if (UNLIKELY((!tree->nodes) ||
881                              (!tree->nodebv) ||
882                              (!tree->nodechild) ||
883                              (!tree->nodearray)))
884                 {
885                         goto fail;
886                 }
887
888                 /* link the dynamic bv and child links */
889                 for (i = 0; i < numnodes; i++) {
890                         tree->nodearray[i].bv = &tree->nodebv[i * axis];
891                         tree->nodearray[i].children = &tree->nodechild[i * tree_type];
892                 }
893
894         }
895         return tree;
896
897
898 fail:
899         MEM_SAFE_FREE(tree->nodes);
900         MEM_SAFE_FREE(tree->nodebv);
901         MEM_SAFE_FREE(tree->nodechild);
902         MEM_SAFE_FREE(tree->nodearray);
903
904         MEM_freeN(tree);
905
906         return NULL;
907 }
908
909 void BLI_bvhtree_free(BVHTree *tree)
910 {
911         if (tree) {
912                 MEM_freeN(tree->nodes);
913                 MEM_freeN(tree->nodearray);
914                 MEM_freeN(tree->nodebv);
915                 MEM_freeN(tree->nodechild);
916                 MEM_freeN(tree);
917         }
918 }
919
920 void BLI_bvhtree_balance(BVHTree *tree)
921 {
922         BVHNode **leafs_array    = tree->nodes;
923
924         /* This function should only be called once
925          * (some big bug goes here if its being called more than once per tree) */
926         BLI_assert(tree->totbranch == 0);
927
928         /* Build the implicit tree */
929         non_recursive_bvh_div_nodes(tree, tree->nodearray + (tree->totleaf - 1), leafs_array, tree->totleaf);
930
931         /* current code expects the branches to be linked to the nodes array
932          * we perform that linkage here */
933         tree->totbranch = implicit_needed_branches(tree->tree_type, tree->totleaf);
934         for (int i = 0; i < tree->totbranch; i++) {
935                 tree->nodes[tree->totleaf + i] = &tree->nodearray[tree->totleaf + i];
936         }
937
938 #ifdef USE_SKIP_LINKS
939         build_skip_links(tree, tree->nodes[tree->totleaf], NULL, NULL);
940 #endif
941
942 #ifdef USE_VERIFY_TREE
943         bvhtree_verify(tree);
944 #endif
945
946 #ifdef USE_PRINT_TREE
947         bvhtree_info(tree);
948 #endif
949 }
950
951 void BLI_bvhtree_insert(BVHTree *tree, int index, const float co[3], int numpoints)
952 {
953         axis_t axis_iter;
954         BVHNode *node = NULL;
955
956         /* insert should only possible as long as tree->totbranch is 0 */
957         BLI_assert(tree->totbranch <= 0);
958         BLI_assert((size_t)tree->totleaf < MEM_allocN_len(tree->nodes) / sizeof(*(tree->nodes)));
959
960         node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]);
961         tree->totleaf++;
962
963         create_kdop_hull(tree, node, co, numpoints, 0);
964         node->index = index;
965
966         /* inflate the bv with some epsilon */
967         for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
968                 node->bv[(2 * axis_iter)] -= tree->epsilon; /* minimum */
969                 node->bv[(2 * axis_iter) + 1] += tree->epsilon; /* maximum */
970         }
971 }
972
973
974 /* call before BLI_bvhtree_update_tree() */
975 bool BLI_bvhtree_update_node(BVHTree *tree, int index, const float co[3], const float co_moving[3], int numpoints)
976 {
977         BVHNode *node = NULL;
978         axis_t axis_iter;
979
980         /* check if index exists */
981         if (index > tree->totleaf)
982                 return false;
983
984         node = tree->nodearray + index;
985
986         create_kdop_hull(tree, node, co, numpoints, 0);
987
988         if (co_moving)
989                 create_kdop_hull(tree, node, co_moving, numpoints, 1);
990
991         /* inflate the bv with some epsilon */
992         for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
993                 node->bv[(2 * axis_iter)]     -= tree->epsilon; /* minimum */
994                 node->bv[(2 * axis_iter) + 1] += tree->epsilon; /* maximum */
995         }
996
997         return true;
998 }
999
1000 /* call BLI_bvhtree_update_node() first for every node/point/triangle */
1001 void BLI_bvhtree_update_tree(BVHTree *tree)
1002 {
1003         /* Update bottom=>top
1004          * TRICKY: the way we build the tree all the childs have an index greater than the parent
1005          * This allows us todo a bottom up update by starting on the bigger numbered branch */
1006
1007         BVHNode **root  = tree->nodes + tree->totleaf;
1008         BVHNode **index = tree->nodes + tree->totleaf + tree->totbranch - 1;
1009
1010         for (; index >= root; index--)
1011                 node_join(tree, *index);
1012 }
1013 /**
1014  * Number of times #BLI_bvhtree_insert has been called.
1015  * mainly useful for asserts functions to check we added the correct number.
1016  */
1017 int BLI_bvhtree_get_len(const BVHTree *tree)
1018 {
1019         return tree->totleaf;
1020 }
1021
1022 /**
1023  * Maximum number of children that a node can have.
1024  */
1025 int BLI_bvhtree_get_tree_type(const BVHTree *tree)
1026 {
1027         return tree->tree_type;
1028 }
1029
1030 float BLI_bvhtree_get_epsilon(const BVHTree *tree)
1031 {
1032         return tree->epsilon;
1033 }
1034
1035 /** \} */
1036
1037
1038 /* -------------------------------------------------------------------- */
1039
1040 /** \name BLI_bvhtree_overlap
1041  * \{ */
1042
1043 /**
1044  * overlap - is it possible for 2 bv's to collide ?
1045  */
1046 static bool tree_overlap_test(const BVHNode *node1, const BVHNode *node2, axis_t start_axis, axis_t stop_axis)
1047 {
1048         const float *bv1     = node1->bv + (start_axis << 1);
1049         const float *bv2     = node2->bv + (start_axis << 1);
1050         const float *bv1_end = node1->bv + (stop_axis  << 1);
1051
1052         /* test all axis if min + max overlap */
1053         for (; bv1 != bv1_end; bv1 += 2, bv2 += 2) {
1054                 if ((bv1[0] > bv2[1]) || (bv2[0] > bv1[1])) {
1055                         return 0;
1056                 }
1057         }
1058
1059         return 1;
1060 }
1061
1062 static void tree_overlap_traverse(
1063         BVHOverlapData_Thread *data_thread,
1064         const BVHNode *node1, const BVHNode *node2)
1065 {
1066         BVHOverlapData_Shared *data = data_thread->shared;
1067         int j;
1068
1069         if (tree_overlap_test(node1, node2, data->start_axis, data->stop_axis)) {
1070                 /* check if node1 is a leaf */
1071                 if (!node1->totnode) {
1072                         /* check if node2 is a leaf */
1073                         if (!node2->totnode) {
1074                                 BVHTreeOverlap *overlap;
1075
1076                                 if (UNLIKELY(node1 == node2)) {
1077                                         return;
1078                                 }
1079
1080                                 /* both leafs, insert overlap! */
1081                                 overlap = BLI_stack_push_r(data_thread->overlap);
1082                                 overlap->indexA = node1->index;
1083                                 overlap->indexB = node2->index;
1084                         }
1085                         else {
1086                                 for (j = 0; j < data->tree2->tree_type; j++) {
1087                                         if (node2->children[j]) {
1088                                                 tree_overlap_traverse(data_thread, node1, node2->children[j]);
1089                                         }
1090                                 }
1091                         }
1092                 }
1093                 else {
1094                         for (j = 0; j < data->tree2->tree_type; j++) {
1095                                 if (node1->children[j]) {
1096                                         tree_overlap_traverse(data_thread, node1->children[j], node2);
1097                                 }
1098                         }
1099                 }
1100         }
1101 }
1102
1103 /**
1104  * a version of #tree_overlap_traverse that runs a callback to check if the nodes really intersect.
1105  */
1106 static void tree_overlap_traverse_cb(
1107         BVHOverlapData_Thread *data_thread,
1108         const BVHNode *node1, const BVHNode *node2)
1109 {
1110         BVHOverlapData_Shared *data = data_thread->shared;
1111         int j;
1112
1113         if (tree_overlap_test(node1, node2, data->start_axis, data->stop_axis)) {
1114                 /* check if node1 is a leaf */
1115                 if (!node1->totnode) {
1116                         /* check if node2 is a leaf */
1117                         if (!node2->totnode) {
1118                                 BVHTreeOverlap *overlap;
1119
1120                                 if (UNLIKELY(node1 == node2)) {
1121                                         return;
1122                                 }
1123
1124                                 /* only difference to tree_overlap_traverse! */
1125                                 if (data->callback(data->userdata, node1->index, node2->index, data_thread->thread)) {
1126                                         /* both leafs, insert overlap! */
1127                                         overlap = BLI_stack_push_r(data_thread->overlap);
1128                                         overlap->indexA = node1->index;
1129                                         overlap->indexB = node2->index;
1130                                 }
1131                         }
1132                         else {
1133                                 for (j = 0; j < data->tree2->tree_type; j++) {
1134                                         if (node2->children[j]) {
1135                                                 tree_overlap_traverse_cb(data_thread, node1, node2->children[j]);
1136                                         }
1137                                 }
1138                         }
1139                 }
1140                 else {
1141                         for (j = 0; j < data->tree2->tree_type; j++) {
1142                                 if (node1->children[j]) {
1143                                         tree_overlap_traverse_cb(data_thread, node1->children[j], node2);
1144                                 }
1145                         }
1146                 }
1147         }
1148 }
1149
1150 /**
1151  * Use to check the total number of threads #BLI_bvhtree_overlap will use.
1152  *
1153  * \warning Must be the first tree passed to #BLI_bvhtree_overlap!
1154  */
1155 int BLI_bvhtree_overlap_thread_num(const BVHTree *tree)
1156 {
1157         return (int)MIN2(tree->tree_type, tree->nodes[tree->totleaf]->totnode);
1158 }
1159
1160 static void bvhtree_overlap_task_cb(
1161         void *__restrict userdata,
1162         const int j,
1163         const ParallelRangeTLS *__restrict UNUSED(tls))
1164 {
1165         BVHOverlapData_Thread *data = &((BVHOverlapData_Thread *)userdata)[j];
1166         BVHOverlapData_Shared *data_shared = data->shared;
1167
1168         if (data_shared->callback) {
1169                 tree_overlap_traverse_cb(
1170                             data, data_shared->tree1->nodes[data_shared->tree1->totleaf]->children[j],
1171                             data_shared->tree2->nodes[data_shared->tree2->totleaf]);
1172         }
1173         else {
1174                 tree_overlap_traverse(
1175                             data, data_shared->tree1->nodes[data_shared->tree1->totleaf]->children[j],
1176                             data_shared->tree2->nodes[data_shared->tree2->totleaf]);
1177         }
1178 }
1179
1180 BVHTreeOverlap *BLI_bvhtree_overlap(
1181         const BVHTree *tree1, const BVHTree *tree2, uint *r_overlap_tot,
1182         /* optional callback to test the overlap before adding (must be thread-safe!) */
1183         BVHTree_OverlapCallback callback, void *userdata)
1184 {
1185         const int thread_num = BLI_bvhtree_overlap_thread_num(tree1);
1186         int j;
1187         size_t total = 0;
1188         BVHTreeOverlap *overlap = NULL, *to = NULL;
1189         BVHOverlapData_Shared data_shared;
1190         BVHOverlapData_Thread *data = BLI_array_alloca(data, (size_t)thread_num);
1191         axis_t start_axis, stop_axis;
1192
1193         /* check for compatibility of both trees (can't compare 14-DOP with 18-DOP) */
1194         if (UNLIKELY((tree1->axis != tree2->axis) &&
1195                      (tree1->axis == 14 || tree2->axis == 14) &&
1196                      (tree1->axis == 18 || tree2->axis == 18)))
1197         {
1198                 BLI_assert(0);
1199                 return NULL;
1200         }
1201
1202         start_axis = min_axis(tree1->start_axis, tree2->start_axis);
1203         stop_axis  = min_axis(tree1->stop_axis,  tree2->stop_axis);
1204
1205         /* fast check root nodes for collision before doing big splitting + traversal */
1206         if (!tree_overlap_test(tree1->nodes[tree1->totleaf], tree2->nodes[tree2->totleaf], start_axis, stop_axis)) {
1207                 return NULL;
1208         }
1209
1210         data_shared.tree1 = tree1;
1211         data_shared.tree2 = tree2;
1212         data_shared.start_axis = start_axis;
1213         data_shared.stop_axis = stop_axis;
1214
1215         /* can be NULL */
1216         data_shared.callback = callback;
1217         data_shared.userdata = userdata;
1218
1219         for (j = 0; j < thread_num; j++) {
1220                 /* init BVHOverlapData_Thread */
1221                 data[j].shared = &data_shared;
1222                 data[j].overlap = BLI_stack_new(sizeof(BVHTreeOverlap), __func__);
1223
1224                 /* for callback */
1225                 data[j].thread = j;
1226         }
1227
1228         ParallelRangeSettings settings;
1229         BLI_parallel_range_settings_defaults(&settings);
1230         settings.use_threading = (tree1->totleaf > KDOPBVH_THREAD_LEAF_THRESHOLD);
1231         BLI_task_parallel_range(
1232                     0, thread_num,
1233                     data,
1234                     bvhtree_overlap_task_cb,
1235                     &settings);
1236
1237         for (j = 0; j < thread_num; j++)
1238                 total += BLI_stack_count(data[j].overlap);
1239
1240         to = overlap = MEM_mallocN(sizeof(BVHTreeOverlap) * total, "BVHTreeOverlap");
1241
1242         for (j = 0; j < thread_num; j++) {
1243                 uint count = (uint)BLI_stack_count(data[j].overlap);
1244                 BLI_stack_pop_n(data[j].overlap, to, count);
1245                 BLI_stack_free(data[j].overlap);
1246                 to += count;
1247         }
1248
1249         *r_overlap_tot = (uint)total;
1250         return overlap;
1251 }
1252
1253 /** \} */
1254
1255
1256 /* -------------------------------------------------------------------- */
1257
1258 /** \name BLI_bvhtree_find_nearest
1259  * \{ */
1260
1261 /* Determines the nearest point of the given node BV. Returns the squared distance to that point. */
1262 static float calc_nearest_point_squared(const float proj[3], BVHNode *node, float nearest[3])
1263 {
1264         int i;
1265         const float *bv = node->bv;
1266
1267         /* nearest on AABB hull */
1268         for (i = 0; i != 3; i++, bv += 2) {
1269                 if (bv[0] > proj[i])
1270                         nearest[i] = bv[0];
1271                 else if (bv[1] < proj[i])
1272                         nearest[i] = bv[1];
1273                 else
1274                         nearest[i] = proj[i];
1275         }
1276
1277         return len_squared_v3v3(proj, nearest);
1278 }
1279
1280 /* TODO: use a priority queue to reduce the number of nodes looked on */
1281 static void dfs_find_nearest_dfs(BVHNearestData *data, BVHNode *node)
1282 {
1283         if (node->totnode == 0) {
1284                 if (data->callback)
1285                         data->callback(data->userdata, node->index, data->co, &data->nearest);
1286                 else {
1287                         data->nearest.index = node->index;
1288                         data->nearest.dist_sq = calc_nearest_point_squared(data->proj, node, data->nearest.co);
1289                 }
1290         }
1291         else {
1292                 /* Better heuristic to pick the closest node to dive on */
1293                 int i;
1294                 float nearest[3];
1295
1296                 if (data->proj[node->main_axis] <= node->children[0]->bv[node->main_axis * 2 + 1]) {
1297
1298                         for (i = 0; i != node->totnode; i++) {
1299                                 if (calc_nearest_point_squared(data->proj, node->children[i], nearest) >= data->nearest.dist_sq)
1300                                         continue;
1301                                 dfs_find_nearest_dfs(data, node->children[i]);
1302                         }
1303                 }
1304                 else {
1305                         for (i = node->totnode - 1; i >= 0; i--) {
1306                                 if (calc_nearest_point_squared(data->proj, node->children[i], nearest) >= data->nearest.dist_sq)
1307                                         continue;
1308                                 dfs_find_nearest_dfs(data, node->children[i]);
1309                         }
1310                 }
1311         }
1312 }
1313
1314 static void dfs_find_nearest_begin(BVHNearestData *data, BVHNode *node)
1315 {
1316         float nearest[3], dist_sq;
1317         dist_sq = calc_nearest_point_squared(data->proj, node, nearest);
1318         if (dist_sq >= data->nearest.dist_sq) {
1319                 return;
1320         }
1321         dfs_find_nearest_dfs(data, node);
1322 }
1323
1324 int BLI_bvhtree_find_nearest(
1325         BVHTree *tree, const float co[3], BVHTreeNearest *nearest,
1326         BVHTree_NearestPointCallback callback, void *userdata)
1327 {
1328         axis_t axis_iter;
1329
1330         BVHNearestData data;
1331         BVHNode *root = tree->nodes[tree->totleaf];
1332
1333         /* init data to search */
1334         data.tree = tree;
1335         data.co = co;
1336
1337         data.callback = callback;
1338         data.userdata = userdata;
1339
1340         for (axis_iter = data.tree->start_axis; axis_iter != data.tree->stop_axis; axis_iter++) {
1341                 data.proj[axis_iter] = dot_v3v3(data.co, bvhtree_kdop_axes[axis_iter]);
1342         }
1343
1344         if (nearest) {
1345                 memcpy(&data.nearest, nearest, sizeof(*nearest));
1346         }
1347         else {
1348                 data.nearest.index = -1;
1349                 data.nearest.dist_sq = FLT_MAX;
1350         }
1351
1352         /* dfs search */
1353         if (root)
1354                 dfs_find_nearest_begin(&data, root);
1355
1356         /* copy back results */
1357         if (nearest) {
1358                 memcpy(nearest, &data.nearest, sizeof(*nearest));
1359         }
1360
1361         return data.nearest.index;
1362 }
1363
1364 /** \} */
1365
1366
1367 /* -------------------------------------------------------------------- */
1368
1369 /** \name BLI_bvhtree_ray_cast
1370  *
1371  * raycast is done by performing a DFS on the BVHTree and saving the closest hit.
1372  *
1373  * \{ */
1374
1375
1376 /* Determines the distance that the ray must travel to hit the bounding volume of the given node */
1377 static float ray_nearest_hit(const BVHRayCastData *data, const float bv[6])
1378 {
1379         int i;
1380
1381         float low = 0, upper = data->hit.dist;
1382
1383         for (i = 0; i != 3; i++, bv += 2) {
1384                 if (data->ray_dot_axis[i] == 0.0f) {
1385                         /* axis aligned ray */
1386                         if (data->ray.origin[i] < bv[0] - data->ray.radius ||
1387                             data->ray.origin[i] > bv[1] + data->ray.radius)
1388                         {
1389                                 return FLT_MAX;
1390                         }
1391                 }
1392                 else {
1393                         float ll = (bv[0] - data->ray.radius - data->ray.origin[i]) / data->ray_dot_axis[i];
1394                         float lu = (bv[1] + data->ray.radius - data->ray.origin[i]) / data->ray_dot_axis[i];
1395
1396                         if (data->ray_dot_axis[i] > 0.0f) {
1397                                 if (ll > low) low = ll;
1398                                 if (lu < upper) upper = lu;
1399                         }
1400                         else {
1401                                 if (lu > low) low = lu;
1402                                 if (ll < upper) upper = ll;
1403                         }
1404
1405                         if (low > upper) return FLT_MAX;
1406                 }
1407         }
1408         return low;
1409 }
1410
1411 /**
1412  * Determines the distance that the ray must travel to hit the bounding volume of the given node
1413  * Based on Tactical Optimization of Ray/Box Intersection, by Graham Fyffe
1414  * [http://tog.acm.org/resources/RTNews/html/rtnv21n1.html#art9]
1415  *
1416  * TODO this doesn't take data->ray.radius into consideration */
1417 static float fast_ray_nearest_hit(const BVHRayCastData *data, const BVHNode *node)
1418 {
1419         const float *bv = node->bv;
1420
1421         float t1x = (bv[data->index[0]] - data->ray.origin[0]) * data->idot_axis[0];
1422         float t2x = (bv[data->index[1]] - data->ray.origin[0]) * data->idot_axis[0];
1423         float t1y = (bv[data->index[2]] - data->ray.origin[1]) * data->idot_axis[1];
1424         float t2y = (bv[data->index[3]] - data->ray.origin[1]) * data->idot_axis[1];
1425         float t1z = (bv[data->index[4]] - data->ray.origin[2]) * data->idot_axis[2];
1426         float t2z = (bv[data->index[5]] - data->ray.origin[2]) * data->idot_axis[2];
1427
1428         if ((t1x > t2y || t2x < t1y || t1x > t2z || t2x < t1z || t1y > t2z || t2y < t1z) ||
1429             (t2x < 0.0f || t2y < 0.0f || t2z < 0.0f) ||
1430             (t1x > data->hit.dist || t1y > data->hit.dist || t1z > data->hit.dist))
1431         {
1432                 return FLT_MAX;
1433         }
1434         else {
1435                 return max_fff(t1x, t1y, t1z);
1436         }
1437 }
1438
1439 static void dfs_raycast(BVHRayCastData *data, BVHNode *node)
1440 {
1441         int i;
1442
1443         /* ray-bv is really fast.. and simple tests revealed its worth to test it
1444          * before calling the ray-primitive functions */
1445         /* XXX: temporary solution for particles until fast_ray_nearest_hit supports ray.radius */
1446         float dist = (data->ray.radius == 0.0f) ? fast_ray_nearest_hit(data, node) : ray_nearest_hit(data, node->bv);
1447         if (dist >= data->hit.dist) {
1448                 return;
1449         }
1450
1451         if (node->totnode == 0) {
1452                 if (data->callback) {
1453                         data->callback(data->userdata, node->index, &data->ray, &data->hit);
1454                 }
1455                 else {
1456                         data->hit.index = node->index;
1457                         data->hit.dist  = dist;
1458                         madd_v3_v3v3fl(data->hit.co, data->ray.origin, data->ray.direction, dist);
1459                 }
1460         }
1461         else {
1462                 /* pick loop direction to dive into the tree (based on ray direction and split axis) */
1463                 if (data->ray_dot_axis[node->main_axis] > 0.0f) {
1464                         for (i = 0; i != node->totnode; i++) {
1465                                 dfs_raycast(data, node->children[i]);
1466                         }
1467                 }
1468                 else {
1469                         for (i = node->totnode - 1; i >= 0; i--) {
1470                                 dfs_raycast(data, node->children[i]);
1471                         }
1472                 }
1473         }
1474 }
1475
1476 /**
1477  * A version of #dfs_raycast with minor changes to reset the index & dist each ray cast.
1478  */
1479 static void dfs_raycast_all(BVHRayCastData *data, BVHNode *node)
1480 {
1481         int i;
1482
1483         /* ray-bv is really fast.. and simple tests revealed its worth to test it
1484          * before calling the ray-primitive functions */
1485         /* XXX: temporary solution for particles until fast_ray_nearest_hit supports ray.radius */
1486         float dist = (data->ray.radius == 0.0f) ? fast_ray_nearest_hit(data, node) : ray_nearest_hit(data, node->bv);
1487         if (dist >= data->hit.dist) {
1488                 return;
1489         }
1490
1491         if (node->totnode == 0) {
1492                 /* no need to check for 'data->callback' (using 'all' only makes sense with a callback). */
1493                 dist = data->hit.dist;
1494                 data->callback(data->userdata, node->index, &data->ray, &data->hit);
1495                 data->hit.index = -1;
1496                 data->hit.dist = dist;
1497         }
1498         else {
1499                 /* pick loop direction to dive into the tree (based on ray direction and split axis) */
1500                 if (data->ray_dot_axis[node->main_axis] > 0.0f) {
1501                         for (i = 0; i != node->totnode; i++) {
1502                                 dfs_raycast_all(data, node->children[i]);
1503                         }
1504                 }
1505                 else {
1506                         for (i = node->totnode - 1; i >= 0; i--) {
1507                                 dfs_raycast_all(data, node->children[i]);
1508                         }
1509                 }
1510         }
1511 }
1512
1513 static void bvhtree_ray_cast_data_precalc(BVHRayCastData *data, int flag)
1514 {
1515         int i;
1516
1517         for (i = 0; i < 3; i++) {
1518                 data->ray_dot_axis[i] = dot_v3v3(data->ray.direction, bvhtree_kdop_axes[i]);
1519                 data->idot_axis[i] = 1.0f / data->ray_dot_axis[i];
1520
1521                 if (fabsf(data->ray_dot_axis[i]) < FLT_EPSILON) {
1522                         data->ray_dot_axis[i] = 0.0;
1523                 }
1524                 data->index[2 * i] = data->idot_axis[i] < 0.0f ? 1 : 0;
1525                 data->index[2 * i + 1] = 1 - data->index[2 * i];
1526                 data->index[2 * i]   += 2 * i;
1527                 data->index[2 * i + 1] += 2 * i;
1528         }
1529
1530 #ifdef USE_KDOPBVH_WATERTIGHT
1531         if (flag & BVH_RAYCAST_WATERTIGHT) {
1532                 isect_ray_tri_watertight_v3_precalc(&data->isect_precalc, data->ray.direction);
1533                 data->ray.isect_precalc = &data->isect_precalc;
1534         }
1535         else {
1536                 data->ray.isect_precalc = NULL;
1537         }
1538 #else
1539         UNUSED_VARS(flag);
1540 #endif
1541 }
1542
1543 int BLI_bvhtree_ray_cast_ex(
1544         BVHTree *tree, const float co[3], const float dir[3], float radius, BVHTreeRayHit *hit,
1545         BVHTree_RayCastCallback callback, void *userdata,
1546         int flag)
1547 {
1548         BVHRayCastData data;
1549         BVHNode *root = tree->nodes[tree->totleaf];
1550
1551         BLI_ASSERT_UNIT_V3(dir);
1552
1553         data.tree = tree;
1554
1555         data.callback = callback;
1556         data.userdata = userdata;
1557
1558         copy_v3_v3(data.ray.origin,    co);
1559         copy_v3_v3(data.ray.direction, dir);
1560         data.ray.radius = radius;
1561
1562         bvhtree_ray_cast_data_precalc(&data, flag);
1563
1564         if (hit) {
1565                 memcpy(&data.hit, hit, sizeof(*hit));
1566         }
1567         else {
1568                 data.hit.index = -1;
1569                 data.hit.dist = BVH_RAYCAST_DIST_MAX;
1570         }
1571
1572         if (root) {
1573                 dfs_raycast(&data, root);
1574 //              iterative_raycast(&data, root);
1575         }
1576
1577
1578         if (hit)
1579                 memcpy(hit, &data.hit, sizeof(*hit));
1580
1581         return data.hit.index;
1582 }
1583
1584 int BLI_bvhtree_ray_cast(
1585         BVHTree *tree, const float co[3], const float dir[3], float radius, BVHTreeRayHit *hit,
1586         BVHTree_RayCastCallback callback, void *userdata)
1587 {
1588         return BLI_bvhtree_ray_cast_ex(tree, co, dir, radius, hit, callback, userdata, BVH_RAYCAST_DEFAULT);
1589 }
1590
1591 float BLI_bvhtree_bb_raycast(const float bv[6], const float light_start[3], const float light_end[3], float pos[3])
1592 {
1593         BVHRayCastData data;
1594         float dist;
1595
1596         data.hit.dist = BVH_RAYCAST_DIST_MAX;
1597
1598         /* get light direction */
1599         sub_v3_v3v3(data.ray.direction, light_end, light_start);
1600
1601         data.ray.radius = 0.0;
1602
1603         copy_v3_v3(data.ray.origin, light_start);
1604
1605         normalize_v3(data.ray.direction);
1606         copy_v3_v3(data.ray_dot_axis, data.ray.direction);
1607
1608         dist = ray_nearest_hit(&data, bv);
1609
1610         madd_v3_v3v3fl(pos, light_start, data.ray.direction, dist);
1611
1612         return dist;
1613
1614 }
1615
1616 /**
1617  * Calls the callback for every ray intersection
1618  *
1619  * \note Using a \a callback which resets or never sets the #BVHTreeRayHit index & dist works too,
1620  * however using this function means existing generic callbacks can be used from custom callbacks without
1621  * having to handle resetting the hit beforehand.
1622  * It also avoid redundant argument and return value which aren't meaningful when collecting multiple hits.
1623  */
1624 void BLI_bvhtree_ray_cast_all_ex(
1625         BVHTree *tree, const float co[3], const float dir[3], float radius, float hit_dist,
1626         BVHTree_RayCastCallback callback, void *userdata,
1627         int flag)
1628 {
1629         BVHRayCastData data;
1630         BVHNode *root = tree->nodes[tree->totleaf];
1631
1632         BLI_ASSERT_UNIT_V3(dir);
1633         BLI_assert(callback != NULL);
1634
1635         data.tree = tree;
1636
1637         data.callback = callback;
1638         data.userdata = userdata;
1639
1640         copy_v3_v3(data.ray.origin,    co);
1641         copy_v3_v3(data.ray.direction, dir);
1642         data.ray.radius = radius;
1643
1644         bvhtree_ray_cast_data_precalc(&data, flag);
1645
1646         data.hit.index = -1;
1647         data.hit.dist = hit_dist;
1648
1649         if (root) {
1650                 dfs_raycast_all(&data, root);
1651         }
1652 }
1653
1654 void BLI_bvhtree_ray_cast_all(
1655         BVHTree *tree, const float co[3], const float dir[3], float radius, float hit_dist,
1656         BVHTree_RayCastCallback callback, void *userdata)
1657 {
1658         BLI_bvhtree_ray_cast_all_ex(tree, co, dir, radius, hit_dist, callback, userdata, BVH_RAYCAST_DEFAULT);
1659 }
1660
1661 /** \} */
1662
1663 /* -------------------------------------------------------------------- */
1664
1665 /** \name BLI_bvhtree_range_query
1666  *
1667  * Allocs and fills an array with the indexs of node that are on the given spherical range (center, radius).
1668  * Returns the size of the array.
1669  *
1670  * \{ */
1671
1672 typedef struct RangeQueryData {
1673         BVHTree *tree;
1674         const float *center;
1675         float radius_sq;  /* squared radius */
1676
1677         int hits;
1678
1679         BVHTree_RangeQuery callback;
1680         void *userdata;
1681 } RangeQueryData;
1682
1683
1684 static void dfs_range_query(RangeQueryData *data, BVHNode *node)
1685 {
1686         if (node->totnode == 0) {
1687 #if 0   /*UNUSED*/
1688                 /* Calculate the node min-coords (if the node was a point then this is the point coordinates) */
1689                 float co[3];
1690                 co[0] = node->bv[0];
1691                 co[1] = node->bv[2];
1692                 co[2] = node->bv[4];
1693 #endif
1694         }
1695         else {
1696                 int i;
1697                 for (i = 0; i != node->totnode; i++) {
1698                         float nearest[3];
1699                         float dist_sq = calc_nearest_point_squared(data->center, node->children[i], nearest);
1700                         if (dist_sq < data->radius_sq) {
1701                                 /* Its a leaf.. call the callback */
1702                                 if (node->children[i]->totnode == 0) {
1703                                         data->hits++;
1704                                         data->callback(data->userdata, node->children[i]->index, data->center, dist_sq);
1705                                 }
1706                                 else
1707                                         dfs_range_query(data, node->children[i]);
1708                         }
1709                 }
1710         }
1711 }
1712
1713 int BLI_bvhtree_range_query(
1714         BVHTree *tree, const float co[3], float radius,
1715         BVHTree_RangeQuery callback, void *userdata)
1716 {
1717         BVHNode *root = tree->nodes[tree->totleaf];
1718
1719         RangeQueryData data;
1720         data.tree = tree;
1721         data.center = co;
1722         data.radius_sq = radius * radius;
1723         data.hits = 0;
1724
1725         data.callback = callback;
1726         data.userdata = userdata;
1727
1728         if (root != NULL) {
1729                 float nearest[3];
1730                 float dist_sq = calc_nearest_point_squared(data.center, root, nearest);
1731                 if (dist_sq < data.radius_sq) {
1732                         /* Its a leaf.. call the callback */
1733                         if (root->totnode == 0) {
1734                                 data.hits++;
1735                                 data.callback(data.userdata, root->index, co, dist_sq);
1736                         }
1737                         else
1738                                 dfs_range_query(&data, root);
1739                 }
1740         }
1741
1742         return data.hits;
1743 }
1744
1745 /** \} */
1746
1747
1748 /* -------------------------------------------------------------------- */
1749
1750 /** \name BLI_bvhtree_nearest_projected
1751 * \{ */
1752
1753 static void bvhtree_nearest_projected_dfs_recursive(
1754         BVHNearestProjectedData *__restrict data, const BVHNode *node)
1755 {
1756         if (node->totnode == 0) {
1757                 if (data->callback) {
1758                         data->callback(
1759                                 data->userdata, node->index, &data->precalc,
1760                                 NULL, 0,
1761                                 &data->nearest);
1762                 }
1763                 else {
1764                         data->nearest.index = node->index;
1765                         data->nearest.dist_sq = dist_squared_to_projected_aabb(
1766                                 &data->precalc,
1767                                 (float[3]) {node->bv[0], node->bv[2], node->bv[4]},
1768                                 (float[3]) {node->bv[1], node->bv[3], node->bv[5]},
1769                                 data->closest_axis);
1770                 }
1771         }
1772         else {
1773                 /* First pick the closest node to recurse into */
1774                 if (data->closest_axis[node->main_axis]) {
1775                         for (int i = 0; i != node->totnode; i++) {
1776                                 const float *bv = node->children[i]->bv;
1777
1778                                 if (dist_squared_to_projected_aabb(
1779                                         &data->precalc,
1780                                         (float[3]) {bv[0], bv[2], bv[4]},
1781                                         (float[3]) {bv[1], bv[3], bv[5]},
1782                                         data->closest_axis) <= data->nearest.dist_sq)
1783                                 {
1784                                         bvhtree_nearest_projected_dfs_recursive(data, node->children[i]);
1785                                 }
1786                         }
1787                 }
1788                 else {
1789                         for (int i = node->totnode; i--;) {
1790                                 const float *bv = node->children[i]->bv;
1791
1792                                 if (dist_squared_to_projected_aabb(
1793                                         &data->precalc,
1794                                         (float[3]) {bv[0], bv[2], bv[4]},
1795                                         (float[3]) {bv[1], bv[3], bv[5]},
1796                                         data->closest_axis) <= data->nearest.dist_sq)
1797                                 {
1798                                         bvhtree_nearest_projected_dfs_recursive(data, node->children[i]);
1799                                 }
1800                         }
1801                 }
1802         }
1803 }
1804
1805 static void bvhtree_nearest_projected_with_clipplane_test_dfs_recursive(
1806         BVHNearestProjectedData *__restrict data, const BVHNode *node)
1807 {
1808         if (node->totnode == 0) {
1809                 if (data->callback) {
1810                         data->callback(
1811                                 data->userdata, node->index, &data->precalc,
1812                                 data->clip_plane, data->clip_plane_len,
1813                                 &data->nearest);
1814                 }
1815                 else {
1816                         data->nearest.index = node->index;
1817                         data->nearest.dist_sq = dist_squared_to_projected_aabb(
1818                                 &data->precalc,
1819                                 (float[3]) {node->bv[0], node->bv[2], node->bv[4]},
1820                                 (float[3]) {node->bv[1], node->bv[3], node->bv[5]},
1821                                 data->closest_axis);
1822                 }
1823         }
1824         else {
1825                 /* First pick the closest node to recurse into */
1826                 if (data->closest_axis[node->main_axis]) {
1827                         for (int i = 0; i != node->totnode; i++) {
1828                                 const float *bv = node->children[i]->bv;
1829                                 const float bb_min[3] = {bv[0], bv[2], bv[4]};
1830                                 const float bb_max[3] = {bv[1], bv[3], bv[5]};
1831
1832                                 int isect_type = isect_aabb_planes_v3(data->clip_plane, data->clip_plane_len, bb_min, bb_max);
1833
1834                                 if ((isect_type != ISECT_AABB_PLANE_BEHIND_ANY) && dist_squared_to_projected_aabb(
1835                                         &data->precalc, bb_min, bb_max,
1836                                         data->closest_axis) <= data->nearest.dist_sq)
1837                                 {
1838                                         if (isect_type == ISECT_AABB_PLANE_CROSS_ANY) {
1839                                                 bvhtree_nearest_projected_with_clipplane_test_dfs_recursive(data, node->children[i]);
1840                                         }
1841                                         else {
1842                                                 /* ISECT_AABB_PLANE_IN_FRONT_ALL */
1843                                                 bvhtree_nearest_projected_dfs_recursive(data, node->children[i]);
1844                                         }
1845                                 }
1846                         }
1847                 }
1848                 else {
1849                         for (int i = node->totnode; i--;) {
1850                                 const float *bv = node->children[i]->bv;
1851                                 const float bb_min[3] = {bv[0], bv[2], bv[4]};
1852                                 const float bb_max[3] = {bv[1], bv[3], bv[5]};
1853
1854                                 int isect_type = isect_aabb_planes_v3(data->clip_plane, data->clip_plane_len, bb_min, bb_max);
1855
1856                                 if (isect_type != ISECT_AABB_PLANE_BEHIND_ANY && dist_squared_to_projected_aabb(
1857                                         &data->precalc, bb_min, bb_max,
1858                                         data->closest_axis) <= data->nearest.dist_sq)
1859                                 {
1860                                         if (isect_type == ISECT_AABB_PLANE_CROSS_ANY) {
1861                                                 bvhtree_nearest_projected_with_clipplane_test_dfs_recursive(data, node->children[i]);
1862                                         }
1863                                         else {
1864                                                 /* ISECT_AABB_PLANE_IN_FRONT_ALL */
1865                                                 bvhtree_nearest_projected_dfs_recursive(data, node->children[i]);
1866                                         }
1867                                 }
1868                         }
1869                 }
1870         }
1871 }
1872
1873 int BLI_bvhtree_find_nearest_projected(
1874         BVHTree *tree, float projmat[4][4], float winsize[2], float mval[2],
1875         float clip_plane[6][4], int clip_plane_len,
1876         BVHTreeNearest *nearest,
1877         BVHTree_NearestProjectedCallback callback, void *userdata)
1878 {
1879         BVHNode *root = tree->nodes[tree->totleaf];
1880         if (root != NULL) {
1881                 BVHNearestProjectedData data;
1882                 dist_squared_to_projected_aabb_precalc(
1883                         &data.precalc, projmat, winsize, mval);
1884
1885                 data.callback = callback;
1886                 data.userdata = userdata;
1887
1888                 if (clip_plane) {
1889                         data.clip_plane_len = clip_plane_len;
1890                         for (int i = 0; i < data.clip_plane_len; i++) {
1891                                 copy_v4_v4(data.clip_plane[i], clip_plane[i]);
1892                         }
1893                 }
1894                 else {
1895                         data.clip_plane_len = 1;
1896                         planes_from_projmat(
1897                                 projmat,
1898                                 NULL, NULL, NULL, NULL,
1899                                 data.clip_plane[0], NULL);
1900                 }
1901
1902                 if (nearest) {
1903                         memcpy(&data.nearest, nearest, sizeof(*nearest));
1904                 }
1905                 else {
1906                         data.nearest.index = -1;
1907                         data.nearest.dist_sq = FLT_MAX;
1908                 }
1909                 {
1910                         const float bb_min[3] = {root->bv[0], root->bv[2], root->bv[4]};
1911                         const float bb_max[3] = {root->bv[1], root->bv[3], root->bv[5]};
1912
1913                         int isect_type = isect_aabb_planes_v3(data.clip_plane, data.clip_plane_len, bb_min, bb_max);
1914
1915                         if (isect_type != 0 && dist_squared_to_projected_aabb(
1916                                 &data.precalc, bb_min, bb_max,
1917                                 data.closest_axis) <= data.nearest.dist_sq)
1918                         {
1919                                 if (isect_type == 1) {
1920                                         bvhtree_nearest_projected_with_clipplane_test_dfs_recursive(&data, root);
1921                                 }
1922                                 else {
1923                                         bvhtree_nearest_projected_dfs_recursive(&data, root);
1924                                 }
1925                         }
1926                 }
1927
1928                 if (nearest) {
1929                         memcpy(nearest, &data.nearest, sizeof(*nearest));
1930                 }
1931
1932                 return data.nearest.index;
1933         }
1934         return -1;
1935 }
1936
1937 /** \} */
1938
1939
1940 /* -------------------------------------------------------------------- */
1941
1942 /** \name BLI_bvhtree_walk_dfs
1943  * \{ */
1944
1945 typedef struct BVHTree_WalkData {
1946         BVHTree_WalkParentCallback walk_parent_cb;
1947         BVHTree_WalkLeafCallback walk_leaf_cb;
1948         BVHTree_WalkOrderCallback walk_order_cb;
1949         void *userdata;
1950 } BVHTree_WalkData;
1951
1952 /**
1953  * Runs first among nodes children of the first node before going to the next node in the same layer.
1954  *
1955  * \return false to break out of the search early.
1956  */
1957 static bool bvhtree_walk_dfs_recursive(
1958         BVHTree_WalkData *walk_data,
1959         const BVHNode *node)
1960 {
1961         if (node->totnode == 0) {
1962                 return walk_data->walk_leaf_cb((const BVHTreeAxisRange *)node->bv, node->index, walk_data->userdata);
1963         }
1964         else {
1965                 /* First pick the closest node to recurse into */
1966                 if (walk_data->walk_order_cb((const BVHTreeAxisRange *)node->bv, node->main_axis, walk_data->userdata)) {
1967                         for (int i = 0; i != node->totnode; i++) {
1968                                 if (walk_data->walk_parent_cb((const BVHTreeAxisRange *)node->children[i]->bv, walk_data->userdata)) {
1969                                         if (!bvhtree_walk_dfs_recursive(walk_data, node->children[i])) {
1970                                                 return false;
1971                                         }
1972                                 }
1973                         }
1974                 }
1975                 else {
1976                         for (int i = node->totnode - 1; i >= 0; i--) {
1977                                 if (walk_data->walk_parent_cb((const BVHTreeAxisRange *)node->children[i]->bv, walk_data->userdata)) {
1978                                         if (!bvhtree_walk_dfs_recursive(walk_data, node->children[i])) {
1979                                                 return false;
1980                                         }
1981                                 }
1982                         }
1983                 }
1984         }
1985         return true;
1986 }
1987
1988 /**
1989  * This is a generic function to perform a depth first search on the BVHTree
1990  * where the search order and nodes traversed depend on callbacks passed in.
1991  *
1992  * \param tree: Tree to walk.
1993  * \param walk_parent_cb: Callback on a parents bound-box to test if it should be traversed.
1994  * \param walk_leaf_cb: Callback to test leaf nodes, callback must store its own result,
1995  * returning false exits early.
1996  * \param walk_order_cb: Callback that indicates which direction to search,
1997  * either from the node with the lower or higher k-dop axis value.
1998  * \param userdata: Argument passed to all callbacks.
1999  */
2000 void BLI_bvhtree_walk_dfs(
2001         BVHTree *tree,
2002         BVHTree_WalkParentCallback walk_parent_cb,
2003         BVHTree_WalkLeafCallback walk_leaf_cb,
2004         BVHTree_WalkOrderCallback walk_order_cb, void *userdata)
2005 {
2006         const BVHNode *root = tree->nodes[tree->totleaf];
2007         if (root != NULL) {
2008                 BVHTree_WalkData walk_data = {walk_parent_cb, walk_leaf_cb, walk_order_cb, userdata};
2009                 /* first make sure the bv of root passes in the test too */
2010                 if (walk_parent_cb((const BVHTreeAxisRange *)root->bv, userdata)) {
2011                         bvhtree_walk_dfs_recursive(&walk_data, root);
2012                 }
2013         }
2014 }
2015
2016 /** \} */