Merged 15170:15635 from trunk (no conflicts or even merges)
[blender.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, char lastaxis)
488 {
489         char laxis;
490         int i, tend;
491         BVHNode *tnode;
492         int slice = (end-start+tree->tree_type-1)/tree->tree_type;      //division rounded up
493         
494         // Determine which axis to split along
495         laxis = get_largest_axis(node->bv);
496         
497         // split nodes along longest axis
498         for (i=0; start < end; start += slice, i++) //i counts the current child
499         {       
500                 tend = start + slice;
501                 
502                 if(tend > end) tend = end;
503                 
504                 if(tend-start == 1)     // ok, we have 1 left for this node
505                 {
506                         node->children[i] = tree->nodes[start];
507                         node->children[i]->parent = node;
508                 }
509                 else
510                 {
511                         tnode = node->children[i] = tree->nodes[tree->totleaf  + tree->totbranch] = &(tree->nodearray[tree->totbranch + tree->totleaf]);
512                         tree->totbranch++;
513                         tnode->parent = node;
514                         
515                         if(tend != end)
516                                 partition_nth_element(tree->nodes, start, end, tend, laxis);
517                         refit_kdop_hull(tree, tnode, start, tend);
518                         bvh_div_nodes(tree, tnode, start, tend, laxis);
519                 }
520                 node->totnode++;
521         }
522         
523         return;
524 }
525
526 #if 0
527 static void verify_tree(BVHTree *tree)
528 {
529         int i, j, check = 0;
530         
531         // check the pointer list
532         for(i = 0; i < tree->totleaf; i++)
533         {
534                 if(tree->nodes[i]->parent == NULL)
535                         printf("Leaf has no parent: %d\n", i);
536                 else
537                 {
538                         for(j = 0; j < tree->tree_type; j++)
539                         {
540                                 if(tree->nodes[i]->parent->children[j] == tree->nodes[i])
541                                         check = 1;
542                         }
543                         if(!check)
544                         {
545                                 printf("Parent child relationship doesn't match: %d\n", i);
546                         }
547                         check = 0;
548                 }
549         }
550         
551         // check the leaf list
552         for(i = 0; i < tree->totleaf; i++)
553         {
554                 if(tree->nodearray[i].parent == NULL)
555                         printf("Leaf has no parent: %d\n", i);
556                 else
557                 {
558                         for(j = 0; j < tree->tree_type; j++)
559                         {
560                                 if(tree->nodearray[i].parent->children[j] == &tree->nodearray[i])
561                                         check = 1;
562                         }
563                         if(!check)
564                         {
565                                 printf("Parent child relationship doesn't match: %d\n", i);
566                         }
567                         check = 0;
568                 }
569         }
570         
571         printf("branches: %d, leafs: %d, total: %d\n", tree->totbranch, tree->totleaf, tree->totbranch + tree->totleaf);
572 }
573 #endif
574         
575 void BLI_bvhtree_balance(BVHTree *tree)
576 {
577         BVHNode *node;
578         
579         if(tree->totleaf == 0)
580                 return;
581         
582         // create root node
583         node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]);
584         tree->totbranch++;
585         
586         // refit root bvh node
587         refit_kdop_hull(tree, tree->nodes[tree->totleaf], 0, tree->totleaf);
588         // create + balance tree
589         bvh_div_nodes(tree, tree->nodes[tree->totleaf], 0, tree->totleaf, 0);
590         
591         // verify_tree(tree);
592 }
593
594 // overlap - is it possbile for 2 bv's to collide ?
595 static int tree_overlap(float *bv1, float *bv2, int start_axis, int stop_axis)
596 {
597         float *bv1_end = bv1 + (stop_axis<<1);
598                 
599         bv1 += start_axis<<1;
600         bv2 += start_axis<<1;
601         
602         // test all axis if min + max overlap
603         for (; bv1 != bv1_end; bv1+=2, bv2+=2)
604         {
605                 if ((*(bv1) > *(bv2 + 1)) || (*(bv2) > *(bv1 + 1))) 
606                         return 0;
607         }
608         
609         return 1;
610 }
611
612 static void traverse(BVHOverlapData *data, BVHNode *node1, BVHNode *node2)
613 {
614         int j;
615         
616         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)))
617         {
618                 // check if node1 is a leaf
619                 if(!node1->totnode)
620                 {
621                         // check if node2 is a leaf
622                         if(!node2->totnode)
623                         {
624                                 
625                                 if(node1 == node2)
626                                 {
627                                         return;
628                                 }
629                                         
630                                 if(data->i >= data->max_overlap)
631                                 {       
632                                         // try to make alloc'ed memory bigger
633                                         data->overlap = realloc(data->overlap, sizeof(BVHTreeOverlap)*data->max_overlap*2);
634                                         
635                                         if(!data->overlap)
636                                         {
637                                                 printf("Out of Memory in traverse\n");
638                                                 return;
639                                         }
640                                         data->max_overlap *= 2;
641                                 }
642                                 
643                                 // both leafs, insert overlap!
644                                 data->overlap[data->i].indexA = node1->index;
645                                 data->overlap[data->i].indexB = node2->index;
646
647                                 data->i++;
648                         }
649                         else
650                         {
651                                 for(j = 0; j < data->tree2->tree_type; j++)
652                                 {
653                                         if(node2->children[j])
654                                                 traverse(data, node1, node2->children[j]);
655                                 }
656                         }
657                 }
658                 else
659                 {
660                         
661                         for(j = 0; j < data->tree2->tree_type; j++)
662                         {
663                                 if(node1->children[j])
664                                         traverse(data, node1->children[j], node2);
665                         }
666                 }
667         }
668         return;
669 }
670
671 BVHTreeOverlap *BLI_bvhtree_overlap(BVHTree *tree1, BVHTree *tree2, int *result)
672 {
673         int j, total = 0;
674         BVHTreeOverlap *overlap = NULL, *to = NULL;
675         BVHOverlapData **data;
676         
677         // check for compatibility of both trees (can't compare 14-DOP with 18-DOP)
678         if((tree1->axis != tree2->axis) && ((tree1->axis == 14) || tree2->axis == 14))
679                 return 0;
680         
681         // fast check root nodes for collision before doing big splitting + traversal
682         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)))
683                 return 0;
684
685         data = MEM_callocN(sizeof(BVHOverlapData *)* tree1->tree_type, "BVHOverlapData_star");
686         
687         for(j = 0; j < tree1->tree_type; j++)
688         {
689                 data[j] = (BVHOverlapData *)MEM_callocN(sizeof(BVHOverlapData), "BVHOverlapData");
690                 
691                 // init BVHOverlapData
692                 data[j]->overlap = (BVHTreeOverlap *)malloc(sizeof(BVHTreeOverlap)*MAX2(tree1->totleaf, tree2->totleaf));
693                 data[j]->tree1 = tree1;
694                 data[j]->tree2 = tree2;
695                 data[j]->max_overlap = MAX2(tree1->totleaf, tree2->totleaf);
696                 data[j]->i = 0;
697         }
698
699 #pragma omp parallel for private(j) schedule(static)
700         for(j = 0; j < MIN2(tree1->tree_type, tree1->nodes[tree1->totleaf]->totnode); j++)
701         {
702                 traverse(data[j], tree1->nodes[tree1->totleaf]->children[j], tree2->nodes[tree2->totleaf]);
703         }
704         
705         for(j = 0; j < tree1->tree_type; j++)
706                 total += data[j]->i;
707         
708         to = overlap = (BVHTreeOverlap *)MEM_callocN(sizeof(BVHTreeOverlap)*total, "BVHTreeOverlap");
709         
710         for(j = 0; j < tree1->tree_type; j++)
711         {
712                 memcpy(to, data[j]->overlap, data[j]->i*sizeof(BVHTreeOverlap));
713                 to+=data[j]->i;
714         }
715         
716         for(j = 0; j < tree1->tree_type; j++)
717         {
718                 free(data[j]->overlap);
719                 MEM_freeN(data[j]);
720         }
721         MEM_freeN(data);
722         
723         (*result) = total;
724         return overlap;
725 }
726
727
728 // bottom up update of bvh tree:
729 // join the 4 children here
730 static void node_join(BVHTree *tree, BVHNode *node)
731 {
732         int i, j;
733         
734         for (i = tree->start_axis; i < tree->stop_axis; i++)
735         {
736                 node->bv[2*i] = FLT_MAX;
737                 node->bv[2*i + 1] = -FLT_MAX;
738         }
739         
740         for (i = 0; i < tree->tree_type; i++)
741         {
742                 if (node->children[i]) 
743                 {
744                         for (j = tree->start_axis; j < tree->stop_axis; j++)
745                         {
746                                 // update minimum 
747                                 if (node->children[i]->bv[(2 * j)] < node->bv[(2 * j)]) 
748                                         node->bv[(2 * j)] = node->children[i]->bv[(2 * j)];
749                                 
750                                 // update maximum 
751                                 if (node->children[i]->bv[(2 * j) + 1] > node->bv[(2 * j) + 1])
752                                         node->bv[(2 * j) + 1] = node->children[i]->bv[(2 * j) + 1];
753                         }
754                 }
755                 else
756                         break;
757         }
758 }
759
760 // call before BLI_bvhtree_update_tree()
761 int BLI_bvhtree_update_node(BVHTree *tree, int index, float *co, float *co_moving, int numpoints)
762 {
763         BVHNode *node= NULL;
764         int i = 0;
765         
766         // check if index exists
767         if(index > tree->totleaf)
768                 return 0;
769         
770         node = tree->nodearray + index;
771         
772         create_kdop_hull(tree, node, co, numpoints, 0);
773         
774         if(co_moving)
775                 create_kdop_hull(tree, node, co_moving, numpoints, 1);
776         
777         // inflate the bv with some epsilon
778         for (i = tree->start_axis; i < tree->stop_axis; i++)
779         {
780                 node->bv[(2 * i)] -= tree->epsilon; // minimum 
781                 node->bv[(2 * i) + 1] += tree->epsilon; // maximum 
782         }
783         
784         return 1;
785 }
786
787 // call BLI_bvhtree_update_node() first for every node/point/triangle
788 void BLI_bvhtree_update_tree(BVHTree *tree)
789 {
790         BVHNode *leaf, *parent;
791         
792         // reset tree traversing flag
793         for (leaf = tree->nodearray + tree->totleaf; leaf != tree->nodearray + tree->totleaf + tree->totbranch; leaf++)
794                 leaf->traversed = 0;
795         
796         for (leaf = tree->nodearray; leaf != tree->nodearray + tree->totleaf; leaf++)
797         {
798                 for (parent = leaf->parent; parent; parent = parent->parent)
799                 {
800                         parent->traversed++;    // we tried to go up in hierarchy 
801                         if (parent->traversed < parent->totnode) 
802                                 break;  // we do not need to check further 
803                         else 
804                                 node_join(tree, parent);
805                 }
806         }
807 }
808
809 float BLI_bvhtree_getepsilon(BVHTree *tree)
810 {
811         return tree->epsilon;
812 }