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