Merged changes in the trunk up to revision 45619.
[blender.git] / source / blender / blenlib / intern / pbvh.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * ***** END GPL LICENSE BLOCK *****
19  */
20
21 /** \file blender/blenlib/intern/pbvh.c
22  *  \ingroup bli
23  */
24
25
26
27 #include "DNA_meshdata_types.h"
28
29 #include "MEM_guardedalloc.h"
30
31 #include "BLI_bitmap.h"
32 #include "BLI_math.h"
33 #include "BLI_utildefines.h"
34 #include "BLI_ghash.h"
35 #include "BLI_pbvh.h"
36
37 #include "BKE_DerivedMesh.h"
38 #include "BKE_mesh.h" /* for mesh_calc_normals */
39 #include "BKE_global.h" /* for mesh_calc_normals */
40 #include "BKE_paint.h"
41 #include "BKE_subsurf.h"
42
43 #include "GPU_buffers.h"
44
45 #define LEAF_LIMIT 10000
46
47 //#define PERFCNTRS
48
49 /* Axis-aligned bounding box */
50 typedef struct {
51         float bmin[3], bmax[3];
52 } BB;
53
54 /* Axis-aligned bounding box with centroid */
55 typedef struct {
56         float bmin[3], bmax[3], bcentroid[3];
57 } BBC;
58
59 struct PBVHNode {
60         /* Opaque handle for drawing code */
61         GPU_Buffers *draw_buffers;
62
63         /* Voxel bounds */
64         BB vb;
65         BB orig_vb;
66
67         /* For internal nodes, the offset of the children in the PBVH
68          * 'nodes' array. */
69         int children_offset;
70
71         /* Pointer into the PBVH prim_indices array and the number of
72          * primitives used by this leaf node.
73          *
74          * Used for leaf nodes in both mesh- and multires-based PBVHs.
75          */
76         int *prim_indices;
77         unsigned int totprim;
78
79         /* Array of indices into the mesh's MVert array. Contains the
80          * indices of all vertices used by faces that are within this
81          * node's bounding box.
82          *
83          * Note that a vertex might be used by a multiple faces, and
84          * these faces might be in different leaf nodes. Such a vertex
85          * will appear in the vert_indices array of each of those leaf
86          * nodes.
87          *
88          * In order to support cases where you want access to multiple
89          * nodes' vertices without duplication, the vert_indices array
90          * is ordered such that the first part of the array, up to
91          * index 'uniq_verts', contains "unique" vertex indices. These
92          * vertices might not be truly unique to this node, but if
93          * they appear in another node's vert_indices array, they will
94          * be above that node's 'uniq_verts' value.
95          *
96          * Used for leaf nodes in a mesh-based PBVH (not multires.)
97          */
98         int *vert_indices;
99         unsigned int uniq_verts, face_verts;
100
101         /* An array mapping face corners into the vert_indices
102          * array. The array is sized to match 'totprim', and each of
103          * the face's corners gets an index into the vert_indices
104          * array, in the same order as the corners in the original
105          * MFace. The fourth value should not be used if the original
106          * face is a triangle.
107          *
108          * Used for leaf nodes in a mesh-based PBVH (not multires.)
109          */
110         int (*face_vert_indices)[4];
111
112         /* Indicates whether this node is a leaf or not; also used for
113          * marking various updates that need to be applied. */
114         PBVHNodeFlags flag : 8;
115
116         /* Used for raycasting: how close bb is to the ray point. */
117         float tmin;
118
119         int proxy_count;
120         PBVHProxyNode* proxies;
121 };
122
123 struct PBVH {
124         PBVHType type;
125
126         PBVHNode *nodes;
127         int node_mem_count, totnode;
128
129         int *prim_indices;
130         int totprim;
131         int totvert;
132
133         int leaf_limit;
134
135         /* Mesh data */
136         MVert *verts;
137         MFace *faces;
138
139         /* Grid Data */
140         DMGridData **grids;
141         DMGridAdjacency *gridadj;
142         void **gridfaces;
143         const DMFlagMat *grid_flag_mats;
144         int totgrid;
145         int gridsize;
146         BLI_bitmap *grid_hidden;
147
148         /* Only used during BVH build and update,
149          * don't need to remain valid after */
150         BLI_bitmap vert_bitmap;
151
152 #ifdef PERFCNTRS
153         int perf_modified;
154 #endif
155
156         /* flag are verts/faces deformed */
157         int deformed;
158 };
159
160 #define STACK_FIXED_DEPTH       100
161
162 typedef struct PBVHStack {
163         PBVHNode *node;
164         int revisiting;
165 } PBVHStack;
166
167 typedef struct PBVHIter {
168         PBVH *bvh;
169         BLI_pbvh_SearchCallback scb;
170         void *search_data;
171
172         PBVHStack *stack;
173         int stacksize;
174
175         PBVHStack stackfixed[STACK_FIXED_DEPTH];
176         int stackspace;
177 } PBVHIter;
178
179 static void BB_reset(BB *bb)
180 {
181         bb->bmin[0] = bb->bmin[1] = bb->bmin[2] = FLT_MAX;
182         bb->bmax[0] = bb->bmax[1] = bb->bmax[2] = -FLT_MAX;
183 }
184
185 /* Expand the bounding box to include a new coordinate */
186 static void BB_expand(BB *bb, float co[3])
187 {
188         int i;
189         for (i = 0; i < 3; ++i) {
190                 bb->bmin[i] = MIN2(bb->bmin[i], co[i]);
191                 bb->bmax[i] = MAX2(bb->bmax[i], co[i]);
192         }
193 }
194
195 /* Expand the bounding box to include another bounding box */
196 static void BB_expand_with_bb(BB *bb, BB *bb2)
197 {
198         int i;
199         for (i = 0; i < 3; ++i) {
200                 bb->bmin[i] = MIN2(bb->bmin[i], bb2->bmin[i]);
201                 bb->bmax[i] = MAX2(bb->bmax[i], bb2->bmax[i]);
202         }
203 }
204
205 /* Return 0, 1, or 2 to indicate the widest axis of the bounding box */
206 static int BB_widest_axis(BB *bb)
207 {
208         float dim[3];
209         int i;
210
211         for (i = 0; i < 3; ++i)
212                 dim[i] = bb->bmax[i] - bb->bmin[i];
213
214         if (dim[0] > dim[1]) {
215                 if (dim[0] > dim[2])
216                         return 0;
217                 else
218                         return 2;
219         }
220         else {
221                 if (dim[1] > dim[2])
222                         return 1;
223                 else
224                         return 2;
225         }
226 }
227
228 static void BBC_update_centroid(BBC *bbc)
229 {
230         int i;
231         for (i = 0; i < 3; ++i)
232                 bbc->bcentroid[i] = (bbc->bmin[i] + bbc->bmax[i]) * 0.5f;
233 }
234
235 /* Not recursive */
236 static void update_node_vb(PBVH *bvh, PBVHNode *node)
237 {
238         BB vb;
239
240         BB_reset(&vb);
241         
242         if (node->flag & PBVH_Leaf) {
243                 PBVHVertexIter vd;
244
245                 BLI_pbvh_vertex_iter_begin(bvh, node, vd, PBVH_ITER_ALL) {
246                         BB_expand(&vb, vd.co);
247                 }
248                 BLI_pbvh_vertex_iter_end;
249         }
250         else {
251                 BB_expand_with_bb(&vb,
252                                   &bvh->nodes[node->children_offset].vb);
253                 BB_expand_with_bb(&vb,
254                                   &bvh->nodes[node->children_offset + 1].vb);
255         }
256
257         node->vb= vb;
258 }
259
260 //void BLI_pbvh_node_BB_reset(PBVHNode* node)
261 //{
262 //      BB_reset(&node->vb);
263 //}
264 //
265 //void BLI_pbvh_node_BB_expand(PBVHNode* node, float co[3])
266 //{
267 //      BB_expand(&node->vb, co);
268 //}
269
270 static int face_materials_match(const MFace *f1, const MFace *f2)
271 {
272         return ((f1->flag & ME_SMOOTH) == (f2->flag & ME_SMOOTH) &&
273                         (f1->mat_nr == f2->mat_nr));
274 }
275
276 static int grid_materials_match(const DMFlagMat *f1, const DMFlagMat *f2)
277 {
278         return ((f1->flag & ME_SMOOTH) == (f2->flag & ME_SMOOTH) &&
279                         (f1->mat_nr == f2->mat_nr));
280 }
281
282 /* Adapted from BLI_kdopbvh.c */
283 /* Returns the index of the first element on the right of the partition */
284 static int partition_indices(int *prim_indices, int lo, int hi, int axis,
285                                  float mid, BBC *prim_bbc)
286 {
287         int i=lo, j=hi;
288         for (;;) {
289                 for (; prim_bbc[prim_indices[i]].bcentroid[axis] < mid; i++);
290                 for (; mid < prim_bbc[prim_indices[j]].bcentroid[axis]; j--);
291                 
292                 if (!(i < j))
293                         return i;
294                 
295                 SWAP(int, prim_indices[i], prim_indices[j]);
296                 i++;
297         }
298 }
299
300 /* Returns the index of the first element on the right of the partition */
301 static int partition_indices_material(PBVH *bvh, int lo, int hi)
302 {
303         const MFace *faces = bvh->faces;
304         const DMFlagMat *flagmats = bvh->grid_flag_mats;
305         const int *indices = bvh->prim_indices;
306         const void *first;
307         int i=lo, j=hi;
308
309         if (bvh->faces)
310                 first = &faces[bvh->prim_indices[lo]];
311         else
312                 first = &flagmats[bvh->prim_indices[lo]];
313
314         for (;;) {
315                 if (bvh->faces) {
316                         for (; face_materials_match(first, &faces[indices[i]]); i++);
317                         for (; !face_materials_match(first, &faces[indices[j]]); j--);
318                 }
319                 else {
320                         for (; grid_materials_match(first, &flagmats[indices[i]]); i++);
321                         for (; !grid_materials_match(first, &flagmats[indices[j]]); j--);
322                 }
323                 
324                 if (!(i < j))
325                         return i;
326
327                 SWAP(int, bvh->prim_indices[i], bvh->prim_indices[j]);
328                 i++;
329         }
330 }
331
332 static void grow_nodes(PBVH *bvh, int totnode)
333 {
334         if (totnode > bvh->node_mem_count) {
335                 PBVHNode *prev = bvh->nodes;
336                 bvh->node_mem_count *= 1.33;
337                 if (bvh->node_mem_count < totnode)
338                         bvh->node_mem_count = totnode;
339                 bvh->nodes = MEM_callocN(sizeof(PBVHNode) * bvh->node_mem_count,
340                                          "bvh nodes");
341                 memcpy(bvh->nodes, prev, bvh->totnode * sizeof(PBVHNode));
342                 MEM_freeN(prev);
343         }
344
345         bvh->totnode = totnode;
346 }
347
348 /* Add a vertex to the map, with a positive value for unique vertices and
349  * a negative value for additional vertices */
350 static int map_insert_vert(PBVH *bvh, GHash *map,
351                                 unsigned int *face_verts,
352                                 unsigned int *uniq_verts, int vertex)
353 {
354         void *value, *key = SET_INT_IN_POINTER(vertex);
355
356         if (!BLI_ghash_haskey(map, key)) {
357                 if (BLI_BITMAP_GET(bvh->vert_bitmap, vertex)) {
358                         value = SET_INT_IN_POINTER(~(*face_verts));
359                         ++(*face_verts);
360                 }
361                 else {
362                         BLI_BITMAP_SET(bvh->vert_bitmap, vertex);
363                         value = SET_INT_IN_POINTER(*uniq_verts);
364                         ++(*uniq_verts);
365                 }
366                 
367                 BLI_ghash_insert(map, key, value);
368                 return GET_INT_FROM_POINTER(value);
369         }
370         else
371                 return GET_INT_FROM_POINTER(BLI_ghash_lookup(map, key));
372 }
373
374 /* Find vertices used by the faces in this node and update the draw buffers */
375 static void build_mesh_leaf_node(PBVH *bvh, PBVHNode *node)
376 {
377         GHashIterator *iter;
378         GHash *map;
379         int i, j, totface;
380
381         map = BLI_ghash_new(BLI_ghashutil_inthash, BLI_ghashutil_intcmp, "build_mesh_leaf_node gh");
382         
383         node->uniq_verts = node->face_verts = 0;
384         totface= node->totprim;
385
386         node->face_vert_indices = MEM_callocN(sizeof(int) * 4*totface,
387                                               "bvh node face vert indices");
388
389         for (i = 0; i < totface; ++i) {
390                 MFace *f = bvh->faces + node->prim_indices[i];
391                 int sides = f->v4 ? 4 : 3;
392
393                 for (j = 0; j < sides; ++j) {
394                         node->face_vert_indices[i][j]= 
395                                 map_insert_vert(bvh, map, &node->face_verts,
396                                                 &node->uniq_verts, (&f->v1)[j]);
397                 }
398         }
399
400         node->vert_indices = MEM_callocN(sizeof(int) *
401                                          (node->uniq_verts + node->face_verts),
402                                          "bvh node vert indices");
403
404         /* Build the vertex list, unique verts first */
405         for (iter = BLI_ghashIterator_new(map), i = 0;
406                 !BLI_ghashIterator_isDone(iter);
407                 BLI_ghashIterator_step(iter), ++i) {
408                 void *value = BLI_ghashIterator_getValue(iter);
409                 int ndx = GET_INT_FROM_POINTER(value);
410
411                 if (ndx < 0)
412                         ndx = -ndx + node->uniq_verts - 1;
413
414                 node->vert_indices[ndx] =
415                         GET_INT_FROM_POINTER(BLI_ghashIterator_getKey(iter));
416         }
417
418         BLI_ghashIterator_free(iter);
419
420         for (i = 0; i < totface; ++i) {
421                 MFace *f = bvh->faces + node->prim_indices[i];
422                 int sides = f->v4 ? 4 : 3;
423                 
424                 for (j = 0; j < sides; ++j) {
425                         if (node->face_vert_indices[i][j] < 0)
426                                 node->face_vert_indices[i][j]=
427                                         -node->face_vert_indices[i][j] +
428                                         node->uniq_verts - 1;
429                 }
430         }
431
432         if (!G.background) {
433                 node->draw_buffers =
434                         GPU_build_mesh_buffers(node->face_vert_indices,
435                                         bvh->faces, bvh->verts,
436                                         node->prim_indices,
437                                         node->totprim);
438         }
439
440         node->flag |= PBVH_UpdateDrawBuffers;
441
442         BLI_ghash_free(map, NULL, NULL);
443 }
444
445 static void build_grids_leaf_node(PBVH *bvh, PBVHNode *node)
446 {
447         if (!G.background) {
448                 node->draw_buffers =
449                         GPU_build_grid_buffers(node->prim_indices,
450                                 node->totprim, bvh->grid_hidden, bvh->gridsize);
451         }
452         node->flag |= PBVH_UpdateDrawBuffers;
453 }
454
455 static void update_vb(PBVH *bvh, PBVHNode *node, BBC *prim_bbc,
456                                           int offset, int count)
457 {
458         int i;
459         
460         BB_reset(&node->vb);
461         for (i = offset + count - 1; i >= offset; --i) {
462                 BB_expand_with_bb(&node->vb, (BB*)(&prim_bbc[bvh->prim_indices[i]]));
463         }
464         node->orig_vb = node->vb;
465 }
466
467 static void build_leaf(PBVH *bvh, int node_index, BBC *prim_bbc,
468                                            int offset, int count)
469 {
470         bvh->nodes[node_index].flag |= PBVH_Leaf;
471
472         bvh->nodes[node_index].prim_indices = bvh->prim_indices + offset;
473         bvh->nodes[node_index].totprim = count;
474
475         /* Still need vb for searches */
476         update_vb(bvh, &bvh->nodes[node_index], prim_bbc, offset, count);
477                 
478         if (bvh->faces)
479                 build_mesh_leaf_node(bvh, bvh->nodes + node_index);
480         else
481                 build_grids_leaf_node(bvh, bvh->nodes + node_index);
482 }
483
484 /* Return zero if all primitives in the node can be drawn with the
485  * same material (including flat/smooth shading), non-zerootherwise */
486 int leaf_needs_material_split(PBVH *bvh, int offset, int count)
487 {
488         int i, prim;
489
490         if (count <= 1)
491                 return 0;
492
493         if (bvh->faces) {
494                 const MFace *first = &bvh->faces[bvh->prim_indices[offset]];
495
496                 for (i = offset + count - 1; i > offset; --i) {
497                         prim = bvh->prim_indices[i];
498                         if (!face_materials_match(first, &bvh->faces[prim]))
499                                 return 1;
500                 }
501         }
502         else {
503                 const DMFlagMat *first = &bvh->grid_flag_mats[bvh->prim_indices[offset]];
504
505                 for (i = offset + count - 1; i > offset; --i) {
506                         prim = bvh->prim_indices[i];
507                         if (!grid_materials_match(first, &bvh->grid_flag_mats[prim]))
508                                 return 1;
509                 }
510         }
511
512         return 0;
513 }
514
515
516 /* Recursively build a node in the tree
517  *
518  * vb is the voxel box around all of the primitives contained in
519  * this node.
520  *
521  * cb is the bounding box around all the centroids of the primitives
522  * contained in this node
523  *
524  * offset and start indicate a range in the array of primitive indices
525  */
526
527 static void build_sub(PBVH *bvh, int node_index, BB *cb, BBC *prim_bbc,
528                                           int offset, int count)
529 {
530         int i, axis, end, below_leaf_limit;
531         BB cb_backing;
532
533         /* Decide whether this is a leaf or not */
534         below_leaf_limit = count <= bvh->leaf_limit;
535         if (below_leaf_limit) {
536                 if (!leaf_needs_material_split(bvh, offset, count)) {
537                         build_leaf(bvh, node_index, prim_bbc, offset, count);
538                         return;
539                 }
540         }
541
542         /* Add two child nodes */
543         bvh->nodes[node_index].children_offset = bvh->totnode;
544         grow_nodes(bvh, bvh->totnode + 2);
545
546         /* Update parent node bounding box */
547         update_vb(bvh, &bvh->nodes[node_index], prim_bbc, offset, count);
548
549         if (!below_leaf_limit) {
550                 /* Find axis with widest range of primitive centroids */
551                 if (!cb) {
552                         cb = &cb_backing;
553                         BB_reset(cb);
554                         for (i = offset + count - 1; i >= offset; --i)
555                                 BB_expand(cb, prim_bbc[bvh->prim_indices[i]].bcentroid);
556                 }
557                 axis = BB_widest_axis(cb);
558
559                 /* Partition primitives along that axis */
560                 end = partition_indices(bvh->prim_indices,
561                                                                 offset, offset + count - 1,
562                                                                 axis,
563                                                                 (cb->bmax[axis] + cb->bmin[axis]) * 0.5f,
564                                                                 prim_bbc);
565         }
566         else {
567                 /* Partition primitives by material */
568                 end = partition_indices_material(bvh, offset, offset + count - 1);
569         }
570
571         /* Build children */
572         build_sub(bvh, bvh->nodes[node_index].children_offset, NULL,
573                           prim_bbc, offset, end - offset);
574         build_sub(bvh, bvh->nodes[node_index].children_offset + 1, NULL,
575                           prim_bbc, end, offset + count - end);
576 }
577
578 static void pbvh_build(PBVH *bvh, BB *cb, BBC *prim_bbc, int totprim)
579 {
580         int i;
581
582         if (totprim != bvh->totprim) {
583                 bvh->totprim = totprim;
584                 if (bvh->nodes) MEM_freeN(bvh->nodes);
585                 if (bvh->prim_indices) MEM_freeN(bvh->prim_indices);
586                 bvh->prim_indices = MEM_callocN(sizeof(int) * totprim,
587                                                 "bvh prim indices");
588                 for (i = 0; i < totprim; ++i)
589                         bvh->prim_indices[i] = i;
590                 bvh->totnode = 0;
591                 if (bvh->node_mem_count < 100) {
592                         bvh->node_mem_count = 100;
593                         bvh->nodes = MEM_callocN(sizeof(PBVHNode) *
594                                                  bvh->node_mem_count,
595                                                  "bvh initial nodes");
596                 }
597         }
598
599         bvh->totnode = 1;
600         build_sub(bvh, 0, cb, prim_bbc, 0, totprim);
601 }
602
603 /* Do a full rebuild with on Mesh data structure */
604 void BLI_pbvh_build_mesh(PBVH *bvh, MFace *faces, MVert *verts, int totface, int totvert)
605 {
606         BBC *prim_bbc = NULL;
607         BB cb;
608         int i, j;
609
610         bvh->type = PBVH_FACES;
611         bvh->faces = faces;
612         bvh->verts = verts;
613         bvh->vert_bitmap = BLI_BITMAP_NEW(totvert, "bvh->vert_bitmap");
614         bvh->totvert = totvert;
615         bvh->leaf_limit = LEAF_LIMIT;
616
617         BB_reset(&cb);
618
619         /* For each face, store the AABB and the AABB centroid */
620         prim_bbc = MEM_mallocN(sizeof(BBC) * totface, "prim_bbc");
621
622         for (i = 0; i < totface; ++i) {
623                 MFace *f = faces + i;
624                 const int sides = f->v4 ? 4 : 3;
625                 BBC *bbc = prim_bbc + i;
626
627                 BB_reset((BB*)bbc);
628
629                 for (j = 0; j < sides; ++j)
630                         BB_expand((BB*)bbc, verts[(&f->v1)[j]].co);
631
632                 BBC_update_centroid(bbc);
633
634                 BB_expand(&cb, bbc->bcentroid);
635         }
636
637         if (totface)
638                 pbvh_build(bvh, &cb, prim_bbc, totface);
639
640         MEM_freeN(prim_bbc);
641         MEM_freeN(bvh->vert_bitmap);
642 }
643
644 /* Do a full rebuild with on Grids data structure */
645 void BLI_pbvh_build_grids(PBVH *bvh, DMGridData **grids, DMGridAdjacency *gridadj,
646         int totgrid, int gridsize, void **gridfaces, DMFlagMat *flagmats, BLI_bitmap *grid_hidden)
647 {
648         BBC *prim_bbc = NULL;
649         BB cb;
650         int i, j;
651
652         bvh->type = PBVH_GRIDS;
653         bvh->grids= grids;
654         bvh->gridadj= gridadj;
655         bvh->gridfaces= gridfaces;
656         bvh->grid_flag_mats= flagmats;
657         bvh->totgrid= totgrid;
658         bvh->gridsize= gridsize;
659         bvh->grid_hidden= grid_hidden;
660         bvh->leaf_limit = MAX2(LEAF_LIMIT/((gridsize-1)*(gridsize-1)), 1);
661
662         BB_reset(&cb);
663
664         /* For each grid, store the AABB and the AABB centroid */
665         prim_bbc = MEM_mallocN(sizeof(BBC) * totgrid, "prim_bbc");
666
667         for (i = 0; i < totgrid; ++i) {
668                 DMGridData *grid= grids[i];
669                 BBC *bbc = prim_bbc + i;
670
671                 BB_reset((BB*)bbc);
672
673                 for (j = 0; j < gridsize*gridsize; ++j)
674                         BB_expand((BB*)bbc, grid[j].co);
675
676                 BBC_update_centroid(bbc);
677
678                 BB_expand(&cb, bbc->bcentroid);
679         }
680
681         if (totgrid)
682                 pbvh_build(bvh, &cb, prim_bbc, totgrid);
683
684         MEM_freeN(prim_bbc);
685 }
686
687 PBVH *BLI_pbvh_new(void)
688 {
689         PBVH *bvh = MEM_callocN(sizeof(PBVH), "pbvh");
690
691         return bvh;
692 }
693
694 void BLI_pbvh_free(PBVH *bvh)
695 {
696         PBVHNode *node;
697         int i;
698
699         for (i = 0; i < bvh->totnode; ++i) {
700                 node= &bvh->nodes[i];
701
702                 if (node->flag & PBVH_Leaf) {
703                         if (node->draw_buffers)
704                                 GPU_free_buffers(node->draw_buffers);
705                         if (node->vert_indices)
706                                 MEM_freeN(node->vert_indices);
707                         if (node->face_vert_indices)
708                                 MEM_freeN(node->face_vert_indices);
709                 }
710         }
711
712         if (bvh->deformed) {
713                 if (bvh->verts) {
714                         /* if pbvh was deformed, new memory was allocated for verts/faces -- free it */
715
716                         MEM_freeN(bvh->verts);
717                         if (bvh->faces)
718                                 MEM_freeN(bvh->faces);
719                 }
720         }
721
722         if (bvh->nodes)
723                 MEM_freeN(bvh->nodes);
724
725         if (bvh->prim_indices)
726                 MEM_freeN(bvh->prim_indices);
727
728         MEM_freeN(bvh);
729 }
730
731 static void pbvh_iter_begin(PBVHIter *iter, PBVH *bvh, BLI_pbvh_SearchCallback scb, void *search_data)
732 {
733         iter->bvh= bvh;
734         iter->scb= scb;
735         iter->search_data= search_data;
736
737         iter->stack= iter->stackfixed;
738         iter->stackspace= STACK_FIXED_DEPTH;
739
740         iter->stack[0].node= bvh->nodes;
741         iter->stack[0].revisiting= 0;
742         iter->stacksize= 1;
743 }
744
745 static void pbvh_iter_end(PBVHIter *iter)
746 {
747         if (iter->stackspace > STACK_FIXED_DEPTH)
748                 MEM_freeN(iter->stack);
749 }
750
751 static void pbvh_stack_push(PBVHIter *iter, PBVHNode *node, int revisiting)
752 {
753         if (iter->stacksize == iter->stackspace) {
754                 PBVHStack *newstack;
755
756                 iter->stackspace *= 2;
757                 newstack= MEM_callocN(sizeof(PBVHStack)*iter->stackspace, "PBVHStack");
758                 memcpy(newstack, iter->stack, sizeof(PBVHStack)*iter->stacksize);
759
760                 if (iter->stackspace > STACK_FIXED_DEPTH)
761                         MEM_freeN(iter->stack);
762                 iter->stack= newstack;
763         }
764
765         iter->stack[iter->stacksize].node= node;
766         iter->stack[iter->stacksize].revisiting= revisiting;
767         iter->stacksize++;
768 }
769
770 static PBVHNode *pbvh_iter_next(PBVHIter *iter)
771 {
772         PBVHNode *node;
773         int revisiting;
774
775         /* purpose here is to traverse tree, visiting child nodes before their
776          * parents, this order is necessary for e.g. computing bounding boxes */
777
778         while (iter->stacksize) {
779                 /* pop node */
780                 iter->stacksize--;
781                 node= iter->stack[iter->stacksize].node;
782
783                 /* on a mesh with no faces this can happen
784                  * can remove this check if we know meshes have at least 1 face */
785                 if (node==NULL)
786                         return NULL;
787
788                 revisiting= iter->stack[iter->stacksize].revisiting;
789
790                 /* revisiting node already checked */
791                 if (revisiting)
792                         return node;
793
794                 if (iter->scb && !iter->scb(node, iter->search_data))
795                         continue; /* don't traverse, outside of search zone */
796
797                 if (node->flag & PBVH_Leaf) {
798                         /* immediately hit leaf node */
799                         return node;
800                 }
801                 else {
802                         /* come back later when children are done */
803                         pbvh_stack_push(iter, node, 1);
804
805                         /* push two child nodes on the stack */
806                         pbvh_stack_push(iter, iter->bvh->nodes+node->children_offset+1, 0);
807                         pbvh_stack_push(iter, iter->bvh->nodes+node->children_offset, 0);
808                 }
809         }
810
811         return NULL;
812 }
813
814 static PBVHNode *pbvh_iter_next_occluded(PBVHIter *iter)
815 {
816         PBVHNode *node;
817
818         while (iter->stacksize) {
819                 /* pop node */
820                 iter->stacksize--;
821                 node= iter->stack[iter->stacksize].node;
822
823                 /* on a mesh with no faces this can happen
824                  * can remove this check if we know meshes have at least 1 face */
825                 if (node==NULL) return NULL;
826
827                 if (iter->scb && !iter->scb(node, iter->search_data)) continue; /* don't traverse, outside of search zone */
828
829                 if (node->flag & PBVH_Leaf) {
830                         /* immediately hit leaf node */
831                         return node;
832                 }
833                 else {
834                         pbvh_stack_push(iter, iter->bvh->nodes+node->children_offset+1, 0);
835                         pbvh_stack_push(iter, iter->bvh->nodes+node->children_offset, 0);
836                 }
837         }
838
839         return NULL;
840 }
841
842 void BLI_pbvh_search_gather(PBVH *bvh,
843         BLI_pbvh_SearchCallback scb, void *search_data,
844         PBVHNode ***r_array, int *r_tot)
845 {
846         PBVHIter iter;
847         PBVHNode **array= NULL, **newarray, *node;
848         int tot= 0, space= 0;
849
850         pbvh_iter_begin(&iter, bvh, scb, search_data);
851
852         while ((node=pbvh_iter_next(&iter))) {
853                 if (node->flag & PBVH_Leaf) {
854                         if (tot == space) {
855                                 /* resize array if needed */
856                                 space= (tot == 0)? 32: space*2;
857                                 newarray= MEM_callocN(sizeof(PBVHNode)*space, "PBVHNodeSearch");
858
859                                 if (array) {
860                                         memcpy(newarray, array, sizeof(PBVHNode)*tot);
861                                         MEM_freeN(array);
862                                 }
863
864                                 array= newarray;
865                         }
866
867                         array[tot]= node;
868                         tot++;
869                 }
870         }
871
872         pbvh_iter_end(&iter);
873
874         if (tot == 0 && array) {
875                 MEM_freeN(array);
876                 array= NULL;
877         }
878
879         *r_array= array;
880         *r_tot= tot;
881 }
882
883 void BLI_pbvh_search_callback(PBVH *bvh,
884         BLI_pbvh_SearchCallback scb, void *search_data,
885         BLI_pbvh_HitCallback hcb, void *hit_data)
886 {
887         PBVHIter iter;
888         PBVHNode *node;
889
890         pbvh_iter_begin(&iter, bvh, scb, search_data);
891
892         while ((node=pbvh_iter_next(&iter)))
893                 if (node->flag & PBVH_Leaf)
894                         hcb(node, hit_data);
895
896         pbvh_iter_end(&iter);
897 }
898
899 typedef struct node_tree {
900         PBVHNode* data;
901
902         struct node_tree* left;
903         struct node_tree* right;
904 } node_tree;
905
906 static void node_tree_insert(node_tree* tree, node_tree* new_node)
907 {
908         if (new_node->data->tmin < tree->data->tmin) {
909                 if (tree->left) {
910                         node_tree_insert(tree->left, new_node);
911                 }
912                 else {
913                         tree->left = new_node;
914                 }
915         }
916         else {
917                 if (tree->right) {
918                         node_tree_insert(tree->right, new_node);
919                 }
920                 else {
921                         tree->right = new_node;
922                 }
923         }
924 }
925
926 static void traverse_tree(node_tree* tree, BLI_pbvh_HitOccludedCallback hcb, void* hit_data, float* tmin)
927 {
928         if (tree->left) traverse_tree(tree->left, hcb, hit_data, tmin);
929
930         hcb(tree->data, hit_data, tmin);
931
932         if (tree->right) traverse_tree(tree->right, hcb, hit_data, tmin);
933 }
934
935 static void free_tree(node_tree* tree)
936 {
937         if (tree->left) {
938                 free_tree(tree->left);
939                 tree->left = 0;
940         }
941
942         if (tree->right) {
943                 free_tree(tree->right);
944                 tree->right = 0;
945         }
946
947         free(tree);
948 }
949
950 float BLI_pbvh_node_get_tmin(PBVHNode* node)
951 {
952         return node->tmin;
953 }
954
955 static void BLI_pbvh_search_callback_occluded(PBVH *bvh,
956         BLI_pbvh_SearchCallback scb, void *search_data,
957         BLI_pbvh_HitOccludedCallback hcb, void *hit_data)
958 {
959         PBVHIter iter;
960         PBVHNode *node;
961         node_tree *tree = 0;
962
963         pbvh_iter_begin(&iter, bvh, scb, search_data);
964
965         while ((node=pbvh_iter_next_occluded(&iter))) {
966                 if (node->flag & PBVH_Leaf) {
967                         node_tree* new_node = malloc(sizeof(node_tree));
968
969                         new_node->data = node;
970
971                         new_node->left  = NULL;
972                         new_node->right = NULL;
973
974                         if (tree) {
975                                 node_tree_insert(tree, new_node);
976                         }
977                         else {
978                                 tree = new_node;
979                         }
980                 }
981         }
982
983         pbvh_iter_end(&iter);
984
985         if (tree) {
986                 float tmin = FLT_MAX;
987                 traverse_tree(tree, hcb, hit_data, &tmin);
988                 free_tree(tree);
989         }
990 }
991
992 static int update_search_cb(PBVHNode *node, void *data_v)
993 {
994         int flag= GET_INT_FROM_POINTER(data_v);
995
996         if (node->flag & PBVH_Leaf)
997                 return (node->flag & flag);
998         
999         return 1;
1000 }
1001
1002 static void pbvh_update_normals(PBVH *bvh, PBVHNode **nodes,
1003         int totnode, float (*face_nors)[3])
1004 {
1005         float (*vnor)[3];
1006         int n;
1007
1008         if (bvh->type != PBVH_FACES)
1009                 return;
1010
1011         /* could be per node to save some memory, but also means
1012          * we have to store for each vertex which node it is in */
1013         vnor= MEM_callocN(sizeof(float)*3*bvh->totvert, "bvh temp vnors");
1014
1015     /* subtle assumptions:
1016      * - We know that for all edited vertices, the nodes with faces
1017      *   adjacent to these vertices have been marked with PBVH_UpdateNormals.
1018      *   This is true because if the vertex is inside the brush radius, the
1019      *   bounding box of it's adjacent faces will be as well.
1020      * - However this is only true for the vertices that have actually been
1021      *   edited, not for all vertices in the nodes marked for update, so we
1022      *   can only update vertices marked with ME_VERT_PBVH_UPDATE.
1023      */
1024
1025         #pragma omp parallel for private(n) schedule(static)
1026         for (n = 0; n < totnode; n++) {
1027                 PBVHNode *node= nodes[n];
1028
1029                 if ((node->flag & PBVH_UpdateNormals)) {
1030                         int i, j, totface, *faces;
1031
1032                         faces= node->prim_indices;
1033                         totface= node->totprim;
1034
1035                         for (i = 0; i < totface; ++i) {
1036                                 MFace *f= bvh->faces + faces[i];
1037                                 float fn[3];
1038                                 unsigned int *fv = &f->v1;
1039                                 int sides= (f->v4)? 4: 3;
1040
1041                                 if (f->v4)
1042                                         normal_quad_v3(fn, bvh->verts[f->v1].co, bvh->verts[f->v2].co,
1043                                                                    bvh->verts[f->v3].co, bvh->verts[f->v4].co);
1044                                 else
1045                                         normal_tri_v3(fn, bvh->verts[f->v1].co, bvh->verts[f->v2].co,
1046                                                                   bvh->verts[f->v3].co);
1047
1048                                 for (j = 0; j < sides; ++j) {
1049                                         int v= fv[j];
1050
1051                                         if (bvh->verts[v].flag & ME_VERT_PBVH_UPDATE) {
1052                                                 /* this seems like it could be very slow but profile
1053                                                  * does not show this, so just leave it for now? */
1054                                                 #pragma omp atomic
1055                                                 vnor[v][0] += fn[0];
1056                                                 #pragma omp atomic
1057                                                 vnor[v][1] += fn[1];
1058                                                 #pragma omp atomic
1059                                                 vnor[v][2] += fn[2];
1060                                         }
1061                                 }
1062
1063                                 if (face_nors)
1064                                         copy_v3_v3(face_nors[faces[i]], fn);
1065                         }
1066                 }
1067         }
1068
1069         #pragma omp parallel for private(n) schedule(static)
1070         for (n = 0; n < totnode; n++) {
1071                 PBVHNode *node= nodes[n];
1072
1073                 if (node->flag & PBVH_UpdateNormals) {
1074                         int i, *verts, totvert;
1075
1076                         verts= node->vert_indices;
1077                         totvert= node->uniq_verts;
1078
1079                         for (i = 0; i < totvert; ++i) {
1080                                 const int v = verts[i];
1081                                 MVert *mvert= &bvh->verts[v];
1082
1083                                 if (mvert->flag & ME_VERT_PBVH_UPDATE) {
1084                                         float no[3];
1085
1086                                         copy_v3_v3(no, vnor[v]);
1087                                         normalize_v3(no);
1088                                         
1089                                         mvert->no[0] = (short)(no[0]*32767.0f);
1090                                         mvert->no[1] = (short)(no[1]*32767.0f);
1091                                         mvert->no[2] = (short)(no[2]*32767.0f);
1092                                         
1093                                         mvert->flag &= ~ME_VERT_PBVH_UPDATE;
1094                                 }
1095                         }
1096
1097                         node->flag &= ~PBVH_UpdateNormals;
1098                 }
1099         }
1100
1101         MEM_freeN(vnor);
1102 }
1103
1104 static void pbvh_update_BB_redraw(PBVH *bvh, PBVHNode **nodes,
1105         int totnode, int flag)
1106 {
1107         int n;
1108
1109         /* update BB, redraw flag */
1110         #pragma omp parallel for private(n) schedule(static)
1111         for (n = 0; n < totnode; n++) {
1112                 PBVHNode *node= nodes[n];
1113
1114                 if ((flag & PBVH_UpdateBB) && (node->flag & PBVH_UpdateBB))
1115                         /* don't clear flag yet, leave it for flushing later */
1116                         update_node_vb(bvh, node);
1117
1118                 if ((flag & PBVH_UpdateOriginalBB) && (node->flag & PBVH_UpdateOriginalBB))
1119                         node->orig_vb= node->vb;
1120
1121                 if ((flag & PBVH_UpdateRedraw) && (node->flag & PBVH_UpdateRedraw))
1122                         node->flag &= ~PBVH_UpdateRedraw;
1123         }
1124 }
1125
1126 static void pbvh_update_draw_buffers(PBVH *bvh, PBVHNode **nodes, int totnode)
1127 {
1128         PBVHNode *node;
1129         int n;
1130
1131         /* can't be done in parallel with OpenGL */
1132         for (n = 0; n < totnode; n++) {
1133                 node= nodes[n];
1134
1135                 if (node->flag & PBVH_RebuildDrawBuffers) {
1136                         GPU_free_buffers(node->draw_buffers);
1137                         if (bvh->grids) {
1138                                 node->draw_buffers =
1139                                         GPU_build_grid_buffers(node->prim_indices,
1140                                                    node->totprim, bvh->grid_hidden, bvh->gridsize);
1141                         }
1142                         else {
1143                                 node->draw_buffers =
1144                                         GPU_build_mesh_buffers(node->face_vert_indices,
1145                                                    bvh->faces, bvh->verts,
1146                                                    node->prim_indices,
1147                                                    node->totprim);
1148                         }
1149  
1150                         node->flag &= ~PBVH_RebuildDrawBuffers;
1151                 }
1152
1153                 if (node->flag & PBVH_UpdateDrawBuffers) {
1154                         switch(bvh->type) {
1155                         case PBVH_GRIDS:
1156                                 GPU_update_grid_buffers(node->draw_buffers,
1157                                                    bvh->grids,
1158                                                    bvh->grid_flag_mats,
1159                                                    node->prim_indices,
1160                                                    node->totprim,
1161                                                    bvh->gridsize);
1162                                 break;
1163                         case PBVH_FACES:
1164                                 GPU_update_mesh_buffers(node->draw_buffers,
1165                                                    bvh->verts,
1166                                                    node->vert_indices,
1167                                                    node->uniq_verts +
1168                                                    node->face_verts);
1169                                 break;
1170                         }
1171
1172                         node->flag &= ~PBVH_UpdateDrawBuffers;
1173                 }
1174         }
1175 }
1176
1177 static int pbvh_flush_bb(PBVH *bvh, PBVHNode *node, int flag)
1178 {
1179         int update= 0;
1180
1181         /* difficult to multithread well, we just do single threaded recursive */
1182         if (node->flag & PBVH_Leaf) {
1183                 if (flag & PBVH_UpdateBB) {
1184                         update |= (node->flag & PBVH_UpdateBB);
1185                         node->flag &= ~PBVH_UpdateBB;
1186                 }
1187
1188                 if (flag & PBVH_UpdateOriginalBB) {
1189                         update |= (node->flag & PBVH_UpdateOriginalBB);
1190                         node->flag &= ~PBVH_UpdateOriginalBB;
1191                 }
1192
1193                 return update;
1194         }
1195         else {
1196                 update |= pbvh_flush_bb(bvh, bvh->nodes + node->children_offset, flag);
1197                 update |= pbvh_flush_bb(bvh, bvh->nodes + node->children_offset + 1, flag);
1198
1199                 if (update & PBVH_UpdateBB)
1200                         update_node_vb(bvh, node);
1201                 if (update & PBVH_UpdateOriginalBB)
1202                         node->orig_vb= node->vb;
1203         }
1204
1205         return update;
1206 }
1207
1208 void BLI_pbvh_update(PBVH *bvh, int flag, float (*face_nors)[3])
1209 {
1210         PBVHNode **nodes;
1211         int totnode;
1212
1213         if (!bvh->nodes)
1214                 return;
1215
1216         BLI_pbvh_search_gather(bvh, update_search_cb, SET_INT_IN_POINTER(flag),
1217                 &nodes, &totnode);
1218
1219         if (flag & PBVH_UpdateNormals)
1220                 pbvh_update_normals(bvh, nodes, totnode, face_nors);
1221
1222         if (flag & (PBVH_UpdateBB|PBVH_UpdateOriginalBB|PBVH_UpdateRedraw))
1223                 pbvh_update_BB_redraw(bvh, nodes, totnode, flag);
1224
1225         if (flag & (PBVH_UpdateBB|PBVH_UpdateOriginalBB))
1226                 pbvh_flush_bb(bvh, bvh->nodes, flag);
1227
1228         if (nodes) MEM_freeN(nodes);
1229 }
1230
1231 void BLI_pbvh_redraw_BB(PBVH *bvh, float bb_min[3], float bb_max[3])
1232 {
1233         PBVHIter iter;
1234         PBVHNode *node;
1235         BB bb;
1236
1237         BB_reset(&bb);
1238
1239         pbvh_iter_begin(&iter, bvh, NULL, NULL);
1240
1241         while ((node=pbvh_iter_next(&iter)))
1242                 if (node->flag & PBVH_UpdateRedraw)
1243                         BB_expand_with_bb(&bb, &node->vb);
1244
1245         pbvh_iter_end(&iter);
1246
1247         copy_v3_v3(bb_min, bb.bmin);
1248         copy_v3_v3(bb_max, bb.bmax);
1249 }
1250
1251 void BLI_pbvh_get_grid_updates(PBVH *bvh, int clear, void ***gridfaces, int *totface)
1252 {
1253         PBVHIter iter;
1254         PBVHNode *node;
1255         GHashIterator *hiter;
1256         GHash *map;
1257         void *face, **faces;
1258         unsigned i;
1259         int tot;
1260
1261         map = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp, "pbvh_get_grid_updates gh");
1262
1263         pbvh_iter_begin(&iter, bvh, NULL, NULL);
1264
1265         while ((node=pbvh_iter_next(&iter))) {
1266                 if (node->flag & PBVH_UpdateNormals) {
1267                         for (i = 0; i < node->totprim; ++i) {
1268                                 face= bvh->gridfaces[node->prim_indices[i]];
1269                                 if (!BLI_ghash_lookup(map, face))
1270                                         BLI_ghash_insert(map, face, face);
1271                         }
1272
1273                         if (clear)
1274                                 node->flag &= ~PBVH_UpdateNormals;
1275                 }
1276         }
1277
1278         pbvh_iter_end(&iter);
1279         
1280         tot= BLI_ghash_size(map);
1281         if (tot == 0) {
1282                 *totface= 0;
1283                 *gridfaces= NULL;
1284                 BLI_ghash_free(map, NULL, NULL);
1285                 return;
1286         }
1287
1288         faces= MEM_callocN(sizeof(void*)*tot, "PBVH Grid Faces");
1289
1290         for (hiter = BLI_ghashIterator_new(map), i = 0;
1291                 !BLI_ghashIterator_isDone(hiter);
1292                 BLI_ghashIterator_step(hiter), ++i)
1293                 faces[i]= BLI_ghashIterator_getKey(hiter);
1294
1295         BLI_ghashIterator_free(hiter);
1296
1297         BLI_ghash_free(map, NULL, NULL);
1298
1299         *totface= tot;
1300         *gridfaces= faces;
1301 }
1302
1303 /***************************** PBVH Access ***********************************/
1304
1305 PBVHType BLI_pbvh_type(const PBVH *bvh)
1306 {
1307         return bvh->type;
1308 }
1309
1310 BLI_bitmap *BLI_pbvh_grid_hidden(const PBVH *bvh)
1311 {
1312         BLI_assert(bvh->type == PBVH_GRIDS);
1313         return bvh->grid_hidden;
1314 }
1315
1316 /***************************** Node Access ***********************************/
1317
1318 void BLI_pbvh_node_mark_update(PBVHNode *node)
1319 {
1320         node->flag |= PBVH_UpdateNormals|PBVH_UpdateBB|PBVH_UpdateOriginalBB|PBVH_UpdateDrawBuffers|PBVH_UpdateRedraw;
1321 }
1322
1323 void BLI_pbvh_node_mark_rebuild_draw(PBVHNode *node)
1324 {
1325         node->flag |= PBVH_RebuildDrawBuffers|PBVH_UpdateDrawBuffers|PBVH_UpdateRedraw;
1326 }
1327
1328 void BLI_pbvh_node_fully_hidden_set(PBVHNode *node, int fully_hidden)
1329 {
1330         BLI_assert(node->flag & PBVH_Leaf);
1331         
1332         if (fully_hidden)
1333                 node->flag |= PBVH_FullyHidden;
1334         else
1335                 node->flag &= ~PBVH_FullyHidden;
1336 }
1337
1338 void BLI_pbvh_node_get_verts(PBVH *bvh, PBVHNode *node, int **vert_indices, MVert **verts)
1339 {
1340         if (vert_indices) *vert_indices= node->vert_indices;
1341         if (verts) *verts= bvh->verts;
1342 }
1343
1344 void BLI_pbvh_node_num_verts(PBVH *bvh, PBVHNode *node, int *uniquevert, int *totvert)
1345 {
1346         int tot;
1347         
1348         switch(bvh->type) {
1349         case PBVH_GRIDS:
1350                 tot= node->totprim*bvh->gridsize*bvh->gridsize;
1351                 if (totvert) *totvert= tot;
1352                 if (uniquevert) *uniquevert= tot;
1353                 break;
1354         case PBVH_FACES:
1355                 if (totvert) *totvert= node->uniq_verts + node->face_verts;
1356                 if (uniquevert) *uniquevert= node->uniq_verts;
1357                 break;
1358         }
1359 }
1360
1361 void BLI_pbvh_node_get_grids(PBVH *bvh, PBVHNode *node, int **grid_indices, int *totgrid, int *maxgrid, int *gridsize, DMGridData ***griddata, DMGridAdjacency **gridadj)
1362 {
1363         switch(bvh->type) {
1364         case PBVH_GRIDS:
1365                 if (grid_indices) *grid_indices= node->prim_indices;
1366                 if (totgrid) *totgrid= node->totprim;
1367                 if (maxgrid) *maxgrid= bvh->totgrid;
1368                 if (gridsize) *gridsize= bvh->gridsize;
1369                 if (griddata) *griddata= bvh->grids;
1370                 if (gridadj) *gridadj= bvh->gridadj;
1371                 break;
1372         case PBVH_FACES:
1373                 if (grid_indices) *grid_indices= NULL;
1374                 if (totgrid) *totgrid= 0;
1375                 if (maxgrid) *maxgrid= 0;
1376                 if (gridsize) *gridsize= 0;
1377                 if (griddata) *griddata= NULL;
1378                 if (gridadj) *gridadj= NULL;
1379                 break;
1380         }
1381 }
1382
1383 void BLI_pbvh_node_get_BB(PBVHNode *node, float bb_min[3], float bb_max[3])
1384 {
1385         copy_v3_v3(bb_min, node->vb.bmin);
1386         copy_v3_v3(bb_max, node->vb.bmax);
1387 }
1388
1389 void BLI_pbvh_node_get_original_BB(PBVHNode *node, float bb_min[3], float bb_max[3])
1390 {
1391         copy_v3_v3(bb_min, node->orig_vb.bmin);
1392         copy_v3_v3(bb_max, node->orig_vb.bmax);
1393 }
1394
1395 void BLI_pbvh_node_get_proxies(PBVHNode* node, PBVHProxyNode** proxies, int* proxy_count)
1396 {
1397         if (node->proxy_count > 0) {
1398                 if (proxies) *proxies = node->proxies;
1399                 if (proxy_count) *proxy_count = node->proxy_count;
1400         }
1401         else {
1402                 if (proxies) *proxies = 0;
1403                 if (proxy_count) *proxy_count = 0;
1404         }
1405 }
1406
1407 /********************************* Raycast ***********************************/
1408
1409 typedef struct {
1410         /* Ray */
1411         float start[3];
1412         int sign[3];
1413         float inv_dir[3];
1414         int original;
1415 } RaycastData;
1416
1417 /* Adapted from here: http://www.gamedev.net/community/forums/topic.asp?topic_id=459973 */
1418 static int ray_aabb_intersect(PBVHNode *node, void *data_v)
1419 {
1420         RaycastData *ray = data_v;
1421         float bbox[2][3];
1422         float tmin, tmax, tymin, tymax, tzmin, tzmax;
1423
1424         if (ray->original)
1425                 BLI_pbvh_node_get_original_BB(node, bbox[0], bbox[1]);
1426         else
1427                 BLI_pbvh_node_get_BB(node, bbox[0], bbox[1]);
1428
1429         tmin = (bbox[ray->sign[0]][0] - ray->start[0]) * ray->inv_dir[0];
1430         tmax = (bbox[1-ray->sign[0]][0] - ray->start[0]) * ray->inv_dir[0];
1431
1432         tymin = (bbox[ray->sign[1]][1] - ray->start[1]) * ray->inv_dir[1];
1433         tymax = (bbox[1-ray->sign[1]][1] - ray->start[1]) * ray->inv_dir[1];
1434
1435         if ((tmin > tymax) || (tymin > tmax))
1436                 return 0;
1437
1438         if (tymin > tmin)
1439                 tmin = tymin;
1440
1441         if (tymax < tmax)
1442                 tmax = tymax;
1443
1444         tzmin = (bbox[ray->sign[2]][2] - ray->start[2]) * ray->inv_dir[2];
1445         tzmax = (bbox[1-ray->sign[2]][2] - ray->start[2]) * ray->inv_dir[2];
1446
1447         if ((tmin > tzmax) || (tzmin > tmax))
1448                 return 0;
1449
1450         if (tzmin > tmin)
1451                 tmin = tzmin;
1452
1453         // XXX jwilkins: tmax does not need to be updated since we don't use it
1454         // keeping this here for future reference
1455         //if(tzmax < tmax) tmax = tzmax; 
1456
1457         node->tmin = tmin;
1458
1459         return 1;
1460 }
1461
1462 void BLI_pbvh_raycast(PBVH *bvh, BLI_pbvh_HitOccludedCallback cb, void *data,
1463                           float ray_start[3], float ray_normal[3], int original)
1464 {
1465         RaycastData rcd;
1466
1467         copy_v3_v3(rcd.start, ray_start);
1468         rcd.inv_dir[0] = 1.0f / ray_normal[0];
1469         rcd.inv_dir[1] = 1.0f / ray_normal[1];
1470         rcd.inv_dir[2] = 1.0f / ray_normal[2];
1471         rcd.sign[0] = rcd.inv_dir[0] < 0;
1472         rcd.sign[1] = rcd.inv_dir[1] < 0;
1473         rcd.sign[2] = rcd.inv_dir[2] < 0;
1474         rcd.original = original;
1475
1476         BLI_pbvh_search_callback_occluded(bvh, ray_aabb_intersect, &rcd, cb, data);
1477 }
1478
1479 static int ray_face_intersection(float ray_start[3], float ray_normal[3],
1480                                  float *t0, float *t1, float *t2, float *t3,
1481                                  float *fdist)
1482 {
1483         float dist;
1484
1485         if ((isect_ray_tri_epsilon_v3(ray_start, ray_normal, t0, t1, t2, &dist, NULL, 0.1f) && dist < *fdist) ||
1486            (t3 && isect_ray_tri_epsilon_v3(ray_start, ray_normal, t0, t2, t3, &dist, NULL, 0.1f) && dist < *fdist))
1487         {
1488                 *fdist = dist;
1489                 return 1;
1490         }
1491         else {
1492                 return 0;
1493         }
1494 }
1495
1496 int BLI_pbvh_node_raycast(PBVH *bvh, PBVHNode *node, float (*origco)[3],
1497         float ray_start[3], float ray_normal[3], float *dist)
1498 {
1499         MVert *vert;
1500         BLI_bitmap gh;
1501         int *faces, totface, gridsize, totgrid;
1502         int i, x, y, hit= 0;
1503
1504         if (node->flag & PBVH_FullyHidden)
1505                 return 0;
1506
1507         switch(bvh->type) {
1508         case PBVH_FACES:
1509                 vert = bvh->verts;
1510                 faces= node->prim_indices;
1511                 totface= node->totprim;
1512
1513                 for (i = 0; i < totface; ++i) {
1514                         const MFace *f = bvh->faces + faces[i];
1515                         int *face_verts = node->face_vert_indices[i];
1516
1517                         if (paint_is_face_hidden(f, vert))
1518                                 continue;
1519
1520                         if (origco) {
1521                                 /* intersect with backuped original coordinates */
1522                                 hit |= ray_face_intersection(ray_start, ray_normal,
1523                                                          origco[face_verts[0]],
1524                                                          origco[face_verts[1]],
1525                                                          origco[face_verts[2]],
1526                                                          f->v4? origco[face_verts[3]]: NULL,
1527                                                          dist);
1528                         }
1529                         else {
1530                                 /* intersect with current coordinates */
1531                                 hit |= ray_face_intersection(ray_start, ray_normal,
1532                                                          vert[f->v1].co,
1533                                                          vert[f->v2].co,
1534                                                          vert[f->v3].co,
1535                                                          f->v4 ? vert[f->v4].co : NULL,
1536                                                          dist);
1537                         }
1538                 }
1539                 break;
1540         case PBVH_GRIDS:
1541                 totgrid= node->totprim;
1542                 gridsize= bvh->gridsize;
1543
1544                 for (i = 0; i < totgrid; ++i) {
1545                         DMGridData *grid= bvh->grids[node->prim_indices[i]];
1546                         if (!grid)
1547                                 continue;
1548
1549                         gh= bvh->grid_hidden[node->prim_indices[i]];
1550
1551                         for (y = 0; y < gridsize-1; ++y) {
1552                                 for (x = 0; x < gridsize-1; ++x) {
1553                                         /* check if grid face is hidden */
1554                                         if (gh) {
1555                                                 if (paint_is_grid_face_hidden(gh, gridsize, x, y))
1556                                                         continue;
1557                                         }
1558
1559                                         if (origco) {
1560                                                 hit |= ray_face_intersection(ray_start, ray_normal,
1561                                                                          origco[y*gridsize + x],
1562                                                                          origco[y*gridsize + x+1],
1563                                                                          origco[(y+1)*gridsize + x+1],
1564                                                                          origco[(y+1)*gridsize + x],
1565                                                                          dist);
1566                                         }
1567                                         else {
1568                                                 hit |= ray_face_intersection(ray_start, ray_normal,
1569                                                                          grid[y*gridsize + x].co,
1570                                                                          grid[y*gridsize + x+1].co,
1571                                                                          grid[(y+1)*gridsize + x+1].co,
1572                                                                          grid[(y+1)*gridsize + x].co,
1573                                                                          dist);
1574                                         }
1575                                 }
1576                         }
1577
1578                         if (origco)
1579                                 origco += gridsize*gridsize;
1580                 }
1581                 break;
1582         }
1583
1584         return hit;
1585 }
1586
1587 //#include <GL/glew.h>
1588
1589 void BLI_pbvh_node_draw(PBVHNode *node, void *setMaterial)
1590 {
1591 #if 0
1592         /* XXX: Just some quick code to show leaf nodes in different colors */
1593         float col[3]; int i;
1594
1595         if (0) { //is_partial) {
1596                 col[0] = (rand() / (float)RAND_MAX); col[1] = col[2] = 0.6;
1597         }
1598         else {
1599                 srand((long long)node);
1600                 for (i = 0; i < 3; ++i)
1601                         col[i] = (rand() / (float)RAND_MAX) * 0.3 + 0.7;
1602         }
1603         glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, col);
1604
1605         glColor3f(1, 0, 0);
1606 #endif
1607
1608         if (!(node->flag & PBVH_FullyHidden))
1609                 GPU_draw_buffers(node->draw_buffers, setMaterial);
1610 }
1611
1612 typedef enum {
1613         ISECT_INSIDE,
1614         ISECT_OUTSIDE,
1615         ISECT_INTERSECT
1616 } PlaneAABBIsect;
1617
1618 /* Adapted from:
1619  * http://www.gamedev.net/community/forums/topic.asp?topic_id=512123
1620  * Returns true if the AABB is at least partially within the frustum
1621  * (ok, not a real frustum), false otherwise.
1622  */
1623 static PlaneAABBIsect test_planes_aabb(const float bb_min[3],
1624                                                                            const float bb_max[3],
1625                                                                            const float (*planes)[4])
1626 {
1627         float vmin[3], vmax[3];
1628         PlaneAABBIsect ret = ISECT_INSIDE;
1629         int i, axis;
1630         
1631         for (i = 0; i < 4; ++i) { 
1632                 for (axis = 0; axis < 3; ++axis) {
1633                         if (planes[i][axis] > 0) { 
1634                                 vmin[axis] = bb_min[axis];
1635                                 vmax[axis] = bb_max[axis];
1636                         }
1637                         else {
1638                                 vmin[axis] = bb_max[axis];
1639                                 vmax[axis] = bb_min[axis];
1640                         }
1641                 }
1642                 
1643                 if (dot_v3v3(planes[i], vmin) + planes[i][3] > 0)
1644                         return ISECT_OUTSIDE;
1645                 else if (dot_v3v3(planes[i], vmax) + planes[i][3] >= 0)
1646                         ret = ISECT_INTERSECT;
1647         } 
1648
1649         return ret;
1650 }
1651
1652 int BLI_pbvh_node_planes_contain_AABB(PBVHNode *node, void *data)
1653 {
1654         float bb_min[3], bb_max[3];
1655         
1656         BLI_pbvh_node_get_BB(node, bb_min, bb_max);
1657         return test_planes_aabb(bb_min, bb_max, data) != ISECT_OUTSIDE;
1658 }
1659
1660 int BLI_pbvh_node_planes_exclude_AABB(PBVHNode *node, void *data)
1661 {
1662         float bb_min[3], bb_max[3];
1663         
1664         BLI_pbvh_node_get_BB(node, bb_min, bb_max);
1665         return test_planes_aabb(bb_min, bb_max, data) != ISECT_INSIDE;
1666 }
1667
1668 void BLI_pbvh_draw(PBVH *bvh, float (*planes)[4], float (*face_nors)[3],
1669                                    DMSetMaterial setMaterial)
1670 {
1671         PBVHNode **nodes;
1672         int totnode;
1673
1674         BLI_pbvh_search_gather(bvh, update_search_cb, SET_INT_IN_POINTER(PBVH_UpdateNormals|PBVH_UpdateDrawBuffers),
1675                 &nodes, &totnode);
1676
1677         pbvh_update_normals(bvh, nodes, totnode, face_nors);
1678         pbvh_update_draw_buffers(bvh, nodes, totnode);
1679
1680         if (nodes) MEM_freeN(nodes);
1681
1682         if (planes) {
1683                 BLI_pbvh_search_callback(bvh, BLI_pbvh_node_planes_contain_AABB,
1684                                 planes, BLI_pbvh_node_draw, setMaterial);
1685         }
1686         else {
1687                 BLI_pbvh_search_callback(bvh, NULL, NULL, BLI_pbvh_node_draw, setMaterial);
1688         }
1689 }
1690
1691 void BLI_pbvh_grids_update(PBVH *bvh, DMGridData **grids, DMGridAdjacency *gridadj, void **gridfaces)
1692 {
1693         bvh->grids= grids;
1694         bvh->gridadj= gridadj;
1695         bvh->gridfaces= gridfaces;
1696 }
1697
1698 float (*BLI_pbvh_get_vertCos(PBVH *pbvh))[3]
1699 {
1700         int a;
1701         float (*vertCos)[3]= NULL;
1702
1703         if (pbvh->verts) {
1704                 float *co;
1705                 MVert *mvert= pbvh->verts;
1706
1707                 vertCos= MEM_callocN(3*pbvh->totvert*sizeof(float), "BLI_pbvh_get_vertCoords");
1708                 co= (float*)vertCos;
1709
1710                 for (a= 0; a<pbvh->totvert; a++, mvert++, co+= 3) {
1711                         copy_v3_v3(co, mvert->co);
1712                 }
1713         }
1714
1715         return vertCos;
1716 }
1717
1718 void BLI_pbvh_apply_vertCos(PBVH *pbvh, float (*vertCos)[3])
1719 {
1720         int a;
1721
1722         if (!pbvh->deformed) {
1723                 if (pbvh->verts) {
1724                         /* if pbvh is not already deformed, verts/faces points to the */
1725                         /* original data and applying new coords to this arrays would lead to */
1726                         /* unneeded deformation -- duplicate verts/faces to avoid this */
1727
1728                         pbvh->verts= MEM_dupallocN(pbvh->verts);
1729                         pbvh->faces= MEM_dupallocN(pbvh->faces);
1730
1731                         pbvh->deformed= 1;
1732                 }
1733         }
1734
1735         if (pbvh->verts) {
1736                 MVert *mvert= pbvh->verts;
1737                 /* copy new verts coords */
1738                 for (a= 0; a < pbvh->totvert; ++a, ++mvert) {
1739                         copy_v3_v3(mvert->co, vertCos[a]);
1740                         mvert->flag |= ME_VERT_PBVH_UPDATE;
1741                 }
1742
1743                 /* coordinates are new -- normals should also be updated */
1744                 mesh_calc_normals_tessface(pbvh->verts, pbvh->totvert, pbvh->faces, pbvh->totprim, NULL);
1745
1746                 for (a= 0; a < pbvh->totnode; ++a)
1747                         BLI_pbvh_node_mark_update(&pbvh->nodes[a]);
1748
1749                 BLI_pbvh_update(pbvh, PBVH_UpdateBB, NULL);
1750                 BLI_pbvh_update(pbvh, PBVH_UpdateOriginalBB, NULL);
1751
1752         }
1753 }
1754
1755 int BLI_pbvh_isDeformed(PBVH *pbvh)
1756 {
1757         return pbvh->deformed;
1758 }
1759 /* Proxies */
1760
1761 PBVHProxyNode* BLI_pbvh_node_add_proxy(PBVH* bvh, PBVHNode* node)
1762 {
1763         int index, totverts;
1764
1765         #pragma omp critical
1766         {
1767
1768                 index = node->proxy_count;
1769
1770                 node->proxy_count++;
1771
1772                 if (node->proxies)
1773                         node->proxies= MEM_reallocN(node->proxies, node->proxy_count*sizeof(PBVHProxyNode));
1774                 else
1775                         node->proxies= MEM_mallocN(sizeof(PBVHProxyNode), "PBVHNodeProxy");
1776
1777                 if (bvh->grids)
1778                         totverts = node->totprim*bvh->gridsize*bvh->gridsize;
1779                 else
1780                         totverts = node->uniq_verts;
1781
1782                 node->proxies[index].co= MEM_callocN(sizeof(float[3])*totverts, "PBVHNodeProxy.co");
1783         }
1784
1785         return node->proxies + index;
1786 }
1787
1788 void BLI_pbvh_node_free_proxies(PBVHNode* node)
1789 {
1790         #pragma omp critical
1791         {
1792                 int p;
1793
1794                 for (p= 0; p < node->proxy_count; p++) {
1795                         MEM_freeN(node->proxies[p].co);
1796                         node->proxies[p].co= 0;
1797                 }
1798
1799                 MEM_freeN(node->proxies);
1800                 node->proxies = 0;
1801
1802                 node->proxy_count= 0;
1803         }
1804 }
1805
1806 void BLI_pbvh_gather_proxies(PBVH* pbvh, PBVHNode*** r_array,  int* r_tot)
1807 {
1808         PBVHNode **array= NULL, **newarray, *node;
1809         int tot= 0, space= 0;
1810         int n;
1811
1812         for (n= 0; n < pbvh->totnode; n++) {
1813                 node = pbvh->nodes + n;
1814
1815                 if (node->proxy_count > 0) {
1816                         if (tot == space) {
1817                                 /* resize array if needed */
1818                                 space= (tot == 0)? 32: space*2;
1819                                 newarray= MEM_callocN(sizeof(PBVHNode)*space, "BLI_pbvh_gather_proxies");
1820
1821                                 if (array) {
1822                                         memcpy(newarray, array, sizeof(PBVHNode)*tot);
1823                                         MEM_freeN(array);
1824                                 }
1825
1826                                 array= newarray;
1827                         }
1828
1829                         array[tot]= node;
1830                         tot++;
1831                 }
1832         }
1833
1834         if (tot == 0 && array) {
1835                 MEM_freeN(array);
1836                 array= NULL;
1837         }
1838
1839         *r_array= array;
1840         *r_tot= tot;
1841 }
1842
1843 void pbvh_vertex_iter_init(PBVH *bvh, PBVHNode *node,
1844                                                    PBVHVertexIter *vi, int mode)
1845 {
1846         struct DMGridData **grids;
1847         struct MVert *verts;
1848         int *grid_indices, *vert_indices;
1849         int totgrid, gridsize, uniq_verts, totvert;
1850         
1851         vi->grid= 0;
1852         vi->no= 0;
1853         vi->fno= 0;
1854         vi->mvert= 0;
1855         
1856         BLI_pbvh_node_get_grids(bvh, node, &grid_indices, &totgrid, NULL, &gridsize, &grids, NULL);
1857         BLI_pbvh_node_num_verts(bvh, node, &uniq_verts, &totvert);
1858         BLI_pbvh_node_get_verts(bvh, node, &vert_indices, &verts);
1859         
1860         vi->grids= grids;
1861         vi->grid_indices= grid_indices;
1862         vi->totgrid= (grids)? totgrid: 1;
1863         vi->gridsize= gridsize;
1864         
1865         if (mode == PBVH_ITER_ALL)
1866                 vi->totvert = totvert;
1867         else
1868                 vi->totvert= uniq_verts;
1869         vi->vert_indices= vert_indices;
1870         vi->mverts= verts;
1871
1872         vi->gh= NULL;
1873         if (vi->grids && mode == PBVH_ITER_UNIQUE)
1874                 vi->grid_hidden= bvh->grid_hidden;
1875 }