← Back to branch summary

~ubuntu-branches/ubuntu/trusty/blender/trusty

« back to all changes in this revision

Viewing changes to source/blender/blenlib/intern/BLI_kdopbvh.c

  • Committer: Package Import Robot
  • Author(s): Jeremy Bicha
  • Date: 2013-03-06 12:08:47 UTC
  • mfrom: (1.5.1) (14.1.8 experimental)
  • Revision ID: package-import@ubuntu.com-20130306120847-frjfaryb2zrotwcg
Tags: 2.66a-1ubuntu1
https://launchpad.net/bugs/1076930
https://launchpad.net/bugs/1089256
https://launchpad.net/bugs/1052743
https://launchpad.net/bugs/999024
* Resynchronize with Debian (LP: #1076930, #1089256, #1052743, #999024,
  #1122888, #1147084)
* debian/control:
  - Lower build-depends on libavcodec-dev since we're not
    doing the libav9 transition in Ubuntu yet

Show diffs side-by-side

added added

removed removed

Lines of Context:
source/blender/blenlib/intern/BLI_kdopbvh.c
29
29
 *  \ingroup bli
30
30
 */
31
31
 
32
 
 
33
32
#include <assert.h>
34
33
 
35
34
#include "MEM_guardedalloc.h"
36
35
 
37
36
#include "BLI_utildefines.h"
38
 
 
39
 
 
40
 
 
41
37
#include "BLI_kdopbvh.h"
42
38
#include "BLI_math.h"
43
39
 
45
41
#include <omp.h>
46
42
#endif
47
43
 
48
 
 
49
 
 
50
44
#define MAX_TREETYPE 32
51
 
#define DEFAULT_FIND_NEAREST_HEAP_SIZE 1024
52
 
 
53
 
typedef struct BVHNode
54
 
{
 
45
 
 
46
typedef unsigned char axis_t;
 
47
 
 
48
typedef struct BVHNode {
55
49
        struct BVHNode **children;
56
 
        struct BVHNode *parent; // some user defined traversed need that
 
50
        struct BVHNode *parent; /* some user defined traversed need that */
57
51
        struct BVHNode *skip[2];
58
 
        float *bv;              // Bounding volume of all nodes, max 13 axis
59
 
        int index;              // face, edge, vertex index
60
 
        char totnode;   // how many nodes are used, used for speedup
61
 
        char main_axis; // Axis used to split this node
 
52
        float *bv;      /* Bounding volume of all nodes, max 13 axis */
 
53
        int index;      /* face, edge, vertex index */
 
54
        char totnode;   /* how many nodes are used, used for speedup */
 
55
        char main_axis; /* Axis used to split this node */
62
56
} BVHNode;
63
57
 
64
 
struct BVHTree
65
 
{
 
58
/* keep under 26 bytes for speed purposes */
 
59
struct BVHTree {
66
60
        BVHNode **nodes;
67
 
        BVHNode *nodearray; /* pre-alloc branch nodes */
68
 
        BVHNode **nodechild;    // pre-alloc childs for nodes
69
 
        float   *nodebv;                // pre-alloc bounding-volumes for nodes
70
 
        float   epsilon; /* epslion is used for inflation of the k-dop     */
71
 
        int     totleaf; // leafs
72
 
        int     totbranch;
73
 
        char    tree_type; // type of tree (4 => quadtree)
74
 
        char    axis; // kdop type (6 => OBB, 7 => AABB, ...)
75
 
        char    start_axis, stop_axis; // KDOP_AXES array indices according to axis
 
61
        BVHNode *nodearray;     /* pre-alloc branch nodes */
 
62
        BVHNode **nodechild;    /* pre-alloc childs for nodes */
 
63
        float   *nodebv;        /* pre-alloc bounding-volumes for nodes */
 
64
        float epsilon;          /* epslion is used for inflation of the k-dop      */
 
65
        int totleaf;            /* leafs */
 
66
        int totbranch;
 
67
        axis_t start_axis, stop_axis;  /* KDOP_AXES array indices according to axis */
 
68
        axis_t axis;              /* kdop type (6 => OBB, 7 => AABB, ...) */
 
69
        char tree_type;         /* type of tree (4 => quadtree) */
76
70
};
77
71
 
78
 
typedef struct BVHOverlapData 
79
 
{  
 
72
/* optimization, ensure we stay small */
 
73
BLI_STATIC_ASSERT((sizeof(void *) == 8 && sizeof(BVHTree) <= 48) ||
 
74
                  (sizeof(void *) == 4 && sizeof(BVHTree) <= 32),
 
75
                  "over sized");
 
76
 
 
77
typedef struct BVHOverlapData {
80
78
        BVHTree *tree1, *tree2; 
81
79
        BVHTreeOverlap *overlap; 
82
80
        int i, max_overlap; /* i is number of overlaps */
83
 
        int start_axis, stop_axis;
 
81
        axis_t start_axis, stop_axis;
84
82
} BVHOverlapData;
85
83
 
86
 
typedef struct BVHNearestData
87
 
{
 
84
typedef struct BVHNearestData {
88
85
        BVHTree *tree;
89
 
        const float     *co;
 
86
        const float *co;
90
87
        BVHTree_NearestPointCallback callback;
91
 
        void    *userdata;
92
 
        float proj[13];                 //coordinates projection over axis
 
88
        void    *userdata;
 
89
        float proj[13];         /* coordinates projection over axis */
93
90
        BVHTreeNearest nearest;
94
91
 
95
92
} BVHNearestData;
96
93
 
97
 
typedef struct BVHRayCastData
98
 
{
 
94
typedef struct BVHRayCastData {
99
95
        BVHTree *tree;
100
96
 
101
97
        BVHTree_RayCastCallback callback;
102
 
        void    *userdata;
103
 
 
104
 
 
105
 
        BVHTreeRay    ray;
 
98
        void    *userdata;
 
99
 
 
100
 
 
101
        BVHTreeRay ray;
106
102
        float ray_dot_axis[13];
107
103
        float idot_axis[13];
108
104
        int index[6];
109
105
 
110
106
        BVHTreeRayHit hit;
111
107
} BVHRayCastData;
112
 
////////////////////////////////////////m
113
 
 
114
 
 
115
 
////////////////////////////////////////////////////////////////////////
116
 
// Bounding Volume Hierarchy Definition
117
 
// 
118
 
// Notes: From OBB until 26-DOP --> all bounding volumes possible, just choose type below
119
 
// Notes: You have to choose the type at compile time ITM
120
 
// Notes: You can choose the tree type --> binary, quad, octree, choose below
121
 
////////////////////////////////////////////////////////////////////////
122
 
 
123
 
static float KDOP_AXES[13][3] =
124
 
{ {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},
125
 
{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},
126
 
{0, 1.0, -1.0}
 
108
 
 
109
 
 
110
/*
 
111
 * Bounding Volume Hierarchy Definition
 
112
 *
 
113
 * Notes: From OBB until 26-DOP --> all bounding volumes possible, just choose type below
 
114
 * Notes: You have to choose the type at compile time ITM
 
115
 * Notes: You can choose the tree type --> binary, quad, octree, choose below
 
116
 */
 
117
 
 
118
static float KDOP_AXES[13][3] = {
 
119
        {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},
 
120
        {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},
 
121
        {0, 1.0, -1.0}
127
122
};
128
123
 
 
124
MINLINE axis_t min_axis(axis_t a, axis_t b)
 
125
{
 
126
        return (a < b) ? a : b;
 
127
}
 
128
MINLINE axis_t max_axis(axis_t a, axis_t b)
 
129
{
 
130
        return (b < a) ? a : b;
 
131
}
 
132
 
 
133
#if 0
 
134
 
129
135
/*
130
136
 * Generic push and pop heap
131
137
 */
132
 
#define PUSH_HEAP_BODY(HEAP_TYPE,PRIORITY,heap,heap_size)       \
133
 
{                                                                                                       \
134
 
        HEAP_TYPE element = heap[heap_size-1];                  \
135
 
        int child = heap_size-1;                                                \
136
 
        while (child != 0)                                                              \
137
 
        {                                                                                               \
138
 
                int parent = (child-1) / 2;                                     \
139
 
                if (PRIORITY(element, heap[parent]))                    \
140
 
                {                                                                                       \
141
 
                        heap[child] = heap[parent];                             \
142
 
                        child = parent;                                                 \
143
 
                }                                                                                       \
144
 
                else break;                                                                     \
145
 
        }                                                                                               \
146
 
        heap[child] = element;                                                  \
147
 
}
148
 
 
149
 
#define POP_HEAP_BODY(HEAP_TYPE, PRIORITY,heap,heap_size)       \
150
 
{                                                                                                       \
151
 
        HEAP_TYPE element = heap[heap_size-1];                  \
152
 
        int parent = 0;                                                                 \
153
 
        while (parent < (heap_size-1)/2 )                               \
154
 
        {                                                                                               \
155
 
                int child2 = (parent+1)*2;                                      \
156
 
                if (PRIORITY(heap[child2-1], heap[child2]))     \
157
 
                        --child2;                                                               \
158
 
                                                                                                        \
159
 
                if (PRIORITY(element, heap[child2]))                    \
160
 
                        break;                                                                  \
161
 
                                                                                                        \
162
 
                heap[parent] = heap[child2];                            \
163
 
                parent = child2;                                                        \
164
 
        }                                                                                               \
165
 
        heap[parent] = element;                                                 \
166
 
}
167
 
 
168
 
#if 0
 
138
#define PUSH_HEAP_BODY(HEAP_TYPE, PRIORITY, heap, heap_size)                  \
 
139
        {                                                                         \
 
140
                HEAP_TYPE element = heap[heap_size - 1];                              \
 
141
                int child = heap_size - 1;                                            \
 
142
                while (child != 0)                                                    \
 
143
                {                                                                     \
 
144
                        int parent = (child - 1) / 2;                                     \
 
145
                        if (PRIORITY(element, heap[parent]))                              \
 
146
                        {                                                                 \
 
147
                                heap[child] = heap[parent];                                   \
 
148
                                child = parent;                                               \
 
149
                        }                                                                 \
 
150
                        else break;                                                       \
 
151
                }                                                                     \
 
152
                heap[child] = element;                                                \
 
153
        } (void)0
 
154
 
 
155
#define POP_HEAP_BODY(HEAP_TYPE, PRIORITY, heap, heap_size)                   \
 
156
        {                                                                         \
 
157
                HEAP_TYPE element = heap[heap_size - 1];                              \
 
158
                int parent = 0;                                                       \
 
159
                while (parent < (heap_size - 1) / 2)                                  \
 
160
                {                                                                     \
 
161
                        int child2 = (parent + 1) * 2;                                    \
 
162
                        if (PRIORITY(heap[child2 - 1], heap[child2])) {                   \
 
163
                                child2--;                                                     \
 
164
                        }                                                                 \
 
165
                        if (PRIORITY(element, heap[child2])) {                            \
 
166
                                break;                                                        \
 
167
                        }                                                                 \
 
168
                        heap[parent] = heap[child2];                                      \
 
169
                        parent = child2;                                                  \
 
170
                }                                                                     \
 
171
                heap[parent] = element;                                               \
 
172
        } (void)0
 
173
 
169
174
static int ADJUST_MEMORY(void *local_memblock, void **memblock, int new_size, int *max_size, int size_per_item)
170
175
{
171
 
        int   new_max_size = *max_size * 2;
 
176
        int new_max_size = *max_size * 2;
172
177
        void *new_memblock = NULL;
173
178
 
174
179
        if (new_size <= *max_size)
176
181
 
177
182
        if (*memblock == local_memblock)
178
183
        {
179
 
                new_memblock = malloc( size_per_item * new_max_size );
180
 
                memcpy( new_memblock, *memblock, size_per_item * *max_size );
 
184
                new_memblock = malloc(size_per_item * new_max_size);
 
185
                memcpy(new_memblock, *memblock, size_per_item * *max_size);
181
186
        }
182
187
        else
183
 
                new_memblock = realloc(*memblock, size_per_item * new_max_size );
 
188
                new_memblock = realloc(*memblock, size_per_item * new_max_size);
184
189
 
185
190
        if (new_memblock)
186
191
        {
193
198
}
194
199
#endif
195
200
 
196
 
//////////////////////////////////////////////////////////////////////////////////////////////////////
197
 
// Introsort 
198
 
// with permission deriven from the following Java code:
199
 
// http://ralphunden.net/content/tutorials/a-guide-to-introsort/
200
 
// and he derived it from the SUN STL 
201
 
//////////////////////////////////////////////////////////////////////////////////////////////////////
 
201
/*
 
202
 * Introsort
 
203
 * with permission deriven from the following Java code:
 
204
 * http://ralphunden.net/content/tutorials/a-guide-to-introsort/
 
205
 * and he derived it from the SUN STL
 
206
 */
 
207
 
202
208
//static int size_threshold = 16;
 
209
 
203
210
/*
204
211
 * Common methods for all algorithms
205
212
 */
206
213
#if 0
207
214
static int floor_lg(int a)
208
215
{
209
 
        return (int)(floor(log(a)/log(2)));
 
216
        return (int)(floor(log(a) / log(2)));
210
217
}
211
218
#endif
212
219
 
215
222
 */
216
223
static void bvh_insertionsort(BVHNode **a, int lo, int hi, int axis)
217
224
{
218
 
        int i,j;
 
225
        int i, j;
219
226
        BVHNode *t;
220
 
        for (i=lo; i < hi; i++)
221
 
        {
222
 
                j=i;
 
227
        for (i = lo; i < hi; i++) {
 
228
                j = i;
223
229
                t = a[i];
224
 
                while ((j!=lo) && (t->bv[axis] < (a[j-1])->bv[axis]))
225
 
                {
226
 
                        a[j] = a[j-1];
 
230
                while ((j != lo) && (t->bv[axis] < (a[j - 1])->bv[axis])) {
 
231
                        a[j] = a[j - 1];
227
232
                        j--;
228
233
                }
229
234
                a[j] = t;
230
235
        }
231
236
}
232
237
 
233
 
static int bvh_partition(BVHNode **a, int lo, int hi, BVHNode * x, int axis)
 
238
static int bvh_partition(BVHNode **a, int lo, int hi, BVHNode *x, int axis)
234
239
{
235
 
        int i=lo, j=hi;
236
 
        while (1)
237
 
        {
 
240
        int i = lo, j = hi;
 
241
        while (1) {
238
242
                while ((a[i])->bv[axis] < x->bv[axis]) i++;
239
243
                j--;
240
244
                while (x->bv[axis] < (a[j])->bv[axis]) j--;
241
245
                if (!(i < j))
242
246
                        return i;
243
 
                SWAP( BVHNode* , a[i], a[j]);
 
247
                SWAP(BVHNode *, a[i], a[j]);
244
248
                i++;
245
249
        }
246
250
}
251
255
#if 0
252
256
static void bvh_downheap(BVHNode **a, int i, int n, int lo, int axis)
253
257
{
254
 
        BVHNode * d = a[lo+i-1];
 
258
        BVHNode *d = a[lo + i - 1];
255
259
        int child;
256
 
        while (i<=n/2)
 
260
        while (i <= n / 2)
257
261
        {
258
 
                child = 2*i;
259
 
                if ((child < n) && ((a[lo+child-1])->bv[axis] < (a[lo+child])->bv[axis]))
 
262
                child = 2 * i;
 
263
                if ((child < n) && ((a[lo + child - 1])->bv[axis] < (a[lo + child])->bv[axis]))
260
264
                {
261
265
                        child++;
262
266
                }
263
 
                if (!(d->bv[axis] < (a[lo+child-1])->bv[axis])) break;
264
 
                a[lo+i-1] = a[lo+child-1];
 
267
                if (!(d->bv[axis] < (a[lo + child - 1])->bv[axis])) break;
 
268
                a[lo + i - 1] = a[lo + child - 1];
265
269
                i = child;
266
270
        }
267
 
        a[lo+i-1] = d;
 
271
        a[lo + i - 1] = d;
268
272
}
269
273
 
270
274
static void bvh_heapsort(BVHNode **a, int lo, int hi, int axis)
271
275
{
272
 
        int n = hi-lo, i;
273
 
        for (i=n/2; i>=1; i=i-1)
 
276
        int n = hi - lo, i;
 
277
        for (i = n / 2; i >= 1; i = i - 1)
274
278
        {
275
 
                bvh_downheap(a, i,n,lo, axis);
 
279
                bvh_downheap(a, i, n, lo, axis);
276
280
        }
277
 
        for (i=n; i>1; i=i-1)
 
281
        for (i = n; i > 1; i = i - 1)
278
282
        {
279
 
                SWAP(BVHNode*, a[lo],a[lo+i-1]);
280
 
                bvh_downheap(a, 1,i-1,lo, axis);
 
283
                SWAP(BVHNode *, a[lo], a[lo + i - 1]);
 
284
                bvh_downheap(a, 1, i - 1, lo, axis);
281
285
        }
282
286
}
283
287
#endif
284
288
 
285
 
static BVHNode *bvh_medianof3(BVHNode **a, int lo, int mid, int hi, int axis) // returns Sortable
 
289
static BVHNode *bvh_medianof3(BVHNode **a, int lo, int mid, int hi, int axis)  /* returns Sortable */
286
290
{
287
 
        if ((a[mid])->bv[axis] < (a[lo])->bv[axis])
288
 
        {
 
291
        if ((a[mid])->bv[axis] < (a[lo])->bv[axis]) {
289
292
                if ((a[hi])->bv[axis] < (a[mid])->bv[axis])
290
293
                        return a[mid];
291
 
                else
292
 
                {
 
294
                else {
293
295
                        if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
294
296
                                return a[hi];
295
297
                        else
296
298
                                return a[lo];
297
299
                }
298
300
        }
299
 
        else
300
 
        {
301
 
                if ((a[hi])->bv[axis] < (a[mid])->bv[axis])
302
 
                {
 
301
        else {
 
302
                if ((a[hi])->bv[axis] < (a[mid])->bv[axis]) {
303
303
                        if ((a[hi])->bv[axis] < (a[lo])->bv[axis])
304
304
                                return a[lo];
305
305
                        else
314
314
/*
315
315
 * Quicksort algorithm modified for Introsort
316
316
 */
317
 
static void bvh_introsort_loop (BVHNode **a, int lo, int hi, int depth_limit, int axis)
 
317
static void bvh_introsort_loop(BVHNode **a, int lo, int hi, int depth_limit, int axis)
318
318
{
319
319
        int p;
320
320
 
321
 
        while (hi-lo > size_threshold)
 
321
        while (hi - lo > size_threshold)
322
322
        {
323
323
                if (depth_limit == 0)
324
324
                {
325
325
                        bvh_heapsort(a, lo, hi, axis);
326
326
                        return;
327
327
                }
328
 
                depth_limit=depth_limit-1;
329
 
                p=bvh_partition(a, lo, hi, bvh_medianof3(a, lo, lo+((hi-lo)/2)+1, hi-1, axis), axis);
 
328
                depth_limit = depth_limit - 1;
 
329
                p = bvh_partition(a, lo, hi, bvh_medianof3(a, lo, lo + ((hi - lo) / 2) + 1, hi - 1, axis), axis);
330
330
                bvh_introsort_loop(a, p, hi, depth_limit, axis);
331
 
                hi=p;
 
331
                hi = p;
332
332
        }
333
333
}
334
334
 
336
336
{
337
337
        if (begin < end)
338
338
        {
339
 
                BVHNode **a=a0;
340
 
                bvh_introsort_loop(a, begin, end, 2*floor_lg(end-begin), axis);
 
339
                BVHNode **a = a0;
 
340
                bvh_introsort_loop(a, begin, end, 2 * floor_lg(end - begin), axis);
341
341
                bvh_insertionsort(a, begin, end, axis);
342
342
        }
343
343
}
348
348
}
349
349
#endif
350
350
 
351
 
//after a call to this function you can expect one of:
352
 
//      every node to left of a[n] are smaller or equal to it
353
 
//      every node to the right of a[n] are greater or equal to it
 
351
/* after a call to this function you can expect one of:
 
352
 * - every node to left of a[n] are smaller or equal to it
 
353
 * - every node to the right of a[n] are greater or equal to it */
354
354
static int partition_nth_element(BVHNode **a, int _begin, int _end, int n, int axis)
355
355
{
356
356
        int begin = _begin, end = _end, cut;
357
 
        while (end-begin > 3)
358
 
        {
359
 
                cut = bvh_partition(a, begin, end, bvh_medianof3(a, begin, (begin+end)/2, end-1, axis), axis );
 
357
        while (end - begin > 3) {
 
358
                cut = bvh_partition(a, begin, end, bvh_medianof3(a, begin, (begin + end) / 2, end - 1, axis), axis);
360
359
                if (cut <= n)
361
360
                        begin = cut;
362
361
                else
367
366
        return n;
368
367
}
369
368
 
370
 
//////////////////////////////////////////////////////////////////////////////////////////////////////
 
369
/* --- */
371
370
static void build_skip_links(BVHTree *tree, BVHNode *node, BVHNode *left, BVHNode *right)
372
371
{
373
372
        int i;
375
374
        node->skip[0] = left;
376
375
        node->skip[1] = right;
377
376
        
378
 
        for (i = 0; i < node->totnode; i++)
379
 
        {
380
 
                if (i+1 < node->totnode)
381
 
                        build_skip_links(tree, node->children[i], left, node->children[i+1] );
 
377
        for (i = 0; i < node->totnode; i++) {
 
378
                if (i + 1 < node->totnode)
 
379
                        build_skip_links(tree, node->children[i], left, node->children[i + 1]);
382
380
                else
383
 
                        build_skip_links(tree, node->children[i], left, right );
 
381
                        build_skip_links(tree, node->children[i], left, right);
384
382
 
385
383
                left = node->children[i];
386
384
        }
393
391
{
394
392
        float newminmax;
395
393
        float *bv = node->bv;
396
 
        int i, k;
 
394
        int k;
 
395
        axis_t axis_iter;
397
396
        
398
 
        // don't init boudings for the moving case
399
 
        if (!moving)
400
 
        {
401
 
                for (i = tree->start_axis; i < tree->stop_axis; i++)
402
 
                {
403
 
                        bv[2*i] = FLT_MAX;
404
 
                        bv[2*i + 1] = -FLT_MAX;
 
397
        /* don't init boudings for the moving case */
 
398
                if (!moving) {
 
399
                for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
 
400
                        bv[2 * axis_iter] = FLT_MAX;
 
401
                        bv[2 * axis_iter + 1] = -FLT_MAX;
405
402
                }
406
403
        }
407
404
        
408
 
        for (k = 0; k < numpoints; k++)
409
 
        {
410
 
                // for all Axes.
411
 
                for (i = tree->start_axis; i < tree->stop_axis; i++)
412
 
                {
413
 
                        newminmax = dot_v3v3(&co[k * 3], KDOP_AXES[i]);
414
 
                        if (newminmax < bv[2 * i])
415
 
                                bv[2 * i] = newminmax;
416
 
                        if (newminmax > bv[(2 * i) + 1])
417
 
                                bv[(2 * i) + 1] = newminmax;
 
405
        for (k = 0; k < numpoints; k++) {
 
406
                /* for all Axes. */
 
407
                for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
 
408
                        newminmax = dot_v3v3(&co[k * 3], KDOP_AXES[axis_iter]);
 
409
                        if (newminmax < bv[2 * axis_iter])
 
410
                                bv[2 * axis_iter] = newminmax;
 
411
                        if (newminmax > bv[(2 * axis_iter) + 1])
 
412
                                bv[(2 * axis_iter) + 1] = newminmax;
418
413
                }
419
414
        }
420
415
}
421
416
 
422
 
// depends on the fact that the BVH's for each face is already build
 
417
/* depends on the fact that the BVH's for each face is already build */
423
418
static void refit_kdop_hull(BVHTree *tree, BVHNode *node, int start, int end)
424
419
{
425
 
        float newmin,newmax;
426
 
        int i, j;
 
420
        float newmin, newmax;
427
421
        float *bv = node->bv;
 
422
        int j;
 
423
        axis_t axis_iter;
428
424
 
429
 
        
430
 
        for (i = tree->start_axis; i < tree->stop_axis; i++)
431
 
        {
432
 
                bv[2*i] = FLT_MAX;
433
 
                bv[2*i + 1] = -FLT_MAX;
 
425
        for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
 
426
                bv[(2 * axis_iter)] = FLT_MAX;
 
427
                bv[(2 * axis_iter) + 1] = -FLT_MAX;
434
428
        }
435
429
 
436
 
        for (j = start; j < end; j++)
437
 
        {
438
 
// for all Axes.
439
 
                for (i = tree->start_axis; i < tree->stop_axis; i++)
440
 
                {
441
 
                        newmin = tree->nodes[j]->bv[(2 * i)];   
442
 
                        if ((newmin < bv[(2 * i)]))
443
 
                                bv[(2 * i)] = newmin;
444
 
 
445
 
                        newmax = tree->nodes[j]->bv[(2 * i) + 1];
446
 
                        if ((newmax > bv[(2 * i) + 1]))
447
 
                                bv[(2 * i) + 1] = newmax;
 
430
        for (j = start; j < end; j++) {
 
431
                /* for all Axes. */
 
432
                for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
 
433
                        newmin = tree->nodes[j]->bv[(2 * axis_iter)];
 
434
                        if ((newmin < bv[(2 * axis_iter)]))
 
435
                                bv[(2 * axis_iter)] = newmin;
 
436
 
 
437
                        newmax = tree->nodes[j]->bv[(2 * axis_iter) + 1];
 
438
                        if ((newmax > bv[(2 * axis_iter) + 1]))
 
439
                                bv[(2 * axis_iter) + 1] = newmax;
448
440
                }
449
441
        }
450
442
 
451
443
}
452
444
 
453
 
// only supports x,y,z axis in the moment
454
 
// but we should use a plain and simple function here for speed sake
 
445
/* only supports x,y,z axis in the moment
 
446
 * but we should use a plain and simple function here for speed sake */
455
447
static char get_largest_axis(float *bv)
456
448
{
457
449
        float middle_point[3];
458
450
 
459
 
        middle_point[0] = (bv[1]) - (bv[0]); // x axis
460
 
        middle_point[1] = (bv[3]) - (bv[2]); // y axis
461
 
        middle_point[2] = (bv[5]) - (bv[4]); // z axis
462
 
        if (middle_point[0] > middle_point[1]) 
463
 
        {
 
451
        middle_point[0] = (bv[1]) - (bv[0]); /* x axis */
 
452
        middle_point[1] = (bv[3]) - (bv[2]); /* y axis */
 
453
        middle_point[2] = (bv[5]) - (bv[4]); /* z axis */
 
454
        if (middle_point[0] > middle_point[1]) {
464
455
                if (middle_point[0] > middle_point[2])
465
 
                        return 1; // max x axis
 
456
                        return 1;  /* max x axis */
466
457
                else
467
 
                        return 5; // max z axis
 
458
                        return 5;  /* max z axis */
468
459
        }
469
 
        else 
470
 
        {
 
460
        else {
471
461
                if (middle_point[1] > middle_point[2])
472
 
                        return 3; // max y axis
 
462
                        return 3;  /* max y axis */
473
463
                else
474
 
                        return 5; // max z axis
 
464
                        return 5;  /* max z axis */
475
465
        }
476
466
}
477
467
 
478
 
// bottom-up update of bvh node BV
479
 
// join the children on the parent BV
 
468
/* bottom-up update of bvh node BV
 
469
 * join the children on the parent BV */
480
470
static void node_join(BVHTree *tree, BVHNode *node)
481
471
{
482
 
        int i, j;
483
 
        
484
 
        for (i = tree->start_axis; i < tree->stop_axis; i++)
485
 
        {
486
 
                node->bv[2*i] = FLT_MAX;
487
 
                node->bv[2*i + 1] = -FLT_MAX;
 
472
        int i;
 
473
        axis_t axis_iter;
 
474
 
 
475
        for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
 
476
                node->bv[(2 * axis_iter)] = FLT_MAX;
 
477
                node->bv[(2 * axis_iter) + 1] = -FLT_MAX;
488
478
        }
489
479
        
490
 
        for (i = 0; i < tree->tree_type; i++)
491
 
        {
492
 
                if (node->children[i]) 
493
 
                {
494
 
                        for (j = tree->start_axis; j < tree->stop_axis; j++)
495
 
                        {
496
 
                                // update minimum 
497
 
                                if (node->children[i]->bv[(2 * j)] < node->bv[(2 * j)]) 
498
 
                                        node->bv[(2 * j)] = node->children[i]->bv[(2 * j)];
499
 
                                
500
 
                                // update maximum 
501
 
                                if (node->children[i]->bv[(2 * j) + 1] > node->bv[(2 * j) + 1])
502
 
                                        node->bv[(2 * j) + 1] = node->children[i]->bv[(2 * j) + 1];
 
480
        for (i = 0; i < tree->tree_type; i++) {
 
481
                if (node->children[i]) {
 
482
                        for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
 
483
                                /* update minimum */
 
484
                                if (node->children[i]->bv[(2 * axis_iter)] < node->bv[(2 * axis_iter)])
 
485
                                        node->bv[(2 * axis_iter)] = node->children[i]->bv[(2 * axis_iter)];
 
486
 
 
487
                                /* update maximum */
 
488
                                if (node->children[i]->bv[(2 * axis_iter) + 1] > node->bv[(2 * axis_iter) + 1])
 
489
                                        node->bv[(2 * axis_iter) + 1] = node->children[i]->bv[(2 * axis_iter) + 1];
503
490
                        }
504
491
                }
505
492
                else
514
501
static void bvhtree_print_tree(BVHTree *tree, BVHNode *node, int depth)
515
502
{
516
503
        int i;
517
 
        for (i=0; i<depth; i++) printf(" ");
 
504
        axis_t axis_iter;
 
505
 
 
506
        for (i = 0; i < depth; i++) printf(" ");
518
507
        printf(" - %d (%ld): ", node->index, node - tree->nodearray);
519
 
        for (i=2*tree->start_axis; i<2*tree->stop_axis; i++)
520
 
                printf("%.3f ", node->bv[i]);
 
508
        for (axis_iter = 2 * tree->start_axis; axis_iter < 2 * tree->stop_axis; axis_iter++)
 
509
                printf("%.3f ", node->bv[axis_iter]);
521
510
        printf("\n");
522
511
 
523
 
        for (i=0; i<tree->tree_type; i++)
 
512
        for (i = 0; i < tree->tree_type; i++)
524
513
                if (node->children[i])
525
 
                        bvhtree_print_tree(tree, node->children[i], depth+1);
 
514
                        bvhtree_print_tree(tree, node->children[i], depth + 1);
526
515
}
527
516
 
528
517
static void bvhtree_info(BVHTree *tree)
530
519
        printf("BVHTree info\n");
531
520
        printf("tree_type = %d, axis = %d, epsilon = %f\n", tree->tree_type, tree->axis, tree->epsilon);
532
521
        printf("nodes = %d, branches = %d, leafs = %d\n", tree->totbranch + tree->totleaf,  tree->totbranch, tree->totleaf);
533
 
        printf("Memory per node = %ldbytes\n", sizeof(BVHNode) + sizeof(BVHNode*)*tree->tree_type + sizeof(float)*tree->axis);
 
522
        printf("Memory per node = %ldbytes\n", sizeof(BVHNode) + sizeof(BVHNode *) * tree->tree_type + sizeof(float) * tree->axis);
534
523
        printf("BV memory = %dbytes\n", MEM_allocN_len(tree->nodebv));
535
524
 
536
 
        printf("Total memory = %ldbytes\n", sizeof(BVHTree)
537
 
                + MEM_allocN_len(tree->nodes)
538
 
                + MEM_allocN_len(tree->nodearray)
539
 
                + MEM_allocN_len(tree->nodechild)
540
 
                + MEM_allocN_len(tree->nodebv)
541
 
                );
 
525
        printf("Total memory = %ldbytes\n", sizeof(BVHTree) +
 
526
               MEM_allocN_len(tree->nodes) +
 
527
               MEM_allocN_len(tree->nodearray) +
 
528
               MEM_allocN_len(tree->nodechild) +
 
529
               MEM_allocN_len(tree->nodebv));
542
530
 
543
531
//      bvhtree_print_tree(tree, tree->nodes[tree->totleaf], 0);
544
532
}
551
539
{
552
540
        int i, j, check = 0;
553
541
        
554
 
        // check the pointer list
 
542
        /* check the pointer list */
555
543
        for (i = 0; i < tree->totleaf; i++)
556
544
        {
557
 
                if (tree->nodes[i]->parent == NULL)
 
545
                if (tree->nodes[i]->parent == NULL) {
558
546
                        printf("Leaf has no parent: %d\n", i);
559
 
                else
560
 
                {
 
547
                }
 
548
                else {
561
549
                        for (j = 0; j < tree->tree_type; j++)
562
550
                        {
563
551
                                if (tree->nodes[i]->parent->children[j] == tree->nodes[i])
571
559
                }
572
560
        }
573
561
        
574
 
        // check the leaf list
 
562
        /* check the leaf list */
575
563
        for (i = 0; i < tree->totleaf; i++)
576
564
        {
577
 
                if (tree->nodearray[i].parent == NULL)
 
565
                if (tree->nodearray[i].parent == NULL) {
578
566
                        printf("Leaf has no parent: %d\n", i);
579
 
                else
580
 
                {
 
567
                }
 
568
                else {
581
569
                        for (j = 0; j < tree->tree_type; j++)
582
570
                        {
583
571
                                if (tree->nodearray[i].parent->children[j] == &tree->nodearray[i])
595
583
}
596
584
#endif
597
585
 
598
 
//Helper data and structures to build a min-leaf generalized implicit tree
599
 
//This code can be easily reduced (basicly this is only method to calculate pow(k, n) in O(1).. and stuff like that)
600
 
typedef struct BVHBuildHelper
601
 
{
602
 
        int tree_type;                          /* */
603
 
        int totleafs;                           /* */
604
 
 
605
 
        int leafs_per_child[32];        /* Min number of leafs that are archievable from a node at depth N */
606
 
        int branches_on_level[32];      /* Number of nodes at depth N (tree_type^N) */
607
 
 
608
 
        int remain_leafs;                       /* Number of leafs that are placed on the level that is not 100% filled */
 
586
/* Helper data and structures to build a min-leaf generalized implicit tree
 
587
 * This code can be easily reduced (basicly this is only method to calculate pow(k, n) in O(1).. and stuff like that) */
 
588
typedef struct BVHBuildHelper {
 
589
        int tree_type;              /* */
 
590
        int totleafs;               /* */
 
591
 
 
592
        int leafs_per_child[32];    /* Min number of leafs that are archievable from a node at depth N */
 
593
        int branches_on_level[32];  /* Number of nodes at depth N (tree_type^N) */
 
594
 
 
595
        int remain_leafs;           /* Number of leafs that are placed on the level that is not 100% filled */
609
596
 
610
597
} BVHBuildHelper;
611
598
 
616
603
        int nnodes;
617
604
 
618
605
        data->totleafs = tree->totleaf;
619
 
        data->tree_type= tree->tree_type;
 
606
        data->tree_type = tree->tree_type;
620
607
 
621
 
        //Calculate the smallest tree_type^n such that tree_type^n >= num_leafs
622
 
        for (
623
 
                data->leafs_per_child[0] = 1;
624
 
                data->leafs_per_child[0] <  data->totleafs;
625
 
                data->leafs_per_child[0] *= data->tree_type
626
 
        );
 
608
        /* Calculate the smallest tree_type^n such that tree_type^n >= num_leafs */
 
609
        for (data->leafs_per_child[0] = 1;
 
610
             data->leafs_per_child[0] <  data->totleafs;
 
611
             data->leafs_per_child[0] *= data->tree_type)
 
612
        {
 
613
                /* pass */
 
614
        }
627
615
 
628
616
        data->branches_on_level[0] = 1;
629
617
 
630
 
        //We could stop the loop first (but I am lazy to find out when)
631
 
        for (depth = 1; depth < 32; depth++)
632
 
        {
633
 
                data->branches_on_level[depth] = data->branches_on_level[depth-1] * data->tree_type;
634
 
                data->leafs_per_child  [depth] = data->leafs_per_child  [depth-1] / data->tree_type;
 
618
        /* We could stop the loop first (but I am lazy to find out when) */
 
619
        for (depth = 1; depth < 32; depth++) {
 
620
                data->branches_on_level[depth] = data->branches_on_level[depth - 1] * data->tree_type;
 
621
                data->leafs_per_child[depth] = data->leafs_per_child[depth - 1] / data->tree_type;
635
622
        }
636
623
 
637
624
        remain = data->totleafs - data->leafs_per_child[1];
642
629
// return the min index of all the leafs archivable with the given branch
643
630
static int implicit_leafs_index(BVHBuildHelper *data, int depth, int child_index)
644
631
{
645
 
        int min_leaf_index = child_index * data->leafs_per_child[depth-1];
 
632
        int min_leaf_index = child_index * data->leafs_per_child[depth - 1];
646
633
        if (min_leaf_index <= data->remain_leafs)
647
634
                return min_leaf_index;
648
635
        else if (data->leafs_per_child[depth])
649
 
                return data->totleafs - (data->branches_on_level[depth-1] - child_index) * data->leafs_per_child[depth];
 
636
                return data->totleafs - (data->branches_on_level[depth - 1] - child_index) * data->leafs_per_child[depth];
650
637
        else
651
638
                return data->remain_leafs;
652
639
}
679
666
 *    (looping elements, knowing if its a leaf or not.. etc...)
680
667
 */
681
668
 
682
 
// This functions returns the number of branches needed to have the requested number of leafs.
 
669
/* This functions returns the number of branches needed to have the requested number of leafs. */
683
670
static int implicit_needed_branches(int tree_type, int leafs)
684
671
{
685
 
        return MAX2(1, (leafs + tree_type - 3) / (tree_type-1) );
 
672
        return max_ii(1, (leafs + tree_type - 3) / (tree_type - 1) );
686
673
}
687
674
 
688
675
/*
693
680
 *  - if all elements are different all partition will get the same subset of elements
694
681
 *    as if the array was sorted.
695
682
 *
696
 
 * partition P is described as the elements in the range ( nth[P] , nth[P+1] ]
 
683
 * partition P is described as the elements in the range ( nth[P], nth[P+1] ]
697
684
 *
698
685
 * TODO: This can be optimized a bit by doing a specialized nth_element instead of K nth_elements
699
686
 */
700
687
static void split_leafs(BVHNode **leafs_array, int *nth, int partitions, int split_axis)
701
688
{
702
689
        int i;
703
 
        for (i=0; i < partitions-1; i++)
704
 
        {
 
690
        for (i = 0; i < partitions - 1; i++) {
705
691
                if (nth[i] >= nth[partitions])
706
692
                        break;
707
693
 
708
 
                partition_nth_element(leafs_array, nth[i], nth[partitions], nth[i+1], split_axis);
 
694
                partition_nth_element(leafs_array, nth[i], nth[partitions], nth[i + 1], split_axis);
709
695
        }
710
696
}
711
697
 
722
708
 * The reason is that we can build level N+1 from level N without any data dependencies.. thus it allows
723
709
 * to use multithread building.
724
710
 *
725
 
 * To archieve this is necessary to find how much leafs are accessible from a certain branch, BVHBuildHelper
 
711
 * To archive this is necessary to find how much leafs are accessible from a certain branch, BVHBuildHelper
726
712
 * implicit_needed_branches and implicit_leafs_index are auxiliary functions to solve that "optimal-split".
727
713
 */
728
714
static void non_recursive_bvh_div_nodes(BVHTree *tree, BVHNode *branches_array, BVHNode **leafs_array, int num_leafs)
730
716
        int i;
731
717
 
732
718
        const int tree_type   = tree->tree_type;
733
 
        const int tree_offset = 2 - tree->tree_type; //this value is 0 (on binary trees) and negative on the others
734
 
        const int num_branches= implicit_needed_branches(tree_type, num_leafs);
 
719
        const int tree_offset = 2 - tree->tree_type; /* this value is 0 (on binary trees) and negative on the others */
 
720
        const int num_branches = implicit_needed_branches(tree_type, num_leafs);
735
721
 
736
722
        BVHBuildHelper data;
737
723
        int depth;
738
724
        
739
 
        // set parent from root node to NULL
740
 
        BVHNode *tmp = branches_array+0;
 
725
        /* set parent from root node to NULL */
 
726
        BVHNode *tmp = branches_array + 0;
741
727
        tmp->parent = NULL;
742
728
 
743
 
        //Most of bvhtree code relies on 1-leaf trees having at least one branch
744
 
        //We handle that special case here
745
 
        if (num_leafs == 1)
746
 
        {
747
 
                BVHNode *root = branches_array+0;
 
729
        /* Most of bvhtree code relies on 1-leaf trees having at least one branch
 
730
         * We handle that special case here */
 
731
        if (num_leafs == 1) {
 
732
                BVHNode *root = branches_array + 0;
748
733
                refit_kdop_hull(tree, root, 0, num_leafs);
749
734
                root->main_axis = get_largest_axis(root->bv) / 2;
750
735
                root->totnode = 1;
753
738
                return;
754
739
        }
755
740
 
756
 
        branches_array--;       //Implicit trees use 1-based indexs
757
 
        
 
741
        branches_array--;  /* Implicit trees use 1-based indexs */
 
742
 
758
743
        build_implicit_tree_helper(tree, &data);
759
744
 
760
 
        //Loop tree levels (log N) loops
761
 
        for (i=1, depth = 1; i <= num_branches; i = i*tree_type + tree_offset, depth++)
762
 
        {
763
 
                const int first_of_next_level = i*tree_type + tree_offset;
764
 
                const int  end_j = MIN2(first_of_next_level, num_branches + 1); //index of last branch on this level
 
745
        /* Loop tree levels (log N) loops */
 
746
        for (i = 1, depth = 1; i <= num_branches; i = i * tree_type + tree_offset, depth++) {
 
747
                const int first_of_next_level = i * tree_type + tree_offset;
 
748
                const int end_j = min_ii(first_of_next_level, num_branches + 1);  /* index of last branch on this level */
765
749
                int j;
766
750
 
767
 
                //Loop all branches on this level
 
751
                /* Loop all branches on this level */
768
752
#pragma omp parallel for private(j) schedule(static)
769
753
                for (j = i; j < end_j; j++) {
770
754
                        int k;
771
 
                        const int parent_level_index= j-i;
772
 
                        BVHNode* parent = branches_array + j;
773
 
                        int nth_positions[ MAX_TREETYPE + 1];
 
755
                        const int parent_level_index = j - i;
 
756
                        BVHNode *parent = branches_array + j;
 
757
                        int nth_positions[MAX_TREETYPE + 1];
774
758
                        char split_axis;
775
759
 
776
760
                        int parent_leafs_begin = implicit_leafs_index(&data, depth, parent_level_index);
777
 
                        int parent_leafs_end   = implicit_leafs_index(&data, depth, parent_level_index+1);
 
761
                        int parent_leafs_end   = implicit_leafs_index(&data, depth, parent_level_index + 1);
778
762
 
779
 
                        //This calculates the bounding box of this branch
780
 
                        //and chooses the largest axis as the axis to divide leafs
 
763
                        /* This calculates the bounding box of this branch
 
764
                         * and chooses the largest axis as the axis to divide leafs */
781
765
                        refit_kdop_hull(tree, parent, parent_leafs_begin, parent_leafs_end);
782
766
                        split_axis = get_largest_axis(parent->bv);
783
767
 
784
 
                        //Save split axis (this can be used on raytracing to speedup the query time)
 
768
                        /* Save split axis (this can be used on raytracing to speedup the query time) */
785
769
                        parent->main_axis = split_axis / 2;
786
770
 
787
 
                        //Split the childs along the split_axis, note: its not needed to sort the whole leafs array
788
 
                        //Only to assure that the elements are partioned on a way that each child takes the elements
789
 
                        //it would take in case the whole array was sorted.
790
 
                        //Split_leafs takes care of that "sort" problem.
791
 
                        nth_positions[        0] = parent_leafs_begin;
 
771
                        /* Split the childs along the split_axis, note: its not needed to sort the whole leafs array
 
772
                         * Only to assure that the elements are partitioned on a way that each child takes the elements
 
773
                         * it would take in case the whole array was sorted.
 
774
                         * Split_leafs takes care of that "sort" problem. */
 
775
                        nth_positions[0] = parent_leafs_begin;
792
776
                        nth_positions[tree_type] = parent_leafs_end;
793
777
                        for (k = 1; k < tree_type; k++) {
794
778
                                int child_index = j * tree_type + tree_offset + k;
795
 
                                int child_level_index = child_index - first_of_next_level; //child level index
796
 
                                nth_positions[k] = implicit_leafs_index(&data, depth+1, child_level_index);
 
779
                                int child_level_index = child_index - first_of_next_level; /* child level index */
 
780
                                nth_positions[k] = implicit_leafs_index(&data, depth + 1, child_level_index);
797
781
                        }
798
782
 
799
783
                        split_leafs(leafs_array, nth_positions, tree_type, split_axis);
800
784
 
801
785
 
802
 
                        //Setup children and totnode counters
803
 
                        //Not really needed but currently most of BVH code relies on having an explicit children structure
 
786
                        /* Setup children and totnode counters
 
787
                         * Not really needed but currently most of BVH code relies on having an explicit children structure */
804
788
                        for (k = 0; k < tree_type; k++) {
805
789
                                int child_index = j * tree_type + tree_offset + k;
806
 
                                int child_level_index = child_index - first_of_next_level; //child level index
 
790
                                int child_level_index = child_index - first_of_next_level; /* child level index */
807
791
 
808
 
                                int child_leafs_begin = implicit_leafs_index(&data, depth+1, child_level_index);
809
 
                                int child_leafs_end   = implicit_leafs_index(&data, depth+1, child_level_index+1);
 
792
                                int child_leafs_begin = implicit_leafs_index(&data, depth + 1, child_level_index);
 
793
                                int child_leafs_end   = implicit_leafs_index(&data, depth + 1, child_level_index + 1);
810
794
 
811
795
                                if (child_leafs_end - child_leafs_begin > 1) {
812
796
                                        parent->children[k] = branches_array + child_index;
813
797
                                        parent->children[k]->parent = parent;
814
798
                                }
815
799
                                else if (child_leafs_end - child_leafs_begin == 1) {
816
 
                                        parent->children[k] = leafs_array[ child_leafs_begin ];
 
800
                                        parent->children[k] = leafs_array[child_leafs_begin];
817
801
                                        parent->children[k]->parent = parent;
818
802
                                }
819
803
                                else {
820
804
                                        break;
821
805
                                }
822
806
 
823
 
                                parent->totnode = k+1;
 
807
                                parent->totnode = k + 1;
824
808
                        }
825
809
                }
826
810
        }
835
819
        BVHTree *tree;
836
820
        int numnodes, i;
837
821
        
838
 
        // theres not support for trees below binary-trees :P
 
822
        /* theres not support for trees below binary-trees :P */
839
823
        if (tree_type < 2)
840
824
                return NULL;
841
 
        
 
825
 
842
826
        if (tree_type > MAX_TREETYPE)
843
827
                return NULL;
844
828
 
845
829
        tree = (BVHTree *)MEM_callocN(sizeof(BVHTree), "BVHTree");
846
830
 
847
 
        //tree epsilon must be >= FLT_EPSILON
848
 
        //so that tangent rays can still hit a bounding volume..
849
 
        //this bug would show up when casting a ray aligned with a kdop-axis and with an edge of 2 faces
850
 
        epsilon = MAX2(FLT_EPSILON, epsilon);
851
 
        
 
831
        /* tree epsilon must be >= FLT_EPSILON
 
832
         * so that tangent rays can still hit a bounding volume..
 
833
         * this bug would show up when casting a ray aligned with a kdop-axis and with an edge of 2 faces */
 
834
        epsilon = max_ff(FLT_EPSILON, epsilon);
 
835
 
852
836
        if (tree) {
853
837
                tree->epsilon = epsilon;
854
838
                tree->tree_type = tree_type; 
880
864
                }
881
865
 
882
866
 
883
 
                //Allocate arrays
 
867
                /* Allocate arrays */
884
868
                numnodes = maxsize + implicit_needed_branches(tree_type, maxsize) + tree_type;
885
869
 
886
 
                tree->nodes = (BVHNode **)MEM_callocN(sizeof(BVHNode *)*numnodes, "BVHNodes");
 
870
                tree->nodes = (BVHNode **)MEM_callocN(sizeof(BVHNode *) * numnodes, "BVHNodes");
887
871
                
888
 
                if (!tree->nodes)
889
 
                {
 
872
                if (!tree->nodes) {
890
873
                        MEM_freeN(tree);
891
874
                        return NULL;
892
875
                }
893
876
                
894
 
                tree->nodebv = (float*)MEM_callocN(sizeof(float)* axis * numnodes, "BVHNodeBV");
895
 
                if (!tree->nodebv)
896
 
                {
 
877
                tree->nodebv = (float *)MEM_callocN(sizeof(float) * axis * numnodes, "BVHNodeBV");
 
878
                if (!tree->nodebv) {
897
879
                        MEM_freeN(tree->nodes);
898
880
                        MEM_freeN(tree);
899
881
                }
900
882
 
901
 
                tree->nodechild = (BVHNode**)MEM_callocN(sizeof(BVHNode*) * tree_type * numnodes, "BVHNodeBV");
902
 
                if (!tree->nodechild)
903
 
                {
 
883
                tree->nodechild = (BVHNode **)MEM_callocN(sizeof(BVHNode *) * tree_type * numnodes, "BVHNodeBV");
 
884
                if (!tree->nodechild) {
904
885
                        MEM_freeN(tree->nodebv);
905
886
                        MEM_freeN(tree->nodes);
906
887
                        MEM_freeN(tree);
907
888
                }
908
889
 
909
 
                tree->nodearray = (BVHNode *)MEM_callocN(sizeof(BVHNode)* numnodes, "BVHNodeArray");
 
890
                tree->nodearray = (BVHNode *)MEM_callocN(sizeof(BVHNode) * numnodes, "BVHNodeArray");
910
891
                
911
 
                if (!tree->nodearray)
912
 
                {
 
892
                if (!tree->nodearray) {
913
893
                        MEM_freeN(tree->nodechild);
914
894
                        MEM_freeN(tree->nodebv);
915
895
                        MEM_freeN(tree->nodes);
917
897
                        return NULL;
918
898
                }
919
899
 
920
 
                //link the dynamic bv and child links
921
 
                for (i=0; i< numnodes; i++)
922
 
                {
 
900
                /* link the dynamic bv and child links */
 
901
                for (i = 0; i < numnodes; i++) {
923
902
                        tree->nodearray[i].bv = tree->nodebv + i * axis;
924
903
                        tree->nodearray[i].children = tree->nodechild + i * tree_type;
925
904
                }
931
910
 
932
911
void BLI_bvhtree_free(BVHTree *tree)
933
912
{       
934
 
        if (tree)
935
 
        {
 
913
        if (tree) {
936
914
                MEM_freeN(tree->nodes);
937
915
                MEM_freeN(tree->nodearray);
938
916
                MEM_freeN(tree->nodebv);
945
923
{
946
924
        int i;
947
925
 
948
 
        BVHNode*  branches_array = tree->nodearray + tree->totleaf;
949
 
        BVHNode** leafs_array    = tree->nodes;
 
926
        BVHNode *branches_array = tree->nodearray + tree->totleaf;
 
927
        BVHNode **leafs_array    = tree->nodes;
950
928
 
951
 
        //This function should only be called once (some big bug goes here if its being called more than once per tree)
 
929
        /* This function should only be called once (some big bug goes here if its being called more than once per tree) */
952
930
        assert(tree->totbranch == 0);
953
931
 
954
 
        //Build the implicit tree
 
932
        /* Build the implicit tree */
955
933
        non_recursive_bvh_div_nodes(tree, branches_array, leafs_array, tree->totleaf);
956
934
 
957
 
        //current code expects the branches to be linked to the nodes array
958
 
        //we perform that linkage here
 
935
        /*current code expects the branches to be linked to the nodes array
 
936
         *we perform that linkage here */
959
937
        tree->totbranch = implicit_needed_branches(tree->tree_type, tree->totleaf);
960
938
        for (i = 0; i < tree->totbranch; i++)
961
939
                tree->nodes[tree->totleaf + i] = branches_array + i;
962
940
 
963
941
        build_skip_links(tree, tree->nodes[tree->totleaf], NULL, NULL);
964
 
        //bvhtree_info(tree);
 
942
        /* bvhtree_info(tree); */
965
943
}
966
944
 
967
 
int BLI_bvhtree_insert(BVHTree *tree, int index, const float *co, int numpoints)
 
945
int BLI_bvhtree_insert(BVHTree *tree, int index, const float co[3], int numpoints)
968
946
{
969
 
        int i;
 
947
        axis_t axis_iter;
970
948
        BVHNode *node = NULL;
971
 
        
972
 
        // insert should only possible as long as tree->totbranch is 0
 
949
 
 
950
        /* insert should only possible as long as tree->totbranch is 0 */
973
951
        if (tree->totbranch > 0)
974
952
                return 0;
975
 
        
976
 
        if (tree->totleaf+1 >= MEM_allocN_len(tree->nodes)/sizeof(*(tree->nodes)))
 
953
 
 
954
        if (tree->totleaf + 1 >= MEM_allocN_len(tree->nodes) / sizeof(*(tree->nodes)))
977
955
                return 0;
978
 
        
979
 
        // TODO check if have enough nodes in array
980
 
        
 
956
 
 
957
        /* TODO check if have enough nodes in array */
 
958
 
981
959
        node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]);
982
960
        tree->totleaf++;
983
 
        
 
961
 
984
962
        create_kdop_hull(tree, node, co, numpoints, 0);
985
 
        node->index= index;
986
 
        
987
 
        // inflate the bv with some epsilon
988
 
        for (i = tree->start_axis; i < tree->stop_axis; i++)
989
 
        {
990
 
                node->bv[(2 * i)] -= tree->epsilon; // minimum 
991
 
                node->bv[(2 * i) + 1] += tree->epsilon; // maximum 
 
963
        node->index = index;
 
964
 
 
965
        /* inflate the bv with some epsilon */
 
966
        for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
 
967
                node->bv[(2 * axis_iter)] -= tree->epsilon; /* minimum */
 
968
                node->bv[(2 * axis_iter) + 1] += tree->epsilon; /* maximum */
992
969
        }
993
970
 
994
971
        return 1;
995
972
}
996
973
 
997
974
 
998
 
// call before BLI_bvhtree_update_tree()
999
 
int BLI_bvhtree_update_node(BVHTree *tree, int index, const float *co, const float *co_moving, int numpoints)
 
975
/* call before BLI_bvhtree_update_tree() */
 
976
int BLI_bvhtree_update_node(BVHTree *tree, int index, const float co[3], const float co_moving[3], int numpoints)
1000
977
{
1001
 
        int i;
1002
 
        BVHNode *node= NULL;
 
978
        BVHNode *node = NULL;
 
979
        axis_t axis_iter;
1003
980
        
1004
 
        // check if index exists
 
981
        /* check if index exists */
1005
982
        if (index > tree->totleaf)
1006
983
                return 0;
1007
984
        
1012
989
        if (co_moving)
1013
990
                create_kdop_hull(tree, node, co_moving, numpoints, 1);
1014
991
        
1015
 
        // inflate the bv with some epsilon
1016
 
        for (i = tree->start_axis; i < tree->stop_axis; i++)
1017
 
        {
1018
 
                node->bv[(2 * i)] -= tree->epsilon; // minimum 
1019
 
                node->bv[(2 * i) + 1] += tree->epsilon; // maximum 
 
992
        /* inflate the bv with some epsilon */
 
993
        for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) {
 
994
                node->bv[(2 * axis_iter)]     -= tree->epsilon; /* minimum */
 
995
                node->bv[(2 * axis_iter) + 1] += tree->epsilon; /* maximum */
1020
996
        }
1021
997
 
1022
998
        return 1;
1023
999
}
1024
1000
 
1025
 
// call BLI_bvhtree_update_node() first for every node/point/triangle
 
1001
/* call BLI_bvhtree_update_node() first for every node/point/triangle */
1026
1002
void BLI_bvhtree_update_tree(BVHTree *tree)
1027
1003
{
1028
 
        //Update bottom=>top
1029
 
        //TRICKY: the way we build the tree all the childs have an index greater than the parent
1030
 
        //This allows us todo a bottom up update by starting on the biger numbered branch
 
1004
        /* Update bottom=>top
 
1005
         * TRICKY: the way we build the tree all the childs have an index greater than the parent
 
1006
         * This allows us todo a bottom up update by starting on the bigger numbered branch */
1031
1007
 
1032
 
        BVHNode** root  = tree->nodes + tree->totleaf;
1033
 
        BVHNode** index = tree->nodes + tree->totleaf + tree->totbranch-1;
 
1008
        BVHNode **root  = tree->nodes + tree->totleaf;
 
1009
        BVHNode **index = tree->nodes + tree->totleaf + tree->totbranch - 1;
1034
1010
 
1035
1011
        for (; index >= root; index--)
1036
1012
                node_join(tree, *index);
1037
1013
}
1038
1014
 
1039
 
float BLI_bvhtree_getepsilon(BVHTree *tree)
 
1015
float BLI_bvhtree_getepsilon(const BVHTree *tree)
1040
1016
{
1041
1017
        return tree->epsilon;
1042
1018
}
1044
1020
 
1045
1021
/*
1046
1022
 * BLI_bvhtree_overlap
1047
 
 */
1048
 
// overlap - is it possbile for 2 bv's to collide ?
1049
 
static int tree_overlap(BVHNode *node1, BVHNode *node2, int start_axis, int stop_axis)
 
1023
 *
 
1024
 * overlap - is it possible for 2 bv's to collide ? */
 
1025
static int tree_overlap(BVHNode *node1, BVHNode *node2, axis_t start_axis, axis_t stop_axis)
1050
1026
{
1051
1027
        float *bv1 = node1->bv;
1052
1028
        float *bv2 = node2->bv;
1053
1029
 
1054
 
        float *bv1_end = bv1 + (stop_axis<<1);
 
1030
        float *bv1_end = bv1 + (stop_axis << 1);
1055
1031
                
1056
 
        bv1 += start_axis<<1;
1057
 
        bv2 += start_axis<<1;
 
1032
        bv1 += start_axis << 1;
 
1033
        bv2 += start_axis << 1;
1058
1034
        
1059
 
        // test all axis if min + max overlap
1060
 
        for (; bv1 != bv1_end; bv1+=2, bv2+=2)
1061
 
        {
1062
 
                if ((*(bv1) > *(bv2 + 1)) || (*(bv2) > *(bv1 + 1))) 
 
1035
        /* test all axis if min + max overlap */
 
1036
        for (; bv1 != bv1_end; bv1 += 2, bv2 += 2) {
 
1037
                if ((*(bv1) > *(bv2 + 1)) || (*(bv2) > *(bv1 + 1)))
1063
1038
                        return 0;
1064
1039
        }
1065
 
        
 
1040
 
1066
1041
        return 1;
1067
1042
}
1068
1043
 
1069
1044
static void traverse(BVHOverlapData *data, BVHNode *node1, BVHNode *node2)
1070
1045
{
1071
1046
        int j;
1072
 
        
1073
 
        if (tree_overlap(node1, node2, data->start_axis, data->stop_axis))
1074
 
        {
1075
 
                // check if node1 is a leaf
1076
 
                if (!node1->totnode)
1077
 
                {
1078
 
                        // check if node2 is a leaf
1079
 
                        if (!node2->totnode)
1080
 
                        {
1081
 
                                
1082
 
                                if (node1 == node2)
1083
 
                                {
 
1047
 
 
1048
        if (tree_overlap(node1, node2, data->start_axis, data->stop_axis)) {
 
1049
                /* check if node1 is a leaf */
 
1050
                if (!node1->totnode) {
 
1051
                        /* check if node2 is a leaf */
 
1052
                        if (!node2->totnode) {
 
1053
 
 
1054
                                if (node1 == node2) {
1084
1055
                                        return;
1085
1056
                                }
1086
 
                                        
1087
 
                                if (data->i >= data->max_overlap)
1088
 
                                {       
1089
 
                                        // try to make alloc'ed memory bigger
1090
 
                                        data->overlap = realloc(data->overlap, sizeof(BVHTreeOverlap)*data->max_overlap*2);
1091
 
                                        
1092
 
                                        if (!data->overlap)
1093
 
                                        {
 
1057
 
 
1058
                                if (data->i >= data->max_overlap) {
 
1059
                                        /* try to make alloc'ed memory bigger */
 
1060
                                        data->overlap = realloc(data->overlap, sizeof(BVHTreeOverlap) * data->max_overlap * 2);
 
1061
                                        
 
1062
                                        if (!data->overlap) {
1094
1063
                                                printf("Out of Memory in traverse\n");
1095
1064
                                                return;
1096
1065
                                        }
1097
1066
                                        data->max_overlap *= 2;
1098
1067
                                }
1099
1068
                                
1100
 
                                // both leafs, insert overlap!
 
1069
                                /* both leafs, insert overlap! */
1101
1070
                                data->overlap[data->i].indexA = node1->index;
1102
1071
                                data->overlap[data->i].indexB = node2->index;
1103
1072
 
1104
1073
                                data->i++;
1105
1074
                        }
1106
 
                        else
1107
 
                        {
1108
 
                                for (j = 0; j < data->tree2->tree_type; j++)
1109
 
                                {
 
1075
                        else {
 
1076
                                for (j = 0; j < data->tree2->tree_type; j++) {
1110
1077
                                        if (node2->children[j])
1111
1078
                                                traverse(data, node1, node2->children[j]);
1112
1079
                                }
1113
1080
                        }
1114
1081
                }
1115
 
                else
1116
 
                {
1117
 
                        
1118
 
                        for (j = 0; j < data->tree2->tree_type; j++)
1119
 
                        {
 
1082
                else {
 
1083
                        for (j = 0; j < data->tree2->tree_type; j++) {
1120
1084
                                if (node1->children[j])
1121
1085
                                        traverse(data, node1->children[j], node2);
1122
1086
                        }
1132
1096
        BVHTreeOverlap *overlap = NULL, *to = NULL;
1133
1097
        BVHOverlapData **data;
1134
1098
        
1135
 
        // check for compatibility of both trees (can't compare 14-DOP with 18-DOP)
 
1099
        /* check for compatibility of both trees (can't compare 14-DOP with 18-DOP) */
1136
1100
        if ((tree1->axis != tree2->axis) && (tree1->axis == 14 || tree2->axis == 14) && (tree1->axis == 18 || tree2->axis == 18))
1137
1101
                return NULL;
1138
1102
        
1139
 
        // fast check root nodes for collision before doing big splitting + traversal
1140
 
        if (!tree_overlap(tree1->nodes[tree1->totleaf], tree2->nodes[tree2->totleaf], MIN2(tree1->start_axis, tree2->start_axis), MIN2(tree1->stop_axis, tree2->stop_axis)))
 
1103
        /* fast check root nodes for collision before doing big splitting + traversal */
 
1104
        if (!tree_overlap(tree1->nodes[tree1->totleaf], tree2->nodes[tree2->totleaf],
 
1105
                          min_axis(tree1->start_axis, tree2->start_axis),
 
1106
                          min_axis(tree1->stop_axis, tree2->stop_axis)))
 
1107
        {
1141
1108
                return NULL;
 
1109
        }
1142
1110
 
1143
 
        data = MEM_callocN(sizeof(BVHOverlapData *)* tree1->tree_type, "BVHOverlapData_star");
 
1111
        data = MEM_callocN(sizeof(BVHOverlapData *) * tree1->tree_type, "BVHOverlapData_star");
1144
1112
        
1145
 
        for (j = 0; j < tree1->tree_type; j++)
1146
 
        {
 
1113
        for (j = 0; j < tree1->tree_type; j++) {
1147
1114
                data[j] = (BVHOverlapData *)MEM_callocN(sizeof(BVHOverlapData), "BVHOverlapData");
1148
1115
                
1149
 
                // init BVHOverlapData
1150
 
                data[j]->overlap = (BVHTreeOverlap *)malloc(sizeof(BVHTreeOverlap)*MAX2(tree1->totleaf, tree2->totleaf));
 
1116
                /* init BVHOverlapData */
 
1117
                data[j]->overlap = (BVHTreeOverlap *)malloc(sizeof(BVHTreeOverlap) * max_ii(tree1->totleaf, tree2->totleaf));
1151
1118
                data[j]->tree1 = tree1;
1152
1119
                data[j]->tree2 = tree2;
1153
 
                data[j]->max_overlap = MAX2(tree1->totleaf, tree2->totleaf);
 
1120
                data[j]->max_overlap = max_ii(tree1->totleaf, tree2->totleaf);
1154
1121
                data[j]->i = 0;
1155
 
                data[j]->start_axis = MIN2(tree1->start_axis, tree2->start_axis);
1156
 
                data[j]->stop_axis  = MIN2(tree1->stop_axis,  tree2->stop_axis );
 
1122
                data[j]->start_axis = min_axis(tree1->start_axis, tree2->start_axis);
 
1123
                data[j]->stop_axis  = min_axis(tree1->stop_axis,  tree2->stop_axis);
1157
1124
        }
1158
1125
 
1159
1126
#pragma omp parallel for private(j) schedule(static)
1160
 
        for (j = 0; j < MIN2(tree1->tree_type, tree1->nodes[tree1->totleaf]->totnode); j++)
1161
 
        {
 
1127
        for (j = 0; j < MIN2(tree1->tree_type, tree1->nodes[tree1->totleaf]->totnode); j++) {
1162
1128
                traverse(data[j], tree1->nodes[tree1->totleaf]->children[j], tree2->nodes[tree2->totleaf]);
1163
1129
        }
1164
1130
        
1165
1131
        for (j = 0; j < tree1->tree_type; j++)
1166
1132
                total += data[j]->i;
1167
1133
        
1168
 
        to = overlap = (BVHTreeOverlap *)MEM_callocN(sizeof(BVHTreeOverlap)*total, "BVHTreeOverlap");
 
1134
        to = overlap = (BVHTreeOverlap *)MEM_callocN(sizeof(BVHTreeOverlap) * total, "BVHTreeOverlap");
1169
1135
        
1170
 
        for (j = 0; j < tree1->tree_type; j++)
1171
 
        {
1172
 
                memcpy(to, data[j]->overlap, data[j]->i*sizeof(BVHTreeOverlap));
1173
 
                to+=data[j]->i;
 
1136
        for (j = 0; j < tree1->tree_type; j++) {
 
1137
                memcpy(to, data[j]->overlap, data[j]->i * sizeof(BVHTreeOverlap));
 
1138
                to += data[j]->i;
1174
1139
        }
1175
1140
        
1176
 
        for (j = 0; j < tree1->tree_type; j++)
1177
 
        {
 
1141
        for (j = 0; j < tree1->tree_type; j++) {
1178
1142
                free(data[j]->overlap);
1179
1143
                MEM_freeN(data[j]);
1180
1144
        }
1184
1148
        return overlap;
1185
1149
}
1186
1150
 
1187
 
//Determines the nearest point of the given node BV. Returns the squared distance to that point.
 
1151
/* Determines the nearest point of the given node BV. Returns the squared distance to that point. */
1188
1152
static float calc_nearest_point(const float proj[3], BVHNode *node, float *nearest)
1189
1153
{
1190
1154
        int i;
1191
1155
        const float *bv = node->bv;
1192
1156
 
1193
 
        //nearest on AABB hull
1194
 
        for (i=0; i != 3; i++, bv += 2)
1195
 
        {
 
1157
        /* nearest on AABB hull */
 
1158
        for (i = 0; i != 3; i++, bv += 2) {
1196
1159
                if (bv[0] > proj[i])
1197
1160
                        nearest[i] = bv[0];
1198
1161
                else if (bv[1] < proj[i])
1202
1165
        }
1203
1166
 
1204
1167
#if 0
1205
 
        //nearest on a general hull
 
1168
        /* nearest on a general hull */
1206
1169
        copy_v3_v3(nearest, data->co);
1207
 
        for (i = data->tree->start_axis; i != data->tree->stop_axis; i++, bv+=2)
 
1170
        for (i = data->tree->start_axis; i != data->tree->stop_axis; i++, bv += 2)
1208
1171
        {
1209
 
                float proj = dot_v3v3( nearest, KDOP_AXES[i]);
 
1172
                float proj = dot_v3v3(nearest, KDOP_AXES[i]);
1210
1173
                float dl = bv[0] - proj;
1211
1174
                float du = bv[1] - proj;
1212
1175
 
1223
1186
}
1224
1187
 
1225
1188
 
1226
 
typedef struct NodeDistance
1227
 
{
 
1189
typedef struct NodeDistance {
1228
1190
        BVHNode *node;
1229
1191
        float dist;
1230
1192
 
1231
1193
} NodeDistance;
1232
1194
 
1233
 
#define NodeDistance_priority(a,b)      ( (a).dist < (b).dist )
1234
 
 
1235
 
// TODO: use a priority queue to reduce the number of nodes looked on
 
1195
/* TODO: use a priority queue to reduce the number of nodes looked on */
1236
1196
static void dfs_find_nearest_dfs(BVHNearestData *data, BVHNode *node)
1237
1197
{
1238
 
        if (node->totnode == 0)
1239
 
        {
 
1198
        if (node->totnode == 0) {
1240
1199
                if (data->callback)
1241
 
                        data->callback(data->userdata , node->index, data->co, &data->nearest);
1242
 
                else
1243
 
                {
1244
 
                        data->nearest.index     = node->index;
1245
 
                        data->nearest.dist      = calc_nearest_point(data->proj, node, data->nearest.co);
 
1200
                        data->callback(data->userdata, node->index, data->co, &data->nearest);
 
1201
                else {
 
1202
                        data->nearest.index = node->index;
 
1203
                        data->nearest.dist  = calc_nearest_point(data->proj, node, data->nearest.co);
1246
1204
                }
1247
1205
        }
1248
 
        else
1249
 
        {
1250
 
                //Better heuristic to pick the closest node to dive on
 
1206
        else {
 
1207
                /* Better heuristic to pick the closest node to dive on */
1251
1208
                int i;
1252
1209
                float nearest[3];
1253
1210
 
1254
 
                if (data->proj[ node->main_axis ] <= node->children[0]->bv[node->main_axis*2+1])
1255
 
                {
 
1211
                if (data->proj[node->main_axis] <= node->children[0]->bv[node->main_axis * 2 + 1]) {
1256
1212
 
1257
 
                        for (i=0; i != node->totnode; i++)
1258
 
                        {
1259
 
                                if ( calc_nearest_point(data->proj, node->children[i], nearest) >= data->nearest.dist) continue;
 
1213
                        for (i = 0; i != node->totnode; i++) {
 
1214
                                if (calc_nearest_point(data->proj, node->children[i], nearest) >= data->nearest.dist) continue;
1260
1215
                                dfs_find_nearest_dfs(data, node->children[i]);
1261
1216
                        }
1262
1217
                }
1263
 
                else
1264
 
                {
1265
 
                        for (i=node->totnode-1; i >= 0 ; i--)
1266
 
                        {
1267
 
                                if ( calc_nearest_point(data->proj, node->children[i], nearest) >= data->nearest.dist) continue;
 
1218
                else {
 
1219
                        for (i = node->totnode - 1; i >= 0; i--) {
 
1220
                                if (calc_nearest_point(data->proj, node->children[i], nearest) >= data->nearest.dist) continue;
1268
1221
                                dfs_find_nearest_dfs(data, node->children[i]);
1269
1222
                        }
1270
1223
                }
1281
1234
 
1282
1235
 
1283
1236
#if 0
 
1237
 
 
1238
#define DEFAULT_FIND_NEAREST_HEAP_SIZE 1024
 
1239
 
 
1240
#define NodeDistance_priority(a, b) ((a).dist < (b).dist)
 
1241
 
1284
1242
static void NodeDistance_push_heap(NodeDistance *heap, int heap_size)
1285
1243
PUSH_HEAP_BODY(NodeDistance, NodeDistance_priority, heap, heap_size)
1286
1244
 
1287
1245
static void NodeDistance_pop_heap(NodeDistance *heap, int heap_size)
1288
1246
POP_HEAP_BODY(NodeDistance, NodeDistance_priority, heap, heap_size)
1289
1247
 
1290
 
//NN function that uses an heap.. this functions leads to an optimal number of min-distance
1291
 
//but for normal tri-faces and BV 6-dop.. a simple dfs with local heuristics (as implemented
1292
 
//in source/blender/blenkernel/intern/shrinkwrap.c) works faster.
1293
 
//
1294
 
//It may make sense to use this function if the callback queries are very slow.. or if its impossible
1295
 
//to get a nice heuristic
1296
 
//
1297
 
//this function uses "malloc/free" instead of the MEM_* because it intends to be openmp safe
 
1248
/* NN function that uses an heap.. this functions leads to an optimal number of min-distance
 
1249
 * but for normal tri-faces and BV 6-dop.. a simple dfs with local heuristics (as implemented
 
1250
 * in source/blender/blenkernel/intern/shrinkwrap.c) works faster.
 
1251
 *
 
1252
 * It may make sense to use this function if the callback queries are very slow.. or if its impossible
 
1253
 * to get a nice heuristic
 
1254
 *
 
1255
 * this function uses "malloc/free" instead of the MEM_* because it intends to be openmp safe */
1298
1256
static void bfs_find_nearest(BVHNearestData *data, BVHNode *node)
1299
1257
{
1300
1258
        int i;
1301
1259
        NodeDistance default_heap[DEFAULT_FIND_NEAREST_HEAP_SIZE];
1302
 
        NodeDistance *heap=default_heap, current;
1303
 
        int heap_size = 0, max_heap_size = sizeof(default_heap)/sizeof(default_heap[0]);
 
1260
        NodeDistance *heap = default_heap, current;
 
1261
        int heap_size = 0, max_heap_size = sizeof(default_heap) / sizeof(default_heap[0]);
1304
1262
        float nearest[3];
1305
1263
 
1306
1264
        int callbacks = 0, push_heaps = 0;
1317
1275
        while (current.dist < data->nearest.dist)
1318
1276
        {
1319
1277
//              printf("%f : %f\n", current.dist, data->nearest.dist);
1320
 
                for (i=0; i< current.node->totnode; i++)
 
1278
                for (i = 0; i < current.node->totnode; i++)
1321
1279
                {
1322
1280
                        BVHNode *child = current.node->children[i];
1323
1281
                        if (child->totnode == 0)
1325
1283
                                callbacks++;
1326
1284
                                dfs_find_nearest_dfs(data, child);
1327
1285
                        }
1328
 
                        else
1329
 
                        {
1330
 
                                //adjust heap size
 
1286
                        else {
 
1287
                                /* adjust heap size */
1331
1288
                                if ((heap_size >= max_heap_size) &&
1332
 
                                    ADJUST_MEMORY(default_heap, (void**)&heap, heap_size+1, &max_heap_size, sizeof(heap[0])) == FALSE)
 
1289
                                    ADJUST_MEMORY(default_heap, (void **)&heap, heap_size + 1, &max_heap_size, sizeof(heap[0])) == FALSE)
1333
1290
                                {
1334
1291
                                        printf("WARNING: bvh_find_nearest got out of memory\n");
1335
1292
 
1346
1303
                                heap_size++;
1347
1304
 
1348
1305
                                NodeDistance_push_heap(heap, heap_size);
1349
 
        //                      PUSH_HEAP_BODY(NodeDistance, NodeDistance_priority, heap, heap_size);
 
1306
                                //                      PUSH_HEAP_BODY(NodeDistance, NodeDistance_priority, heap, heap_size);
1350
1307
                                push_heaps++;
1351
1308
                        }
1352
1309
                }
1367
1324
#endif
1368
1325
 
1369
1326
 
1370
 
int BLI_bvhtree_find_nearest(BVHTree *tree, const float co[3], BVHTreeNearest *nearest, BVHTree_NearestPointCallback callback, void *userdata)
 
1327
int BLI_bvhtree_find_nearest(BVHTree *tree, const float co[3], BVHTreeNearest *nearest,
 
1328
                             BVHTree_NearestPointCallback callback, void *userdata)
1371
1329
{
1372
 
        int i;
 
1330
        axis_t axis_iter;
1373
1331
 
1374
1332
        BVHNearestData data;
1375
 
        BVHNode* root = tree->nodes[tree->totleaf];
 
1333
        BVHNode *root = tree->nodes[tree->totleaf];
1376
1334
 
1377
 
        //init data to search
 
1335
        /* init data to search */
1378
1336
        data.tree = tree;
1379
1337
        data.co = co;
1380
1338
 
1381
1339
        data.callback = callback;
1382
1340
        data.userdata = userdata;
1383
1341
 
1384
 
        for (i = data.tree->start_axis; i != data.tree->stop_axis; i++)
1385
 
        {
1386
 
                data.proj[i] = dot_v3v3(data.co, KDOP_AXES[i]);
 
1342
        for (axis_iter = data.tree->start_axis; axis_iter != data.tree->stop_axis; axis_iter++) {
 
1343
                data.proj[axis_iter] = dot_v3v3(data.co, KDOP_AXES[axis_iter]);
1387
1344
        }
1388
1345
 
1389
 
        if (nearest)
1390
 
        {
1391
 
                memcpy( &data.nearest , nearest, sizeof(*nearest) );
 
1346
        if (nearest) {
 
1347
                memcpy(&data.nearest, nearest, sizeof(*nearest));
1392
1348
        }
1393
 
        else
1394
 
        {
 
1349
        else {
1395
1350
                data.nearest.index = -1;
1396
1351
                data.nearest.dist = FLT_MAX;
1397
1352
        }
1398
1353
 
1399
 
        //dfs search
 
1354
        /* dfs search */
1400
1355
        if (root)
1401
1356
                dfs_find_nearest_begin(&data, root);
1402
1357
 
1403
 
        //copy back results
1404
 
        if (nearest)
1405
 
        {
 
1358
        /* copy back results */
 
1359
        if (nearest) {
1406
1360
                memcpy(nearest, &data.nearest, sizeof(*nearest));
1407
1361
        }
1408
1362
 
1417
1371
 */
1418
1372
 
1419
1373
 
1420
 
//Determines the distance that the ray must travel to hit the bounding volume of the given node
1421
 
static float ray_nearest_hit(BVHRayCastData *data, float *bv)
 
1374
/* Determines the distance that the ray must travel to hit the bounding volume of the given node */
 
1375
static float ray_nearest_hit(BVHRayCastData *data, const float bv[6])
1422
1376
{
1423
1377
        int i;
1424
1378
 
1425
1379
        float low = 0, upper = data->hit.dist;
1426
1380
 
1427
 
        for (i=0; i != 3; i++, bv += 2)
1428
 
        {
1429
 
                if (data->ray_dot_axis[i] == 0.0f)
1430
 
                {
1431
 
                        //axis aligned ray
 
1381
        for (i = 0; i != 3; i++, bv += 2) {
 
1382
                if (data->ray_dot_axis[i] == 0.0f) {
 
1383
                        /* axis aligned ray */
1432
1384
                        if (data->ray.origin[i] < bv[0] - data->ray.radius ||
1433
1385
                            data->ray.origin[i] > bv[1] + data->ray.radius)
1434
1386
                        {
1435
1387
                                return FLT_MAX;
1436
1388
                        }
1437
1389
                }
1438
 
                else
1439
 
                {
 
1390
                else {
1440
1391
                        float ll = (bv[0] - data->ray.radius - data->ray.origin[i]) / data->ray_dot_axis[i];
1441
1392
                        float lu = (bv[1] + data->ray.radius - data->ray.origin[i]) / data->ray_dot_axis[i];
1442
1393
 
1443
 
                        if (data->ray_dot_axis[i] > 0.0f)
1444
 
                        {
1445
 
                                if (ll > low)   low = ll;
 
1394
                        if (data->ray_dot_axis[i] > 0.0f) {
 
1395
                                if (ll > low) low = ll;
1446
1396
                                if (lu < upper) upper = lu;
1447
1397
                        }
1448
 
                        else
1449
 
                        {
1450
 
                                if (lu > low)   low = lu;
 
1398
                        else {
 
1399
                                if (lu > low) low = lu;
1451
1400
                                if (ll < upper) upper = ll;
1452
1401
                        }
1453
1402
        
1457
1406
        return low;
1458
1407
}
1459
1408
 
1460
 
//Determines the distance that the ray must travel to hit the bounding volume of the given node
1461
 
//Based on Tactical Optimization of Ray/Box Intersection, by Graham Fyffe
1462
 
//[http://tog.acm.org/resources/RTNews/html/rtnv21n1.html#art9]
1463
 
//
1464
 
//TODO this doens't has data->ray.radius in consideration
 
1409
/* Determines the distance that the ray must travel to hit the bounding volume of the given node
 
1410
 * Based on Tactical Optimization of Ray/Box Intersection, by Graham Fyffe
 
1411
 * [http://tog.acm.org/resources/RTNews/html/rtnv21n1.html#art9]
 
1412
 *
 
1413
 * TODO this doesn't take data->ray.radius into consideration */
1465
1414
static float fast_ray_nearest_hit(const BVHRayCastData *data, const BVHNode *node)
1466
1415
{
1467
1416
        const float *bv = node->bv;
1488
1437
{
1489
1438
        int i;
1490
1439
 
1491
 
        //ray-bv is really fast.. and simple tests revealed its worth to test it
1492
 
        //before calling the ray-primitive functions
 
1440
        /* ray-bv is really fast.. and simple tests revealed its worth to test it
 
1441
         * before calling the ray-primitive functions */
1493
1442
        /* XXX: temporary solution for particles until fast_ray_nearest_hit supports ray.radius */
1494
1443
        float dist = (data->ray.radius > 0.0f) ? ray_nearest_hit(data, node->bv) : fast_ray_nearest_hit(data, node);
1495
1444
        if (dist >= data->hit.dist) return;
1496
1445
 
1497
 
        if (node->totnode == 0)
1498
 
        {
1499
 
                if (data->callback)
 
1446
        if (node->totnode == 0) {
 
1447
                if (data->callback) {
1500
1448
                        data->callback(data->userdata, node->index, &data->ray, &data->hit);
1501
 
                else
1502
 
                {
1503
 
                        data->hit.index = node->index;
 
1449
                }
 
1450
                else {
 
1451
                        data->hit.index = node->index;
1504
1452
                        data->hit.dist  = dist;
1505
1453
                        madd_v3_v3v3fl(data->hit.co, data->ray.origin, data->ray.direction, dist);
1506
1454
                }
1507
1455
        }
1508
 
        else
1509
 
        {
1510
 
                //pick loop direction to dive into the tree (based on ray direction and split axis)
1511
 
                if (data->ray_dot_axis[ (int)node->main_axis ] > 0.0f)
1512
 
                {
1513
 
                        for (i=0; i != node->totnode; i++)
1514
 
                        {
 
1456
        else {
 
1457
                /* pick loop direction to dive into the tree (based on ray direction and split axis) */
 
1458
                if (data->ray_dot_axis[(int)node->main_axis] > 0.0f) {
 
1459
                        for (i = 0; i != node->totnode; i++) {
1515
1460
                                dfs_raycast(data, node->children[i]);
1516
1461
                        }
1517
1462
                }
1518
 
                else
1519
 
                {
1520
 
                        for (i=node->totnode-1; i >= 0; i--)
1521
 
                        {
 
1463
                else {
 
1464
                        for (i = node->totnode - 1; i >= 0; i--) {
1522
1465
                                dfs_raycast(data, node->children[i]);
1523
1466
                        }
1524
1467
                }
1539
1482
 
1540
1483
                if (node->totnode == 0)
1541
1484
                {
1542
 
                        if (data->callback)
 
1485
                        if (data->callback) {
1543
1486
                                data->callback(data->userdata, node->index, &data->ray, &data->hit);
1544
 
                        else
1545
 
                        {
1546
 
                                data->hit.index = node->index;
 
1487
                        }
 
1488
                        else {
 
1489
                                data->hit.index = node->index;
1547
1490
                                data->hit.dist  = dist;
1548
1491
                                madd_v3_v3v3fl(data->hit.co, data->ray.origin, data->ray.direction, dist);
1549
1492
                        }
1550
1493
                        
1551
1494
                        node = node->skip[1];
1552
1495
                }
1553
 
                else
1554
 
                {
 
1496
                else {
1555
1497
                        node = node->children[0];
1556
 
                }       
 
1498
                }
1557
1499
        }
1558
1500
}
1559
1501
#endif
1560
1502
 
1561
 
int BLI_bvhtree_ray_cast(BVHTree *tree, const float co[3], const float dir[3], float radius, BVHTreeRayHit *hit, BVHTree_RayCastCallback callback, void *userdata)
 
1503
int BLI_bvhtree_ray_cast(BVHTree *tree, const float co[3], const float dir[3], float radius, BVHTreeRayHit *hit,
 
1504
                         BVHTree_RayCastCallback callback, void *userdata)
1562
1505
{
1563
1506
        int i;
1564
1507
        BVHRayCastData data;
1565
 
        BVHNode * root = tree->nodes[tree->totleaf];
 
1508
        BVHNode *root = tree->nodes[tree->totleaf];
1566
1509
 
1567
1510
        data.tree = tree;
1568
1511
 
1575
1518
 
1576
1519
        normalize_v3(data.ray.direction);
1577
1520
 
1578
 
        for (i=0; i<3; i++)
1579
 
        {
 
1521
        for (i = 0; i < 3; i++) {
1580
1522
                data.ray_dot_axis[i] = dot_v3v3(data.ray.direction, KDOP_AXES[i]);
1581
1523
                data.idot_axis[i] = 1.0f / data.ray_dot_axis[i];
1582
1524
 
1583
 
                if (fabsf(data.ray_dot_axis[i]) < FLT_EPSILON)
1584
 
                {
 
1525
                if (fabsf(data.ray_dot_axis[i]) < FLT_EPSILON) {
1585
1526
                        data.ray_dot_axis[i] = 0.0;
1586
1527
                }
1587
 
                data.index[2*i] = data.idot_axis[i] < 0.0f ? 1 : 0;
1588
 
                data.index[2*i+1] = 1 - data.index[2*i];
1589
 
                data.index[2*i]   += 2*i;
1590
 
                data.index[2*i+1] += 2*i;
 
1528
                data.index[2 * i] = data.idot_axis[i] < 0.0f ? 1 : 0;
 
1529
                data.index[2 * i + 1] = 1 - data.index[2 * i];
 
1530
                data.index[2 * i]   += 2 * i;
 
1531
                data.index[2 * i + 1] += 2 * i;
1591
1532
        }
1592
1533
 
1593
1534
 
1594
1535
        if (hit)
1595
 
                memcpy( &data.hit, hit, sizeof(*hit) );
1596
 
        else
1597
 
        {
 
1536
                memcpy(&data.hit, hit, sizeof(*hit));
 
1537
        else {
1598
1538
                data.hit.index = -1;
1599
1539
                data.hit.dist = FLT_MAX;
1600
1540
        }
1601
1541
 
1602
 
        if (root)
1603
 
        {
 
1542
        if (root) {
1604
1543
                dfs_raycast(&data, root);
1605
1544
//              iterative_raycast(&data, root);
1606
 
         }
 
1545
        }
1607
1546
 
1608
1547
 
1609
1548
        if (hit)
1610
 
                memcpy( hit, &data.hit, sizeof(*hit) );
 
1549
                memcpy(hit, &data.hit, sizeof(*hit));
1611
1550
 
1612
1551
        return data.hit.index;
1613
1552
}
1614
1553
 
1615
 
float BLI_bvhtree_bb_raycast(float *bv, const float light_start[3], const float light_end[3], float pos[3])
 
1554
float BLI_bvhtree_bb_raycast(const float bv[6], const float light_start[3], const float light_end[3], float pos[3])
1616
1555
{
1617
1556
        BVHRayCastData data;
1618
 
        float dist = 0.0;
 
1557
        float dist;
1619
1558
 
1620
1559
        data.hit.dist = FLT_MAX;
1621
1560
        
1622
 
        // get light direction
 
1561
        /* get light direction */
1623
1562
        data.ray.direction[0] = light_end[0] - light_start[0];
1624
1563
        data.ray.direction[1] = light_end[1] - light_start[1];
1625
1564
        data.ray.direction[2] = light_end[2] - light_start[2];
1635
1574
        
1636
1575
        dist = ray_nearest_hit(&data, bv);
1637
1576
        
1638
 
        if (dist > 0.0f)
1639
 
        {
 
1577
        if (dist > 0.0f) {
1640
1578
                madd_v3_v3v3fl(pos, light_start, data.ray.direction, dist);
1641
1579
        }
1642
1580
        return dist;
1649
1587
 * Allocs and fills an array with the indexs of node that are on the given spherical range (center, radius) 
1650
1588
 * Returns the size of the array.
1651
1589
 */
1652
 
typedef struct RangeQueryData
1653
 
{
 
1590
typedef struct RangeQueryData {
1654
1591
        BVHTree *tree;
1655
1592
        const float *center;
1656
 
        float radius;                   //squared radius
 
1593
        float radius;  /* squared radius */
1657
1594
 
1658
1595
        int hits;
1659
1596
 
1666
1603
 
1667
1604
static void dfs_range_query(RangeQueryData *data, BVHNode *node)
1668
1605
{
1669
 
        if (node->totnode == 0)
1670
 
        {
1671
 
#if 0   /*UNUSED*/
1672
 
                //Calculate the node min-coords (if the node was a point then this is the point coordinates)
 
1606
        if (node->totnode == 0) {
 
1607
#if 0   /*UNUSED*/
 
1608
                /* Calculate the node min-coords (if the node was a point then this is the point coordinates) */
1673
1609
                float co[3];
1674
1610
                co[0] = node->bv[0];
1675
1611
                co[1] = node->bv[2];
1676
1612
                co[2] = node->bv[4];
1677
1613
#endif
1678
1614
        }
1679
 
        else
1680
 
        {
 
1615
        else {
1681
1616
                int i;
1682
 
                for (i=0; i != node->totnode; i++)
1683
 
                {
 
1617
                for (i = 0; i != node->totnode; i++) {
1684
1618
                        float nearest[3];
1685
1619
                        float dist = calc_nearest_point(data->center, node->children[i], nearest);
1686
 
                        if (dist < data->radius)
1687
 
                        {
1688
 
                                //Its a leaf.. call the callback
1689
 
                                if (node->children[i]->totnode == 0)
1690
 
                                {
 
1620
                        if (dist < data->radius) {
 
1621
                                /* Its a leaf.. call the callback */
 
1622
                                if (node->children[i]->totnode == 0) {
1691
1623
                                        data->hits++;
1692
1624
                                        data->callback(data->userdata, node->children[i]->index, dist);
1693
1625
                                }
1694
1626
                                else
1695
 
                                        dfs_range_query( data, node->children[i] );
 
1627
                                        dfs_range_query(data, node->children[i]);
1696
1628
                        }
1697
1629
                }
1698
1630
        }
1700
1632
 
1701
1633
int BLI_bvhtree_range_query(BVHTree *tree, const float co[3], float radius, BVHTree_RangeQuery callback, void *userdata)
1702
1634
{
1703
 
        BVHNode * root = tree->nodes[tree->totleaf];
 
1635
        BVHNode *root = tree->nodes[tree->totleaf];
1704
1636
 
1705
1637
        RangeQueryData data;
1706
1638
        data.tree = tree;
1707
1639
        data.center = co;
1708
 
        data.radius = radius*radius;
 
1640
        data.radius = radius * radius;
1709
1641
        data.hits = 0;
1710
1642
 
1711
1643
        data.callback = callback;
1712
1644
        data.userdata = userdata;
1713
1645
 
1714
 
        if (root != NULL)
1715
 
        {
 
1646
        if (root != NULL) {
1716
1647
                float nearest[3];
1717
1648
                float dist = calc_nearest_point(data.center, root, nearest);
1718
 
                if (dist < data.radius)
1719
 
                {
1720
 
                        //Its a leaf.. call the callback
1721
 
                        if (root->totnode == 0)
1722
 
                        {
 
1649
                if (dist < data.radius) {
 
1650
                        /* Its a leaf.. call the callback */
 
1651
                        if (root->totnode == 0) {
1723
1652
                                data.hits++;
1724
1653
                                data.callback(data.userdata, root->index, dist);
1725
1654
                        }
1726
1655
                        else
1727
 
                                dfs_range_query( &data, root );
 
1656
                                dfs_range_query(&data, root);
1728
1657
                }
1729
1658
        }
1730
1659