b3f11039ce15bfc5491b32e1d5d6f0a262f6e90b
[blender.git] / source / blender / blenlib / intern / BLI_kdopbvh.c
1 /*  kdop.c      
2
3 *
4 * ***** BEGIN GPL/BL DUAL LICENSE BLOCK *****
5 *
6 * This program is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU General Public License
8 * as published by the Free Software Foundation; either version 2
9 * of the License, or (at your option) any later version. The Blender
10 * Foundation also sells licenses for use in proprietary software under
11 * the Blender License.  See http://www.blender.org/BL/ for information
12 * about this.
13 *
14 * This program is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17 * GNU General Public License for more details.
18 *
19 * You should have received a copy of the GNU General Public License
20 * along with this program; if not, write to the Free Software Foundation,
21 * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
22 *
23 * The Original Code is Copyright (C) Blender Foundation
24 * All rights reserved.
25 *
26 * The Original Code is: all of this file.
27 *
28 * Contributor(s): Daniel Genrich, Andre Pinto
29 *
30 * ***** END GPL/BL DUAL LICENSE BLOCK *****
31 */
32
33 #include "math.h"
34 #include <stdio.h>
35 #include <stdlib.h> 
36 #include <string.h>
37
38 #include "MEM_guardedalloc.h"
39
40 #include "BKE_utildefines.h"
41
42 #include "BLI_kdopbvh.h"
43 #include "BLI_arithb.h"
44
45 #ifdef _OPENMP
46 #include <omp.h>
47 #endif
48
49 #include <time.h>
50
51 /* Util macros */
52 #define TO_STR(a)       #a
53 #define JOIN(a,b)       a##b
54
55 /* Benchmark macros */
56 #if 1
57
58 #define BENCH(a)        \
59         do {                    \
60                 clock_t _clock_init = clock();  \
61                 (a);                                                    \
62                 printf("%s: %fms\n", #a, (float)(clock()-_clock_init)*1000/CLOCKS_PER_SEC);     \
63 } while(0)
64
65 #define BENCH_VAR(name)         clock_t JOIN(_bench_step,name) = 0, JOIN(_bench_total,name) = 0
66 #define BENCH_BEGIN(name)       JOIN(_bench_step, name) = clock()
67 #define BENCH_END(name)         JOIN(_bench_total,name) += clock() - JOIN(_bench_step,name)
68 #define BENCH_RESET(name)       JOIN(_bench_total, name) = 0
69 #define BENCH_REPORT(name)      printf("%s: %fms\n", TO_STR(name), JOIN(_bench_total,name)*1000.0f/CLOCKS_PER_SEC)
70
71 #else
72
73 #define BENCH(a)        (a)
74 #define BENCH_VAR(name)
75 #define BENCH_BEGIN(name)
76 #define BENCH_END(name)
77 #define BENCH_RESET(name)
78 #define BENCH_REPORT(name)
79
80 #endif
81
82
83
84 typedef struct BVHNode
85 {
86         struct BVHNode *children[8]; // max 8 children
87         struct BVHNode *parent; // needed for bottom - top update
88         float bv[26]; // Bounding volume of all nodes, max 13 axis
89         int index; /* face, edge, vertex index */
90         char totnode; // how many nodes are used, used for speedup
91         char traversed;  // how many nodes already traversed until this level?
92         char main_axis;
93 } BVHNode;
94
95 struct BVHTree
96 {
97         BVHNode **nodes;
98         BVHNode *nodearray; /* pre-alloc branch nodes */
99         float   epsilon; /* epsilon is used for inflation of the k-dop     */
100         int     totleaf; // leafs
101         int     totbranch;
102         char    tree_type; // type of tree (4 => quadtree)
103         char    axis; // kdop type (6 => OBB, 7 => AABB, ...)
104         char    start_axis, stop_axis; // KDOP_AXES array indices according to axis
105 };
106
107 typedef struct BVHOverlapData 
108 {  
109         BVHTree *tree1, *tree2; 
110         BVHTreeOverlap *overlap; 
111         int i, max_overlap; /* i is number of overlaps */
112 } BVHOverlapData;
113 ////////////////////////////////////////
114
115
116 ////////////////////////////////////////////////////////////////////////
117 // Bounding Volume Hierarchy Definition
118 // 
119 // Notes: From OBB until 26-DOP --> all bounding volumes possible, just choose type below
120 // Notes: You have to choose the type at compile time ITM
121 // Notes: You can choose the tree type --> binary, quad, octree, choose below
122 ////////////////////////////////////////////////////////////////////////
123
124 static float KDOP_AXES[13][3] =
125 { {1.0, 0, 0}, {0, 1.0, 0}, {0, 0, 1.0}, {1.0, 1.0, 1.0}, {1.0, -1.0, 1.0}, {1.0, 1.0, -1.0},
126 {1.0, -1.0, -1.0}, {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},
127 {0, 1.0, -1.0}
128 };
129
130 //////////////////////////////////////////////////////////////////////////////////////////////////////
131 // Introsort 
132 // with permission deriven from the following Java code:
133 // http://ralphunden.net/content/tutorials/a-guide-to-introsort/
134 // and he derived it from the SUN STL 
135 //////////////////////////////////////////////////////////////////////////////////////////////////////
136 static int size_threshold = 16;
137 /*
138 * Common methods for all algorithms
139 */
140 static int floor_lg(int a)
141 {
142         return (int)(floor(log(a)/log(2)));
143 }
144
145 /*
146 * Insertion sort algorithm
147 */
148 static void bvh_insertionsort(BVHNode **a, int lo, int hi, int axis)
149 {
150         int i,j;
151         BVHNode *t;
152         for (i=lo; i < hi; i++)
153         {
154                 j=i;
155                 t = a[i];
156                 while((j!=lo) && (t->bv[axis] < (a[j-1])->bv[axis]))
157                 {
158                         a[j] = a[j-1];
159                         j--;
160                 }
161                 a[j] = t;
162         }
163 }
164
165 static int bvh_partition(BVHNode **a, int lo, int hi, BVHNode * x, int axis)
166 {
167         int i=lo, j=hi;
168         while (1)
169         {
170                 while ((a[i])->bv[axis] < x->bv[axis]) i++;
171                 j--;
172                 while (x->bv[axis] < (a[j])->bv[axis]) j--;
173                 if(!(i < j))
174                         return i;
175                 SWAP( BVHNode* , a[i], a[j]);
176                 i++;
177         }
178 }
179
180 /*
181 * Heapsort algorithm
182 */
183 static void bvh_downheap(BVHNode **a, int i, int n, int lo, int axis)
184 {
185         BVHNode * d = a[lo+i-1];
186         int child;
187         while (i<=n/2)
188         {
189                 child = 2*i;
190                 if ((child < n) && ((a[lo+child-1])->bv[axis] < (a[lo+child])->bv[axis]))
191                 {
192                         child++;
193                 }
194                 if (!(d->bv[axis] < (a[lo+child-1])->bv[axis])) break;
195                 a[lo+i-1] = a[lo+child-1];
196                 i = child;
197         }
198         a[lo+i-1] = d;
199 }
200
201 static void bvh_heapsort(BVHNode **a, int lo, int hi, int axis)
202 {
203         int n = hi-lo, i;
204         for (i=n/2; i>=1; i=i-1)
205         {
206                 bvh_downheap(a, i,n,lo, axis);
207         }
208         for (i=n; i>1; i=i-1)
209         {
210                 SWAP(BVHNode*, a[lo],a[lo+i-1]);
211                 bvh_downheap(a, 1,i-1,lo, axis);
212         }
213 }
214
215 static BVHNode *bvh_medianof3(BVHNode **a, int lo, int mid, int hi, int axis) // returns Sortable
216 {
217         if ((a[mid])->bv[axis] < (a[lo])->bv[axis])
218         {
219                 if ((a[hi])->bv[axis] < (a[mid])->bv[axis])
220                         return a[mid];
221                 else
222                 {
223                         if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
224                                 return a[hi];
225                         else
226                                 return a[lo];
227                 }
228         }
229         else
230         {
231                 if ((a[hi])->bv[axis] < (a[mid])->bv[axis])
232                 {
233                         if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
234                                 return a[lo];
235                         else
236                                 return a[hi];
237                 }
238                 else
239                         return a[mid];
240         }
241 }
242 /*
243 * Quicksort algorithm modified for Introsort
244 */
245 static void bvh_introsort_loop (BVHNode **a, int lo, int hi, int depth_limit, int axis)
246 {
247         int p;
248
249         while (hi-lo > size_threshold)
250         {
251                 if (depth_limit == 0)
252                 {
253                         bvh_heapsort(a, lo, hi, axis);
254                         return;
255                 }
256                 depth_limit=depth_limit-1;
257                 p=bvh_partition(a, lo, hi, bvh_medianof3(a, lo, lo+((hi-lo)/2)+1, hi-1, axis), axis);
258                 bvh_introsort_loop(a, p, hi, depth_limit, axis);
259                 hi=p;
260         }
261 }
262
263 static void sort(BVHNode **a0, int begin, int end, int axis)
264 {
265         if (begin < end)
266         {
267                 BVHNode **a=a0;
268                 bvh_introsort_loop(a, begin, end, 2*floor_lg(end-begin), axis);
269                 bvh_insertionsort(a, begin, end, axis);
270         }
271 }
272 void sort_along_axis(BVHTree *tree, int start, int end, int axis)
273 {
274         sort(tree->nodes, start, end, axis);
275 }
276
277 //after a call to this function you can expect one of:
278 //      every node to left of a[n] are smaller than it
279 //      every node to the right of a[n-1] are greater than it
280 void partition_nth_element(BVHNode **a, int _begin, int _end, int n, int axis)
281 {
282         int begin = _begin, end = _end;
283         while(begin < n && end >= n)
284         {
285                 int mid = bvh_partition(a, begin, end, bvh_medianof3(a, begin, (begin+end-1)/2, end-1, axis), axis );
286
287                 if(mid >= n)
288                         end = n-1;
289                 else
290                         begin = n+1;
291         }
292
293 }
294
295
296 //////////////////////////////////////////////////////////////////////////////////////////////////////
297
298 void BLI_bvhtree_free(BVHTree *tree)
299 {       
300         if(tree)
301         {
302                 MEM_freeN(tree->nodes);
303                 MEM_freeN(tree->nodearray);
304                 MEM_freeN(tree);
305         }
306 }
307
308 BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis)
309 {
310         BVHTree *tree;
311         int numbranches=0, i;
312         
313         // only support up to octree
314         if(tree_type > 8)
315                 return NULL;
316
317         tree = (BVHTree *)MEM_callocN(sizeof(BVHTree), "BVHTree");
318         
319         if(tree)
320         {
321                 // calculate max number of branches, our bvh kdop is "almost perfect"
322                 for(i = 1; i <= (int)ceil((float)((float)log(maxsize)/(float)log(tree_type))); i++)
323                         numbranches += (pow(tree_type, i) / tree_type);
324                 
325                 tree->nodes = (BVHNode **)MEM_callocN(sizeof(BVHNode *)*(numbranches+maxsize + tree_type), "BVHNodes");
326                 
327                 if(!tree->nodes)
328                 {
329                         MEM_freeN(tree);
330                         return NULL;
331                 }
332                 
333                 tree->nodearray = (BVHNode *)MEM_callocN(sizeof(BVHNode)*(numbranches+maxsize + tree_type), "BVHNodeArray");
334                 
335                 if(!tree->nodearray)
336                 {
337                         MEM_freeN(tree);
338                         MEM_freeN(tree->nodes);
339                         return NULL;
340                 }
341                 
342                 tree->epsilon = epsilon;
343                 tree->tree_type = tree_type; 
344                 tree->axis = axis;
345                 
346                 if(axis == 26)
347                 {
348                         tree->start_axis = 0;
349                         tree->stop_axis = 13;
350                 }
351                 else if(axis == 18)
352                 {
353                         tree->start_axis = 7;
354                         tree->stop_axis = 13;
355                 }
356                 else if(axis == 14)
357                 {
358                         tree->start_axis = 0;
359                         tree->stop_axis = 7;
360                 }
361                 else if(axis == 8) // AABB
362                 {
363                         tree->start_axis = 0;
364                         tree->stop_axis = 4;
365                 }
366                 else if(axis == 6) // OBB
367                 {
368                         tree->start_axis = 0;
369                         tree->stop_axis = 3;
370                 }
371                 else
372                 {
373                         BLI_bvhtree_free(tree);
374                         return NULL;
375                 }
376         }
377
378         return tree;
379 }
380
381 static void create_kdop_hull(BVHTree *tree, BVHNode *node, float *co, int numpoints, int moving)
382 {
383         float newminmax;
384         int i, k;
385         
386         // don't init boudings for the moving case
387         if(!moving)
388         {
389                 for (i = tree->start_axis; i < tree->stop_axis; i++)
390                 {
391                         node->bv[2*i] = FLT_MAX;
392                         node->bv[2*i + 1] = -FLT_MAX;
393                 }
394         }
395         
396         for(k = 0; k < numpoints; k++)
397         {
398                 // for all Axes.
399                 for (i = tree->start_axis; i < tree->stop_axis; i++)
400                 {
401                         newminmax = INPR(&co[k * 3], KDOP_AXES[i]);
402                         if (newminmax < node->bv[2 * i])
403                                 node->bv[2 * i] = newminmax;
404                         if (newminmax > node->bv[(2 * i) + 1])
405                                 node->bv[(2 * i) + 1] = newminmax;
406                 }
407         }
408 }
409
410 // depends on the fact that the BVH's for each face is already build
411 static void refit_kdop_hull(BVHTree *tree, BVHNode *node, int start, int end)
412 {
413         float newmin,newmax;
414         int i, j;
415         float *bv = node->bv;
416         
417         for (i = tree->start_axis; i < tree->stop_axis; i++)
418         {
419                 bv[2*i] = FLT_MAX;
420                 bv[2*i + 1] = -FLT_MAX;
421         }
422
423         for (j = start; j < end; j++)
424         {
425                 // for all Axes.
426                 for (i = tree->start_axis; i < tree->stop_axis; i++)
427                 {
428                         newmin = tree->nodes[j]->bv[(2 * i)];   
429                         if ((newmin < bv[(2 * i)]))
430                                 bv[(2 * i)] = newmin;
431  
432                         newmax = tree->nodes[j]->bv[(2 * i) + 1];
433                         if ((newmax > bv[(2 * i) + 1]))
434                                 bv[(2 * i) + 1] = newmax;
435                 }
436         }
437 }
438
439 int BLI_bvhtree_insert(BVHTree *tree, int index, float *co, int numpoints)
440 {
441         BVHNode *node= NULL;
442         int i;
443         
444         // insert should only possible as long as tree->totbranch is 0
445         if(tree->totbranch > 0)
446                 return 0;
447         
448         if(tree->totleaf+1 >= MEM_allocN_len(tree->nodes))
449                 return 0;
450         
451         // TODO check if have enough nodes in array
452         
453         node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]);
454         tree->totleaf++;
455         
456         create_kdop_hull(tree, node, co, numpoints, 0);
457         
458         // inflate the bv with some epsilon
459         for (i = tree->start_axis; i < tree->stop_axis; i++)
460         {
461                 node->bv[(2 * i)] -= tree->epsilon; // minimum 
462                 node->bv[(2 * i) + 1] += tree->epsilon; // maximum 
463         }
464
465         node->index= index;
466         
467         return 1;
468 }
469
470 // only supports x,y,z axis in the moment
471 // but we should use a plain and simple function here for speed sake
472 static char get_largest_axis(float *bv)
473 {
474         float middle_point[3];
475
476         middle_point[0] = (bv[1]) - (bv[0]); // x axis
477         middle_point[1] = (bv[3]) - (bv[2]); // y axis
478         middle_point[2] = (bv[5]) - (bv[4]); // z axis
479         if (middle_point[0] > middle_point[1]) 
480         {
481                 if (middle_point[0] > middle_point[2])
482                         return 1; // max x axis
483                 else
484                         return 5; // max z axis
485         }
486         else 
487         {
488                 if (middle_point[1] > middle_point[2])
489                         return 3; // max y axis
490                 else
491                         return 5; // max z axis
492         }
493 }
494
495 static void bvh_div_nodes(BVHTree *tree, BVHNode *node, int start, int end, char lastaxis)
496 {
497         char laxis;
498         int i, tend;
499         BVHNode *tnode;
500         int slice = (end-start+tree->tree_type-1)/tree->tree_type;      //division rounded up
501         
502         // Determine which axis to split along
503         laxis = get_largest_axis(node->bv);
504         
505         // split nodes along longest axis
506         for (i=0; start < end; start += slice, i++) //i counts the current child
507         {       
508                 tend = start + slice;
509                 
510                 partition_nth_element(tree->nodes, start, end, tend, laxis);
511                 
512                 if(tend > end) tend = end;
513                 
514                 if(tend-start == 1)     // ok, we have 1 left for this node
515                 {
516                         node->children[i] = tree->nodes[start];
517                         node->children[i]->parent = node;
518                 }
519                 else
520                 {
521                         tnode = node->children[i] = tree->nodes[tree->totleaf  + tree->totbranch] = &(tree->nodearray[tree->totbranch + tree->totleaf]);
522                         tree->totbranch++;
523                         tnode->parent = node;
524
525                         refit_kdop_hull(tree, tnode, start, tend);
526                         bvh_div_nodes(tree, tnode, start, tend, laxis);
527                 }
528                 node->totnode++;
529         }
530         
531         return;
532 }
533
534 static void verify_tree(BVHTree *tree)
535 {
536         int i, j, check = 0;
537         
538         // check the pointer list
539         for(i = 0; i < tree->totleaf; i++)
540         {
541                 if(tree->nodes[i]->parent == NULL)
542                         printf("Leaf has no parent: %d\n", i);
543                 else
544                 {
545                         for(j = 0; j < tree->tree_type; j++)
546                         {
547                                 if(tree->nodes[i]->parent->children[j] == tree->nodes[i])
548                                         check = 1;
549                         }
550                         if(!check)
551                         {
552                                 printf("Parent child relationship doesn't match: %d\n", i);
553                         }
554                         check = 0;
555                 }
556         }
557         
558         // check the leaf list
559         for(i = 0; i < tree->totleaf; i++)
560         {
561                 if(tree->nodearray[i].parent == NULL)
562                         printf("Leaf has no parent: %d\n", i);
563                 else
564                 {
565                         for(j = 0; j < tree->tree_type; j++)
566                         {
567                                 if(tree->nodearray[i].parent->children[j] == &tree->nodearray[i])
568                                         check = 1;
569                         }
570                         if(!check)
571                         {
572                                 printf("Parent child relationship doesn't match: %d\n", i);
573                         }
574                         check = 0;
575                 }
576         }
577         
578         printf("branches: %d, leafs: %d, total: %d\n", tree->totbranch, tree->totleaf, tree->totbranch + tree->totleaf);
579 }
580         
581 void BLI_bvhtree_balance(BVHTree *tree)
582 {
583         BVHNode *node;
584         
585         if(tree->totleaf == 0)
586                 return;
587         
588         // create root node
589         node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]);
590         tree->totbranch++;
591         
592         // refit root bvh node
593         refit_kdop_hull(tree, tree->nodes[tree->totleaf], 0, tree->totleaf);
594         // create + balance tree
595         bvh_div_nodes(tree, tree->nodes[tree->totleaf], 0, tree->totleaf, 0);
596         
597         // verify_tree(tree);
598 }
599
600 // overlap - is it possbile for 2 bv's to collide ?
601 static int tree_overlap(float *bv1, float *bv2, int start_axis, int stop_axis)
602 {
603         float *bv1_end = bv1 + (stop_axis<<1);
604                 
605         bv1 += start_axis<<1;
606         bv2 += start_axis<<1;
607         
608         // test all axis if min + max overlap
609         for (; bv1 != bv1_end; bv1+=2, bv2+=2)
610         {
611                 if ((*(bv1) > *(bv2 + 1)) || (*(bv2) > *(bv1 + 1))) 
612                         return 0;
613         }
614         
615         return 1;
616 }
617
618 static void traverse(BVHOverlapData *data, BVHNode *node1, BVHNode *node2)
619 {
620         int j;
621         
622         if(tree_overlap(node1->bv, node2->bv, MIN2(data->tree1->start_axis, data->tree2->start_axis), MIN2(data->tree1->stop_axis, data->tree2->stop_axis)))
623         {
624                 // check if node1 is a leaf
625                 if(node1->index)
626                 {
627                         // check if node2 is a leaf
628                         if(node2->index)
629                         {
630                                 
631                                 if(node1 == node2)
632                                 {
633                                         return;
634                                 }
635                                         
636                                 if(data->i >= data->max_overlap)
637                                 {       
638                                         // try to make alloc'ed memory bigger
639                                         data->overlap = realloc(data->overlap, sizeof(BVHTreeOverlap)*data->max_overlap*2);
640                                         
641                                         if(!data->overlap)
642                                         {
643                                                 printf("Out of Memory in traverse\n");
644                                                 return;
645                                         }
646                                         data->max_overlap *= 2;
647                                 }
648                                 
649                                 // both leafs, insert overlap!
650                                 data->overlap[data->i].indexA = node1->index;
651                                 data->overlap[data->i].indexB = node2->index;
652
653                                 data->i++;
654                         }
655                         else
656                         {
657                                 for(j = 0; j < data->tree2->tree_type; j++)
658                                 {
659                                         if(node2->children[j])
660                                                 traverse(data, node1, node2->children[j]);
661                                 }
662                         }
663                 }
664                 else
665                 {
666                         
667                         for(j = 0; j < data->tree2->tree_type; j++)
668                         {
669                                 if(node1->children[j])
670                                         traverse(data, node1->children[j], node2);
671                         }
672                 }
673         }
674         return;
675 }
676
677 BVHTreeOverlap *BLI_bvhtree_overlap(BVHTree *tree1, BVHTree *tree2, int *result)
678 {
679         int j, total = 0;
680         BVHTreeOverlap *overlap = NULL, *to = NULL;
681         BVHOverlapData *data[tree1->tree_type];
682         
683         // check for compatibility of both trees (can't compare 14-DOP with 18-DOP)
684         if((tree1->axis != tree2->axis) && ((tree1->axis == 14) || tree2->axis == 14))
685                 return 0;
686         
687         // fast check root nodes for collision before doing big splitting + traversal
688         if(!tree_overlap(tree1->nodes[tree1->totleaf]->bv, tree2->nodes[tree2->totleaf]->bv, MIN2(tree1->start_axis, tree2->start_axis), MIN2(tree1->stop_axis, tree2->stop_axis)))
689                 return 0;
690         
691         for(j = 0; j < tree1->tree_type; j++)
692         {
693                 data[j] = (BVHOverlapData *)MEM_callocN(sizeof(BVHOverlapData), "BVHOverlapData");
694                 
695                 // init BVHOverlapData
696                 data[j]->overlap = (BVHTreeOverlap *)malloc(sizeof(BVHTreeOverlap)*MAX2(tree1->totleaf, tree2->totleaf));
697                 data[j]->tree1 = tree1;
698                 data[j]->tree2 = tree2;
699                 data[j]->max_overlap = MAX2(tree1->totleaf, tree2->totleaf);
700                 data[j]->i = 0;
701         }
702
703 #pragma omp parallel for private(j) schedule(static)
704         for(j = 0; j < tree1->tree_type; j++)
705         {
706                 traverse(data[j], tree1->nodes[tree1->totleaf]->children[j], tree2->nodes[tree2->totleaf]);
707         }
708         
709         for(j = 0; j < tree1->tree_type; j++)
710                 total += data[j]->i;
711         
712         to = overlap = (BVHTreeOverlap *)MEM_callocN(sizeof(BVHTreeOverlap)*total, "BVHTreeOverlap");
713         
714         for(j = 0; j < tree1->tree_type; j++)
715         {
716                 memcpy(to, data[j]->overlap, data[j]->i*sizeof(BVHTreeOverlap));
717                 to+=data[j]->i;
718         }
719         
720         for(j = 0; j < tree1->tree_type; j++)
721         {
722                 free(data[j]->overlap);
723                 MEM_freeN(data[j]);
724         }
725         
726         (*result) = total;
727         return overlap;
728 }
729
730
731 // bottom up update of bvh tree:
732 // join the 4 children here
733 static void node_join(BVHTree *tree, BVHNode *node)
734 {
735         int i, j;
736         
737         for (i = tree->start_axis; i < tree->stop_axis; i++)
738         {
739                 node->bv[2*i] = FLT_MAX;
740                 node->bv[2*i + 1] = -FLT_MAX;
741         }
742         
743         for (i = 0; i < tree->tree_type; i++)
744         {
745                 if (node->children[i]) 
746                 {
747                         for (j = tree->start_axis; j < tree->stop_axis; j++)
748                         {
749                                 // update minimum 
750                                 if (node->children[i]->bv[(2 * j)] < node->bv[(2 * j)]) 
751                                         node->bv[(2 * j)] = node->children[i]->bv[(2 * j)];
752                                 
753                                 // update maximum 
754                                 if (node->children[i]->bv[(2 * j) + 1] > node->bv[(2 * j) + 1])
755                                         node->bv[(2 * j) + 1] = node->children[i]->bv[(2 * j) + 1];
756                         }
757                 }
758                 else
759                         break;
760         }
761 }
762
763 // call before BLI_bvhtree_update_tree()
764 int BLI_bvhtree_update_node(BVHTree *tree, int index, float *co, float *co_moving, int numpoints)
765 {
766         BVHNode *node= NULL;
767         int i = 0;
768         
769         // check if index exists
770         if(index > tree->totleaf)
771                 return 0;
772         
773         node = tree->nodearray + index;
774         
775         create_kdop_hull(tree, node, co, numpoints, 0);
776         
777         if(co_moving)
778                 create_kdop_hull(tree, node, co_moving, numpoints, 1);
779         
780         // inflate the bv with some epsilon
781         for (i = tree->start_axis; i < tree->stop_axis; i++)
782         {
783                 node->bv[(2 * i)] -= tree->epsilon; // minimum 
784                 node->bv[(2 * i) + 1] += tree->epsilon; // maximum 
785         }
786         
787         return 1;
788 }
789
790 // call BLI_bvhtree_update_node() first for every node/point/triangle
791 void BLI_bvhtree_update_tree(BVHTree *tree)
792 {
793         BVHNode *leaf, *parent;
794         
795         // reset tree traversing flag
796         for (leaf = tree->nodearray + tree->totleaf; leaf != tree->nodearray + tree->totleaf + tree->totbranch; leaf++)
797                 leaf->traversed = 0;
798         
799         for (leaf = tree->nodearray; leaf != tree->nodearray + tree->totleaf; leaf++)
800         {
801                 for (parent = leaf->parent; parent; parent = parent->parent)
802                 {
803                         parent->traversed++;    // we tried to go up in hierarchy 
804                         if (parent->traversed < parent->totnode) 
805                                 break;  // we do not need to check further 
806                         else 
807                                 node_join(tree, parent);
808                 }
809         }
810 }
811
812 float BLI_bvhtree_getepsilon(BVHTree *tree)
813 {
814         return tree->epsilon;
815 }