5c0b09f0ff0684503666907429b73dd8421c6889
[blender.git] / source / blender / blenkernel / intern / mball_tessellate.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  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * Contributor(s): Jiri Hnidek <jiri.hnidek@vslib.cz>.
22  *
23  * ***** END GPL LICENSE BLOCK *****
24  */
25
26 /** \file blender/blenkernel/intern/mball_tessellate.c
27  *  \ingroup bke
28  */
29
30 #include <stdio.h>
31 #include <string.h>
32 #include <math.h>
33 #include <stdlib.h>
34 #include <ctype.h>
35 #include <float.h>
36
37 #include "MEM_guardedalloc.h"
38
39 #include "DNA_object_types.h"
40 #include "DNA_meta_types.h"
41 #include "DNA_scene_types.h"
42
43 #include "BLI_listbase.h"
44 #include "BLI_math.h"
45 #include "BLI_string_utils.h"
46 #include "BLI_utildefines.h"
47 #include "BLI_memarena.h"
48
49 #include "BKE_global.h"
50
51 #include "BKE_depsgraph.h"
52 #include "BKE_scene.h"
53 #include "BKE_displist.h"
54 #include "BKE_mball_tessellate.h"  /* own include */
55
56 #include "BLI_strict_flags.h"
57
58 /* experimental (faster) normal calculation */
59 // #define USE_ACCUM_NORMAL
60
61 /* Data types */
62
63 typedef struct corner {         /* corner of a cube */
64         int i, j, k;                /* (i, j, k) is index within lattice */
65         float co[3], value;       /* location and function value */
66         struct corner *next;
67 } CORNER;
68
69 typedef struct cube {           /* partitioning cell (cube) */
70         int i, j, k;                /* lattice location of cube */
71         CORNER *corners[8];         /* eight corners */
72 } CUBE;
73
74 typedef struct cubes {          /* linked list of cubes acting as stack */
75         CUBE cube;                  /* a single cube */
76         struct cubes *next;         /* remaining elements */
77 } CUBES;
78
79 typedef struct centerlist {     /* list of cube locations */
80         int i, j, k;                /* cube location */
81         struct centerlist *next;    /* remaining elements */
82 } CENTERLIST;
83
84 typedef struct edgelist {       /* list of edges */
85         int i1, j1, k1, i2, j2, k2; /* edge corner ids */
86         int vid;                    /* vertex id */
87         struct edgelist *next;      /* remaining elements */
88 } EDGELIST;
89
90 typedef struct intlist {        /* list of integers */
91         int i;                      /* an integer */
92         struct intlist *next;       /* remaining elements */
93 } INTLIST;
94
95 typedef struct intlists {       /* list of list of integers */
96         INTLIST *list;              /* a list of integers */
97         struct intlists *next;      /* remaining elements */
98 } INTLISTS;
99
100 typedef struct Box {                    /* an AABB with pointer to metalelem */
101         float min[3], max[3];
102         const MetaElem *ml;
103 } Box;
104
105 typedef struct MetaballBVHNode {        /* BVH node */
106         Box bb[2];                                              /* AABB of children */
107         struct MetaballBVHNode *child[2];
108 } MetaballBVHNode;
109
110 typedef struct process {        /* parameters, storage */
111         float thresh, size;                     /* mball threshold, single cube size */
112         float delta;                            /* small delta for calculating normals */
113         unsigned int converge_res;      /* converge procedure resolution (more = slower) */
114
115         MetaElem **mainb;                       /* array of all metaelems */
116         unsigned int totelem, mem;      /* number of metaelems */
117
118         MetaballBVHNode metaball_bvh; /* The simplest bvh */
119         Box allbb;                   /* Bounding box of all metaelems */
120
121         MetaballBVHNode **bvh_queue; /* Queue used during bvh traversal */
122         unsigned int bvh_queue_size;
123
124         CUBES *cubes;               /* stack of cubes waiting for polygonization */
125         CENTERLIST **centers;       /* cube center hash table */
126         CORNER **corners;           /* corner value hash table */
127         EDGELIST **edges;           /* edge and vertex id hash table */
128
129         int (*indices)[4];          /* output indices */
130         unsigned int totindex;          /* size of memory allocated for indices */
131         unsigned int curindex;          /* number of currently added indices */
132
133         float (*co)[3], (*no)[3];   /* surface vertices - positions and normals */
134         unsigned int totvertex;         /* memory size */
135         unsigned int curvertex;         /* currently added vertices */
136
137         /* memory allocation from common pool */
138         MemArena *pgn_elements;
139 } PROCESS;
140
141 /* Forward declarations */
142 static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2);
143 static void add_cube(PROCESS *process, int i, int j, int k);
144 static void make_face(PROCESS *process, int i1, int i2, int i3, int i4);
145 static void converge(PROCESS *process, const CORNER *c1, const CORNER *c2, float r_p[3]);
146
147 /* ******************* SIMPLE BVH ********************* */
148
149 static void make_box_union(const BoundBox *a, const Box *b, Box *r_out)
150 {
151         r_out->min[0] = min_ff(a->vec[0][0], b->min[0]);
152         r_out->min[1] = min_ff(a->vec[0][1], b->min[1]);
153         r_out->min[2] = min_ff(a->vec[0][2], b->min[2]);
154
155         r_out->max[0] = max_ff(a->vec[6][0], b->max[0]);
156         r_out->max[1] = max_ff(a->vec[6][1], b->max[1]);
157         r_out->max[2] = max_ff(a->vec[6][2], b->max[2]);
158 }
159
160 static void make_box_from_metaelem(Box *r, const MetaElem *ml)
161 {
162         copy_v3_v3(r->max, ml->bb->vec[6]);
163         copy_v3_v3(r->min, ml->bb->vec[0]);
164         r->ml = ml;
165 }
166
167 /**
168  * Partitions part of mainb array [start, end) along axis s. Returns i,
169  * where centroids of elements in the [start, i) segment lie "on the right side" of div,
170  * and elements in the [i, end) segment lie "on the left"
171  */
172 static unsigned int partition_mainb(MetaElem **mainb, unsigned int start, unsigned int end, unsigned int s, float div)
173 {
174         unsigned int i = start, j = end - 1;
175         div *= 2.0f;
176
177         while (1) {
178                 while (i < j && div > (mainb[i]->bb->vec[6][s] + mainb[i]->bb->vec[0][s])) i++;
179                 while (j > i && div < (mainb[j]->bb->vec[6][s] + mainb[j]->bb->vec[0][s])) j--;
180
181                 if (i >= j)
182                         break;
183
184                 SWAP(MetaElem *, mainb[i], mainb[j]);
185                 i++;
186                 j--;
187         }
188
189         if (i == start) {
190                 i++;
191         }
192
193         return i;
194 }
195
196 /**
197  * Recursively builds a BVH, dividing elements along the middle of the longest axis of allbox.
198  */
199 static void build_bvh_spatial(
200         PROCESS *process, MetaballBVHNode *node,
201         unsigned int start, unsigned int end, const Box *allbox)
202 {
203         unsigned int part, j, s;
204         float dim[3], div;
205
206         /* Maximum bvh queue size is number of nodes which are made, equals calls to this function. */
207         process->bvh_queue_size++;
208
209         dim[0] = allbox->max[0] - allbox->min[0];
210         dim[1] = allbox->max[1] - allbox->min[1];
211         dim[2] = allbox->max[2] - allbox->min[2];
212
213         s = 0;
214         if (dim[1] > dim[0] && dim[1] > dim[2]) s = 1;
215         else if (dim[2] > dim[1] && dim[2] > dim[0]) s = 2;
216
217         div = allbox->min[s] + (dim[s] / 2.0f);
218
219         part = partition_mainb(process->mainb, start, end, s, div);
220
221         make_box_from_metaelem(&node->bb[0], process->mainb[start]);
222         node->child[0] = NULL;
223
224         if (part > start + 1) {
225                 for (j = start; j < part; j++) {
226                         make_box_union(process->mainb[j]->bb, &node->bb[0], &node->bb[0]);
227                 }
228
229                 node->child[0] = BLI_memarena_alloc(process->pgn_elements, sizeof(MetaballBVHNode));
230                 build_bvh_spatial(process, node->child[0], start, part, &node->bb[0]);
231         }
232
233         node->child[1] = NULL;
234         if (part < end) {
235                 make_box_from_metaelem(&node->bb[1], process->mainb[part]);
236
237                 if (part < end - 1) {
238                         for (j = part; j < end; j++) {
239                                 make_box_union(process->mainb[j]->bb, &node->bb[1], &node->bb[1]);
240                         }
241
242                         node->child[1] = BLI_memarena_alloc(process->pgn_elements, sizeof(MetaballBVHNode));
243                         build_bvh_spatial(process, node->child[1], part, end, &node->bb[1]);
244                 }
245         }
246         else {
247                 INIT_MINMAX(node->bb[1].min, node->bb[1].max);
248         }
249 }
250
251 /* ******************** ARITH ************************* */
252
253 /**
254  * BASED AT CODE (but mostly rewritten) :
255  * C code from the article
256  * "An Implicit Surface Polygonizer"
257  * by Jules Bloomenthal, jbloom@beauty.gmu.edu
258  * in "Graphics Gems IV", Academic Press, 1994
259  *
260  * Authored by Jules Bloomenthal, Xerox PARC.
261  * Copyright (c) Xerox Corporation, 1991.  All rights reserved.
262  * Permission is granted to reproduce, use and distribute this code for
263  * any and all purposes, provided that this notice appears in all copies.
264  */
265
266 #define L   0  /* left direction:       -x, -i */
267 #define R   1  /* right direction:      +x, +i */
268 #define B   2  /* bottom direction: -y, -j */
269 #define T   3  /* top direction:        +y, +j */
270 #define N   4  /* near direction:       -z, -k */
271 #define F   5  /* far direction:        +z, +k */
272 #define LBN 0  /* left bottom near corner  */
273 #define LBF 1  /* left bottom far corner   */
274 #define LTN 2  /* left top near corner     */
275 #define LTF 3  /* left top far corner      */
276 #define RBN 4  /* right bottom near corner */
277 #define RBF 5  /* right bottom far corner  */
278 #define RTN 6  /* right top near corner    */
279 #define RTF 7  /* right top far corner     */
280
281 /**
282  * the LBN corner of cube (i, j, k), corresponds with location
283  * (i-0.5)*size, (j-0.5)*size, (k-0.5)*size)
284  */
285
286 #define HASHBIT     (5)
287 #define HASHSIZE    (size_t)(1 << (3 * HASHBIT))   /*! < hash table size (32768) */
288
289 #define HASH(i, j, k) ((((( (i) & 31) << 5) | ( (j) & 31)) << 5) | ( (k) & 31) )
290
291 #define MB_BIT(i, bit) (((i) >> (bit)) & 1)
292 // #define FLIP(i, bit) ((i) ^ 1 << (bit)) /* flip the given bit of i */
293
294 /* ******************** DENSITY COPMPUTATION ********************* */
295
296 /**
297  * Computes density from given metaball at given position.
298  * Metaball equation is: ``(1 - r^2 / R^2)^3 * s``
299  *
300  * r = distance from center
301  * R = metaball radius
302  * s - metaball stiffness
303  */
304 static float densfunc(const MetaElem *ball, float x, float y, float z)
305 {
306         float dist2;
307         float dvec[3] = {x, y, z};
308
309         mul_m4_v3((float (*)[4])ball->imat, dvec);
310
311         switch (ball->type) {
312                 case MB_BALL:
313                         /* do nothing */
314                         break;
315                 case MB_CUBE:
316                         if      (dvec[2] > ball->expz)  dvec[2] -= ball->expz;
317                         else if (dvec[2] < -ball->expz) dvec[2] += ball->expz;
318                         else                            dvec[2] = 0.0;
319                         /* fall through */
320                 case MB_PLANE:
321                         if      (dvec[1] >  ball->expy) dvec[1] -= ball->expy;
322                         else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
323                         else                            dvec[1] = 0.0;
324                         /* fall through */
325                 case MB_TUBE:
326                         if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
327                         else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
328                         else                            dvec[0] = 0.0;
329                         break;
330                 case MB_ELIPSOID:
331                         dvec[0] /= ball->expx;
332                         dvec[1] /= ball->expy;
333                         dvec[2] /= ball->expz;
334                         break;
335
336                 /* *** deprecated, could be removed?, do-versioned at least *** */
337                 case MB_TUBEX:
338                         if      (dvec[0] >  ball->len) dvec[0] -= ball->len;
339                         else if (dvec[0] < -ball->len) dvec[0] += ball->len;
340                         else                           dvec[0] = 0.0;
341                         break;
342                 case MB_TUBEY:
343                         if      (dvec[1] >  ball->len) dvec[1] -= ball->len;
344                         else if (dvec[1] < -ball->len) dvec[1] += ball->len;
345                         else                           dvec[1] = 0.0;
346                         break;
347                 case MB_TUBEZ:
348                         if      (dvec[2] >  ball->len) dvec[2] -= ball->len;
349                         else if (dvec[2] < -ball->len) dvec[2] += ball->len;
350                         else                           dvec[2] = 0.0;
351                         break;
352                         /* *** end deprecated *** */
353         }
354
355         /* ball->rad2 is inverse of squared rad */
356         dist2 = 1.0f - (len_squared_v3(dvec) * ball->rad2);
357
358         /* ball->s is negative if metaball is negative */
359         return (dist2 < 0.0f) ? 0.0f : (ball->s * dist2 * dist2 * dist2);
360 }
361
362 /**
363  * Computes density at given position form all metaballs which contain this point in their box.
364  * Traverses BVH using a queue.
365  */
366 static float metaball(PROCESS *process, float x, float y, float z)
367 {
368         int i;
369         float dens = 0.0f;
370         unsigned int front = 0, back = 0;
371         MetaballBVHNode *node;
372
373         process->bvh_queue[front++] = &process->metaball_bvh;
374
375         while (front != back) {
376                 node = process->bvh_queue[back++];
377
378                 for (i = 0; i < 2; i++) {
379                         if ((node->bb[i].min[0] <= x) && (node->bb[i].max[0] >= x) &&
380                             (node->bb[i].min[1] <= y) && (node->bb[i].max[1] >= y) &&
381                             (node->bb[i].min[2] <= z) && (node->bb[i].max[2] >= z))
382                         {
383                                 if (node->child[i])     process->bvh_queue[front++] = node->child[i];
384                                 else dens += densfunc(node->bb[i].ml, x, y, z);
385                         }
386                 }
387         }
388
389         return process->thresh - dens;
390 }
391
392 /**
393  * Adds face to indices, expands memory if needed.
394  */
395 static void make_face(PROCESS *process, int i1, int i2, int i3, int i4)
396 {
397         int *cur;
398
399 #ifdef USE_ACCUM_NORMAL
400         float n[3];
401 #endif
402
403         if (UNLIKELY(process->totindex == process->curindex)) {
404                 process->totindex += 4096;
405                 process->indices = MEM_reallocN(process->indices, sizeof(int[4]) * process->totindex);
406         }
407
408         cur = process->indices[process->curindex++];
409
410         /* displists now support array drawing, we treat tri's as fake quad */
411
412         cur[0] = i1;
413         cur[1] = i2;
414         cur[2] = i3;
415
416         if (i4 == 0) {
417                 cur[3] = i3;
418         }
419         else {
420                 cur[3] = i4;
421         }
422
423 #ifdef USE_ACCUM_NORMAL
424         if (i4 == 0) {
425                 normal_tri_v3(n, process->co[i1], process->co[i2], process->co[i3]);
426                 accumulate_vertex_normals(
427                         process->no[i1], process->no[i2], process->no[i3], NULL, n,
428                         process->co[i1], process->co[i2], process->co[i3], NULL);
429         }
430         else {
431                 normal_quad_v3(n, process->co[i1], process->co[i2], process->co[i3], process->co[i4]);
432                 accumulate_vertex_normals(
433                         process->no[i1], process->no[i2], process->no[i3], process->no[i4], n,
434                         process->co[i1], process->co[i2], process->co[i3], process->co[i4]);
435         }
436 #endif
437
438 }
439
440 /* Frees allocated memory */
441 static void freepolygonize(PROCESS *process)
442 {
443         if (process->corners) MEM_freeN(process->corners);
444         if (process->edges) MEM_freeN(process->edges);
445         if (process->centers) MEM_freeN(process->centers);
446         if (process->mainb) MEM_freeN(process->mainb);
447         if (process->bvh_queue) MEM_freeN(process->bvh_queue);
448         if (process->pgn_elements) BLI_memarena_free(process->pgn_elements);
449 }
450
451 /* **************** POLYGONIZATION ************************ */
452
453 /**** Cubical Polygonization (optional) ****/
454
455 #define LB  0  /* left bottom edge      */
456 #define LT  1  /* left top edge */
457 #define LN  2  /* left near edge        */
458 #define LF  3  /* left far edge */
459 #define RB  4  /* right bottom edge */
460 #define RT  5  /* right top edge        */
461 #define RN  6  /* right near edge       */
462 #define RF  7  /* right far edge        */
463 #define BN  8  /* bottom near edge      */
464 #define BF  9  /* bottom far edge       */
465 #define TN  10 /* top near edge */
466 #define TF  11 /* top far edge  */
467
468 static INTLISTS *cubetable[256];
469 static char faces[256];
470
471 /* edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
472 static int corner1[12] = {
473         LBN, LTN, LBN, LBF, RBN, RTN, RBN, RBF, LBN, LBF, LTN, LTF
474 };
475 static int corner2[12] = {
476         LBF, LTF, LTN, LTF, RBF, RTF, RTN, RTF, RBN, RBF, RTN, RTF
477 };
478 static int leftface[12] = {
479         B, L, L, F, R, T, N, R, N, B, T, F
480 };
481 /* face on left when going corner1 to corner2 */
482 static int rightface[12] = {
483         L, T, N, L, B, R, R, F, B, F, N, T
484 };
485 /* face on right when going corner1 to corner2 */
486
487 /**
488  * triangulate the cube directly, without decomposition
489  */
490 static void docube(PROCESS *process, CUBE *cube)
491 {
492         INTLISTS *polys;
493         CORNER *c1, *c2;
494         int i, index = 0, count, indexar[8];
495
496         /* Determine which case cube falls into. */
497         for (i = 0; i < 8; i++) {
498                 if (cube->corners[i]->value > 0.0f) {
499                         index += (1 << i);
500                 }
501         }
502
503         /* Using faces[] table, adds neighbouring cube if surface intersects face in this direction. */
504         if (MB_BIT(faces[index], 0)) add_cube(process, cube->i - 1, cube->j, cube->k);
505         if (MB_BIT(faces[index], 1)) add_cube(process, cube->i + 1, cube->j, cube->k);
506         if (MB_BIT(faces[index], 2)) add_cube(process, cube->i, cube->j - 1, cube->k);
507         if (MB_BIT(faces[index], 3)) add_cube(process, cube->i, cube->j + 1, cube->k);
508         if (MB_BIT(faces[index], 4)) add_cube(process, cube->i, cube->j, cube->k - 1);
509         if (MB_BIT(faces[index], 5)) add_cube(process, cube->i, cube->j, cube->k + 1);
510
511         /* Using cubetable[], determines polygons for output. */
512         for (polys = cubetable[index]; polys; polys = polys->next) {
513                 INTLIST *edges;
514
515                 count = 0;
516                 /* Sets needed vertex id's lying on the edges. */
517                 for (edges = polys->list; edges; edges = edges->next) {
518                         c1 = cube->corners[corner1[edges->i]];
519                         c2 = cube->corners[corner2[edges->i]];
520
521                         indexar[count] = vertid(process, c1, c2);
522                         count++;
523                 }
524
525                 /* Adds faces to output. */
526                 if (count > 2) {
527                         switch (count) {
528                                 case 3:
529                                         make_face(process, indexar[2], indexar[1], indexar[0], 0);
530                                         break;
531                                 case 4:
532                                         if (indexar[0] == 0) make_face(process, indexar[0], indexar[3], indexar[2], indexar[1]);
533                                         else make_face(process, indexar[3], indexar[2], indexar[1], indexar[0]);
534                                         break;
535                                 case 5:
536                                         if (indexar[0] == 0) make_face(process, indexar[0], indexar[3], indexar[2], indexar[1]);
537                                         else make_face(process, indexar[3], indexar[2], indexar[1], indexar[0]);
538
539                                         make_face(process, indexar[4], indexar[3], indexar[0], 0);
540                                         break;
541                                 case 6:
542                                         if (indexar[0] == 0) {
543                                                 make_face(process, indexar[0], indexar[3], indexar[2], indexar[1]);
544                                                 make_face(process, indexar[0], indexar[5], indexar[4], indexar[3]);
545                                         }
546                                         else {
547                                                 make_face(process, indexar[3], indexar[2], indexar[1], indexar[0]);
548                                                 make_face(process, indexar[5], indexar[4], indexar[3], indexar[0]);
549                                         }
550                                         break;
551                                 case 7:
552                                         if (indexar[0] == 0) {
553                                                 make_face(process, indexar[0], indexar[3], indexar[2], indexar[1]);
554                                                 make_face(process, indexar[0], indexar[5], indexar[4], indexar[3]);
555                                         }
556                                         else {
557                                                 make_face(process, indexar[3], indexar[2], indexar[1], indexar[0]);
558                                                 make_face(process, indexar[5], indexar[4], indexar[3], indexar[0]);
559                                         }
560
561                                         make_face(process, indexar[6], indexar[5], indexar[0], 0);
562
563                                         break;
564                         }
565                 }
566         }
567 }
568
569 /**
570  * return corner with the given lattice location
571  * set (and cache) its function value
572  */
573 static CORNER *setcorner(PROCESS *process, int i, int j, int k)
574 {
575         /* for speed, do corner value caching here */
576         CORNER *c;
577         int index;
578
579         /* does corner exist? */
580         index = HASH(i, j, k);
581         c = process->corners[index];
582
583         for (; c != NULL; c = c->next) {
584                 if (c->i == i && c->j == j && c->k == k) {
585                         return c;
586                 }
587         }
588
589         c = BLI_memarena_alloc(process->pgn_elements, sizeof(CORNER));
590
591         c->i = i;
592         c->co[0] = ((float)i - 0.5f) * process->size;
593         c->j = j;
594         c->co[1] = ((float)j - 0.5f) * process->size;
595         c->k = k;
596         c->co[2] = ((float)k - 0.5f) * process->size;
597
598         c->value = metaball(process, c->co[0], c->co[1], c->co[2]);
599
600         c->next = process->corners[index];
601         process->corners[index] = c;
602
603         return c;
604 }
605
606 /**
607  * return next clockwise edge from given edge around given face
608  */
609 static int nextcwedge(int edge, int face)
610 {
611         switch (edge) {
612                 case LB:
613                         return (face == L) ? LF : BN;
614                 case LT:
615                         return (face == L) ? LN : TF;
616                 case LN:
617                         return (face == L) ? LB : TN;
618                 case LF:
619                         return (face == L) ? LT : BF;
620                 case RB:
621                         return (face == R) ? RN : BF;
622                 case RT:
623                         return (face == R) ? RF : TN;
624                 case RN:
625                         return (face == R) ? RT : BN;
626                 case RF:
627                         return (face == R) ? RB : TF;
628                 case BN:
629                         return (face == B) ? RB : LN;
630                 case BF:
631                         return (face == B) ? LB : RF;
632                 case TN:
633                         return (face == T) ? LT : RN;
634                 case TF:
635                         return (face == T) ? RT : LF;
636         }
637         return 0;
638 }
639
640 /**
641  * \return the face adjoining edge that is not the given face
642  */
643 static int otherface(int edge, int face)
644 {
645         int other = leftface[edge];
646         return face == other ? rightface[edge] : other;
647 }
648
649 /**
650  * create the 256 entry table for cubical polygonization
651  */
652 static void makecubetable(void)
653 {
654         static bool is_done = false;
655         int i, e, c, done[12], pos[8];
656
657         if (is_done) return;
658         is_done = true;
659
660         for (i = 0; i < 256; i++) {
661                 for (e = 0; e < 12; e++) done[e] = 0;
662                 for (c = 0; c < 8; c++) pos[c] = MB_BIT(i, c);
663                 for (e = 0; e < 12; e++)
664                         if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
665                                 INTLIST *ints = NULL;
666                                 INTLISTS *lists = MEM_callocN(sizeof(INTLISTS), "mball_intlist");
667                                 int start = e, edge = e;
668
669                                 /* get face that is to right of edge from pos to neg corner: */
670                                 int face = pos[corner1[e]] ? rightface[e] : leftface[e];
671
672                                 while (1) {
673                                         edge = nextcwedge(edge, face);
674                                         done[edge] = 1;
675                                         if (pos[corner1[edge]] != pos[corner2[edge]]) {
676                                                 INTLIST *tmp = ints;
677
678                                                 ints = MEM_callocN(sizeof(INTLIST), "mball_intlist");
679                                                 ints->i = edge;
680                                                 ints->next = tmp; /* add edge to head of list */
681
682                                                 if (edge == start) break;
683                                                 face = otherface(edge, face);
684                                         }
685                                 }
686                                 lists->list = ints; /* add ints to head of table entry */
687                                 lists->next = cubetable[i];
688                                 cubetable[i] = lists;
689                         }
690         }
691
692         for (i = 0; i < 256; i++) {
693                 INTLISTS *polys;
694                 faces[i] = 0;
695                 for (polys = cubetable[i]; polys; polys = polys->next) {
696                         INTLIST *edges;
697
698                         for (edges = polys->list; edges; edges = edges->next) {
699                                 if (edges->i == LB || edges->i == LT || edges->i == LN || edges->i == LF) faces[i] |= 1 << L;
700                                 if (edges->i == RB || edges->i == RT || edges->i == RN || edges->i == RF) faces[i] |= 1 << R;
701                                 if (edges->i == LB || edges->i == RB || edges->i == BN || edges->i == BF) faces[i] |= 1 << B;
702                                 if (edges->i == LT || edges->i == RT || edges->i == TN || edges->i == TF) faces[i] |= 1 << T;
703                                 if (edges->i == LN || edges->i == RN || edges->i == BN || edges->i == TN) faces[i] |= 1 << N;
704                                 if (edges->i == LF || edges->i == RF || edges->i == BF || edges->i == TF) faces[i] |= 1 << F;
705                         }
706                 }
707         }
708 }
709
710 void BKE_mball_cubeTable_free(void)
711 {
712         int i;
713         INTLISTS *lists, *nlists;
714         INTLIST *ints, *nints;
715
716         for (i = 0; i < 256; i++) {
717                 lists = cubetable[i];
718                 while (lists) {
719                         nlists = lists->next;
720
721                         ints = lists->list;
722                         while (ints) {
723                                 nints = ints->next;
724                                 MEM_freeN(ints);
725                                 ints = nints;
726                         }
727
728                         MEM_freeN(lists);
729                         lists = nlists;
730                 }
731                 cubetable[i] = NULL;
732         }
733 }
734
735 /**** Storage ****/
736
737 /**
738  * Inserts cube at lattice i, j, k into hash table, marking it as "done"
739  */
740 static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const int j, const int k)
741 {
742         int index;
743         CENTERLIST *newc, *l, *q;
744
745         index = HASH(i, j, k);
746         q = table[index];
747
748         for (l = q; l != NULL; l = l->next) {
749                 if (l->i == i && l->j == j && l->k == k) return 1;
750         }
751
752         newc = BLI_memarena_alloc(process->pgn_elements, sizeof(CENTERLIST));
753         newc->i = i;
754         newc->j = j;
755         newc->k = k;
756         newc->next = q;
757         table[index] = newc;
758
759         return 0;
760 }
761
762 /**
763  * Sets vid of vertex lying on given edge.
764  */
765 static void setedge(
766         PROCESS *process,
767         int i1, int j1, int k1,
768         int i2, int j2, int k2,
769         int vid)
770 {
771         int index;
772         EDGELIST *newe;
773
774         if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
775                 int t = i1;
776                 i1 = i2;
777                 i2 = t;
778                 t = j1;
779                 j1 = j2;
780                 j2 = t;
781                 t = k1;
782                 k1 = k2;
783                 k2 = t;
784         }
785         index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
786         newe = BLI_memarena_alloc(process->pgn_elements, sizeof(EDGELIST));
787
788         newe->i1 = i1;
789         newe->j1 = j1;
790         newe->k1 = k1;
791         newe->i2 = i2;
792         newe->j2 = j2;
793         newe->k2 = k2;
794         newe->vid = vid;
795         newe->next = process->edges[index];
796         process->edges[index] = newe;
797 }
798
799 /**
800  * \return vertex id for edge; return -1 if not set
801  */
802 static int getedge(EDGELIST *table[],
803                    int i1, int j1, int k1,
804                    int i2, int j2, int k2)
805 {
806         EDGELIST *q;
807
808         if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
809                 int t = i1;
810                 i1 = i2;
811                 i2 = t;
812                 t = j1;
813                 j1 = j2;
814                 j2 = t;
815                 t = k1;
816                 k1 = k2;
817                 k2 = t;
818         }
819         q = table[HASH(i1, j1, k1) + HASH(i2, j2, k2)];
820         for (; q != NULL; q = q->next) {
821                 if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 &&
822                     q->i2 == i2 && q->j2 == j2 && q->k2 == k2)
823                 {
824                         return q->vid;
825                 }
826         }
827         return -1;
828 }
829
830 /**
831  * Adds a vertex, expands memory if needed.
832  */
833 static void addtovertices(PROCESS *process, const float v[3], const float no[3])
834 {
835         if (process->curvertex == process->totvertex) {
836                 process->totvertex += 4096;
837                 process->co = MEM_reallocN(process->co, process->totvertex * sizeof(float[3]));
838                 process->no = MEM_reallocN(process->no, process->totvertex * sizeof(float[3]));
839         }
840
841         copy_v3_v3(process->co[process->curvertex], v);
842         copy_v3_v3(process->no[process->curvertex], no);
843
844         process->curvertex++;
845 }
846
847 #ifndef USE_ACCUM_NORMAL
848 /**
849  * Computes normal from density field at given point.
850  *
851  * \note Doesn't do normalization!
852  */
853 static void vnormal(PROCESS *process, const float point[3], float r_no[3])
854 {
855         const float delta = process->delta;
856         const float f = metaball(process, point[0], point[1], point[2]);
857
858         r_no[0] = metaball(process, point[0] + delta, point[1], point[2]) - f;
859         r_no[1] = metaball(process, point[0], point[1] + delta, point[2]) - f;
860         r_no[2] = metaball(process, point[0], point[1], point[2] + delta) - f;
861
862 #if 0
863         f = normalize_v3(r_no);
864
865         if (0) {
866                 float tvec[3];
867
868                 delta *= 2.0f;
869
870                 f = process->function(process, point[0], point[1], point[2]);
871
872                 tvec[0] = process->function(process, point[0] + delta, point[1], point[2]) - f;
873                 tvec[1] = process->function(process, point[0], point[1] + delta, point[2]) - f;
874                 tvec[2] = process->function(process, point[0], point[1], point[2] + delta) - f;
875
876                 if (normalize_v3(tvec) != 0.0f) {
877                         add_v3_v3(r_no, tvec);
878                         normalize_v3(r_no);
879                 }
880         }
881 #endif
882 }
883 #endif  /* USE_ACCUM_NORMAL */
884
885 /**
886  * \return the id of vertex between two corners.
887  *
888  * If it wasn't previously computed, does #converge() and adds vertex to process.
889  */
890 static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2)
891 {
892         float v[3], no[3];
893         int vid = getedge(process->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
894
895         if (vid != -1) return vid;  /* previously computed */
896
897         converge(process, c1, c2, v);  /* position */
898
899 #ifdef USE_ACCUM_NORMAL
900         zero_v3(no);
901 #else
902         vnormal(process, v, no);
903 #endif
904
905         addtovertices(process, v, no);            /* save vertex */
906         vid = (int)process->curvertex - 1;
907         setedge(process, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
908
909         return vid;
910 }
911
912 /**
913  * Given two corners, computes approximation of surface intersection point between them.
914  * In case of small threshold, do bisection.
915  */
916 static void converge(PROCESS *process, const CORNER *c1, const CORNER *c2, float r_p[3])
917 {
918         float tmp, dens;
919         unsigned int i;
920         float c1_value, c1_co[3];
921         float c2_value, c2_co[3];
922
923         if (c1->value < c2->value) {
924                 c1_value = c2->value;
925                 copy_v3_v3(c1_co, c2->co);
926                 c2_value = c1->value;
927                 copy_v3_v3(c2_co, c1->co);
928         }
929         else {
930                 c1_value = c1->value;
931                 copy_v3_v3(c1_co, c1->co);
932                 c2_value = c2->value;
933                 copy_v3_v3(c2_co, c2->co);
934         }
935
936
937         for (i = 0; i < process->converge_res; i++) {
938                 interp_v3_v3v3(r_p, c1_co, c2_co, 0.5f);
939                 dens = metaball(process, r_p[0], r_p[1], r_p[2]);
940
941                 if (dens > 0.0f) {
942                         c1_value = dens;
943                         copy_v3_v3(c1_co, r_p);
944                 }
945                 else {
946                         c2_value = dens;
947                         copy_v3_v3(c2_co, r_p);
948                 }
949         }
950
951         tmp = -c1_value / (c2_value - c1_value);
952         interp_v3_v3v3(r_p, c1_co, c2_co, tmp);
953 }
954
955 /**
956  * Adds cube at given lattice position to cube stack of process.
957  */
958 static void add_cube(PROCESS *process, int i, int j, int k)
959 {
960         CUBES *ncube;
961         int n;
962
963         /* test if cube has been found before */
964         if (setcenter(process, process->centers, i, j, k) == 0) {
965                 /* push cube on stack: */
966                 ncube = BLI_memarena_alloc(process->pgn_elements, sizeof(CUBES));
967                 ncube->next = process->cubes;
968                 process->cubes = ncube;
969
970                 ncube->cube.i = i;
971                 ncube->cube.j = j;
972                 ncube->cube.k = k;
973
974                 /* set corners of initial cube: */
975                 for (n = 0; n < 8; n++)
976                         ncube->cube.corners[n] = setcorner(process, i + MB_BIT(n, 2), j + MB_BIT(n, 1), k + MB_BIT(n, 0));
977         }
978 }
979
980 static void next_lattice(int r[3], const float pos[3], const float size)
981 {
982         r[0] = (int)ceil((pos[0] / size) + 0.5f);
983         r[1] = (int)ceil((pos[1] / size) + 0.5f);
984         r[2] = (int)ceil((pos[2] / size) + 0.5f);
985 }
986 static void prev_lattice(int r[3], const float pos[3], const float size)
987 {
988         next_lattice(r, pos, size);
989         r[0]--; r[1]--; r[2]--;
990 }
991 static void closest_latice(int r[3], const float pos[3], const float size)
992 {
993         r[0] = (int)floorf(pos[0] / size + 1.0f);
994         r[1] = (int)floorf(pos[1] / size + 1.0f);
995         r[2] = (int)floorf(pos[2] / size + 1.0f);
996 }
997
998 /**
999  * Find at most 26 cubes to start polygonization from.
1000  */
1001 static void find_first_points(PROCESS *process, const unsigned int em)
1002 {
1003         const MetaElem *ml;
1004         int center[3], lbn[3], rtf[3], it[3], dir[3], add[3];
1005         float tmp[3], a, b;
1006
1007         ml = process->mainb[em];
1008
1009         mid_v3_v3v3(tmp, ml->bb->vec[0], ml->bb->vec[6]);
1010         closest_latice(center, tmp, process->size);
1011         prev_lattice(lbn, ml->bb->vec[0], process->size);
1012         next_lattice(rtf, ml->bb->vec[6], process->size);
1013
1014         for (dir[0] = -1; dir[0] <= 1; dir[0]++) {
1015                 for (dir[1] = -1; dir[1] <= 1; dir[1]++) {
1016                         for (dir[2] = -1; dir[2] <= 1; dir[2]++) {
1017                                 if (dir[0] == 0 && dir[1] == 0 && dir[2] == 0) {
1018                                         continue;
1019                                 }
1020
1021                                 copy_v3_v3_int(it, center);
1022
1023                                 b = setcorner(process, it[0], it[1], it[2])->value;
1024                                 do {
1025                                         it[0] += dir[0];
1026                                         it[1] += dir[1];
1027                                         it[2] += dir[2];
1028                                         a = b;
1029                                         b = setcorner(process, it[0], it[1], it[2])->value;
1030
1031                                         if (a * b < 0.0f) {
1032                                                 add[0] = it[0] - dir[0];
1033                                                 add[1] = it[1] - dir[1];
1034                                                 add[2] = it[2] - dir[2];
1035                                                 DO_MIN(it, add);
1036                                                 add_cube(process, add[0], add[1], add[2]);
1037                                                 break;
1038                                         }
1039                                 } while ((it[0] > lbn[0]) && (it[1] > lbn[1]) && (it[2] > lbn[2]) &&
1040                                          (it[0] < rtf[0]) && (it[1] < rtf[1]) && (it[2] < rtf[2]));
1041                         }
1042                 }
1043         }
1044 }
1045
1046 /**
1047  * The main polygonization proc.
1048  * Allocates memory, makes cubetable,
1049  * finds starting surface points
1050  * and processes cubes on the stack until none left.
1051  */
1052 static void polygonize(PROCESS *process)
1053 {
1054         CUBE c;
1055         unsigned int i;
1056
1057         process->centers = MEM_callocN(HASHSIZE * sizeof(CENTERLIST *), "mbproc->centers");
1058         process->corners = MEM_callocN(HASHSIZE * sizeof(CORNER *), "mbproc->corners");
1059         process->edges = MEM_callocN(2 * HASHSIZE * sizeof(EDGELIST *), "mbproc->edges");
1060         process->bvh_queue = MEM_callocN(sizeof(MetaballBVHNode *) * process->bvh_queue_size, "Metaball BVH Queue");
1061
1062         makecubetable();
1063
1064         for (i = 0; i < process->totelem; i++) {
1065                 find_first_points(process, i);
1066         }
1067
1068         while (process->cubes != NULL) {
1069                 c = process->cubes->cube;
1070                 process->cubes = process->cubes->next;
1071
1072                 docube(process, &c);
1073         }
1074 }
1075
1076 /**
1077  * Iterates over ALL objects in the scene and all of its sets, including
1078  * making all duplis(not only metas). Copies metas to mainb array.
1079  * Computes bounding boxes for building BVH. */
1080 static void init_meta(EvaluationContext *eval_ctx, PROCESS *process, Scene *scene, Object *ob)
1081 {
1082         Scene *sce_iter = scene;
1083         Base *base;
1084         Object *bob;
1085         MetaBall *mb;
1086         const MetaElem *ml;
1087         float obinv[4][4], obmat[4][4];
1088         unsigned int i;
1089         int obnr, zero_size = 0;
1090         char obname[MAX_ID_NAME];
1091         SceneBaseIter iter;
1092
1093         copy_m4_m4(obmat, ob->obmat);   /* to cope with duplicators from BKE_scene_base_iter_next */
1094         invert_m4_m4(obinv, ob->obmat);
1095
1096         BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
1097
1098         /* make main array */
1099         BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 0, NULL, NULL);
1100         while (BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 1, &base, &bob)) {
1101                 if (bob->type == OB_MBALL) {
1102                         zero_size = 0;
1103                         ml = NULL;
1104
1105                         if (bob == ob && (base->flag & OB_FROMDUPLI) == 0) {
1106                                 mb = ob->data;
1107
1108                                 if (mb->editelems) ml = mb->editelems->first;
1109                                 else ml = mb->elems.first;
1110                         }
1111                         else {
1112                                 char name[MAX_ID_NAME];
1113                                 int nr;
1114
1115                                 BLI_split_name_num(name, &nr, bob->id.name + 2, '.');
1116                                 if (STREQ(obname, name)) {
1117                                         mb = bob->data;
1118
1119                                         if (mb->editelems) ml = mb->editelems->first;
1120                                         else ml = mb->elems.first;
1121                                 }
1122                         }
1123
1124                         /* when metaball object has zero scale, then MetaElem to this MetaBall
1125                          * will not be put to mainb array */
1126                         if (has_zero_axis_m4(bob->obmat)) {
1127                                 zero_size = 1;
1128                         }
1129                         else if (bob->parent) {
1130                                 struct Object *pob = bob->parent;
1131                                 while (pob) {
1132                                         if (has_zero_axis_m4(pob->obmat)) {
1133                                                 zero_size = 1;
1134                                                 break;
1135                                         }
1136                                         pob = pob->parent;
1137                                 }
1138                         }
1139
1140                         if (zero_size) {
1141                                 while (ml) {
1142                                         ml = ml->next;
1143                                 }
1144                         }
1145                         else {
1146                                 while (ml) {
1147                                         if (!(ml->flag & MB_HIDE)) {
1148                                                 float pos[4][4], rot[4][4];
1149                                                 float expx, expy, expz;
1150                                                 float tempmin[3], tempmax[3];
1151
1152                                                 MetaElem *new_ml;
1153
1154                                                 /* make a copy because of duplicates */
1155                                                 new_ml = BLI_memarena_alloc(process->pgn_elements, sizeof(MetaElem));
1156                                                 *(new_ml) = *ml;
1157                                                 new_ml->bb = BLI_memarena_alloc(process->pgn_elements, sizeof(BoundBox));
1158                                                 new_ml->mat = BLI_memarena_alloc(process->pgn_elements, 4 * 4 * sizeof(float));
1159                                                 new_ml->imat = BLI_memarena_alloc(process->pgn_elements, 4 * 4 * sizeof(float));
1160
1161                                                 /* too big stiffness seems only ugly due to linear interpolation
1162                                                  * no need to have possibility for too big stiffness */
1163                                                 if (ml->s > 10.0f) new_ml->s = 10.0f;
1164                                                 else new_ml->s = ml->s;
1165
1166                                                 /* if metaball is negative, set stiffness negative */
1167                                                 if (new_ml->flag & MB_NEGATIVE) new_ml->s = -new_ml->s;
1168
1169                                                 /* Translation of MetaElem */
1170                                                 unit_m4(pos);
1171                                                 pos[3][0] = ml->x;
1172                                                 pos[3][1] = ml->y;
1173                                                 pos[3][2] = ml->z;
1174
1175                                                 /* Rotation of MetaElem is stored in quat */
1176                                                 quat_to_mat4(rot, ml->quat);
1177
1178                                                 /* basis object space -> world -> ml object space -> position -> rotation -> ml local space */
1179                                                 mul_m4_series((float(*)[4])new_ml->mat, obinv, bob->obmat, pos, rot);
1180                                                 /* ml local space -> basis object space */
1181                                                 invert_m4_m4((float(*)[4])new_ml->imat, (float(*)[4])new_ml->mat);
1182
1183                                                 /* rad2 is inverse of squared radius */
1184                                                 new_ml->rad2 = 1 / (ml->rad * ml->rad);
1185
1186                                                 /* initial dimensions = radius */
1187                                                 expx = ml->rad;
1188                                                 expy = ml->rad;
1189                                                 expz = ml->rad;
1190
1191                                                 switch (ml->type) {
1192                                                         case MB_BALL:
1193                                                                 break;
1194                                                         case MB_CUBE: /* cube is "expanded" by expz, expy and expx */
1195                                                                 expz += ml->expz;
1196                                                                 /* fall through */
1197                                                         case MB_PLANE: /* plane is "expanded" by expy and expx */
1198                                                                 expy += ml->expy;
1199                                                                 /* fall through */
1200                                                         case MB_TUBE: /* tube is "expanded" by expx */
1201                                                                 expx += ml->expx;
1202                                                                 break;
1203                                                         case MB_ELIPSOID: /* ellipsoid is "stretched" by exp* */
1204                                                                 expx *= ml->expx;
1205                                                                 expy *= ml->expy;
1206                                                                 expz *= ml->expz;
1207                                                                 break;
1208                                                 }
1209
1210                                                 /* untransformed Bounding Box of MetaElem */
1211                                                 /* TODO, its possible the elem type has been changed and the exp* values can use a fallback */
1212                                                 copy_v3_fl3(new_ml->bb->vec[0], -expx, -expy, -expz);  /* 0 */
1213                                                 copy_v3_fl3(new_ml->bb->vec[1], +expx, -expy, -expz);  /* 1 */
1214                                                 copy_v3_fl3(new_ml->bb->vec[2], +expx, +expy, -expz);  /* 2 */
1215                                                 copy_v3_fl3(new_ml->bb->vec[3], -expx, +expy, -expz);  /* 3 */
1216                                                 copy_v3_fl3(new_ml->bb->vec[4], -expx, -expy, +expz);  /* 4 */
1217                                                 copy_v3_fl3(new_ml->bb->vec[5], +expx, -expy, +expz);  /* 5 */
1218                                                 copy_v3_fl3(new_ml->bb->vec[6], +expx, +expy, +expz);  /* 6 */
1219                                                 copy_v3_fl3(new_ml->bb->vec[7], -expx, +expy, +expz);  /* 7 */
1220
1221                                                 /* transformation of Metalem bb */
1222                                                 for (i = 0; i < 8; i++)
1223                                                         mul_m4_v3((float(*)[4])new_ml->mat, new_ml->bb->vec[i]);
1224
1225                                                 /* find max and min of transformed bb */
1226                                                 INIT_MINMAX(tempmin, tempmax);
1227                                                 for (i = 0; i < 8; i++) {
1228                                                         DO_MINMAX(new_ml->bb->vec[i], tempmin, tempmax);
1229                                                 }
1230
1231                                                 /* set only point 0 and 6 - AABB of Metaelem */
1232                                                 copy_v3_v3(new_ml->bb->vec[0], tempmin);
1233                                                 copy_v3_v3(new_ml->bb->vec[6], tempmax);
1234
1235                                                 /* add new_ml to mainb[] */
1236                                                 if (UNLIKELY(process->totelem == process->mem)) {
1237                                                         process->mem = process->mem * 2 + 10;
1238                                                         process->mainb = MEM_reallocN(process->mainb, sizeof(MetaElem *) * process->mem);
1239                                                 }
1240                                                 process->mainb[process->totelem++] = new_ml;
1241                                         }
1242                                         ml = ml->next;
1243                                 }
1244                         }
1245                 }
1246         }
1247
1248         /* compute AABB of all Metaelems */
1249         if (process->totelem > 0) {
1250                 copy_v3_v3(process->allbb.min, process->mainb[0]->bb->vec[0]);
1251                 copy_v3_v3(process->allbb.max, process->mainb[0]->bb->vec[6]);
1252                 for (i = 1; i < process->totelem; i++)
1253                         make_box_union(process->mainb[i]->bb, &process->allbb, &process->allbb);
1254         }
1255 }
1256
1257 void BKE_mball_polygonize(EvaluationContext *eval_ctx, Scene *scene, Object *ob, ListBase *dispbase)
1258 {
1259         MetaBall *mb;
1260         DispList *dl;
1261         unsigned int a;
1262         PROCESS process = {0};
1263
1264         mb = ob->data;
1265
1266         process.thresh = mb->thresh;
1267
1268         if      (process.thresh < 0.001f) process.converge_res = 16;
1269         else if (process.thresh < 0.01f)  process.converge_res = 8;
1270         else if (process.thresh < 0.1f)   process.converge_res = 4;
1271         else                              process.converge_res = 2;
1272
1273         if ((eval_ctx->mode != DAG_EVAL_RENDER) && (mb->flag == MB_UPDATE_NEVER)) return;
1274         if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_FAST) return;
1275
1276         if (eval_ctx->mode == DAG_EVAL_RENDER) {
1277                 process.size = mb->rendersize;
1278         }
1279         else {
1280                 process.size = mb->wiresize;
1281                 if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_HALFRES) {
1282                         process.size *= 2.0f;
1283                 }
1284         }
1285
1286         process.delta = process.size * 0.001f;
1287
1288         process.pgn_elements = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "Metaball memarena");
1289
1290         /* initialize all mainb (MetaElems) */
1291         init_meta(eval_ctx, &process, scene, ob);
1292
1293         if (process.totelem > 0) {
1294                 build_bvh_spatial(&process, &process.metaball_bvh, 0, process.totelem, &process.allbb);
1295
1296                 /* don't polygonize metaballs with too high resolution (base mball to small)
1297                  * note: Eps was 0.0001f but this was giving problems for blood animation for durian, using 0.00001f */
1298                 if (ob->size[0] > 0.00001f * (process.allbb.max[0] - process.allbb.min[0]) ||
1299                     ob->size[1] > 0.00001f * (process.allbb.max[1] - process.allbb.min[1]) ||
1300                     ob->size[2] > 0.00001f * (process.allbb.max[2] - process.allbb.min[2]))
1301                 {
1302                         polygonize(&process);
1303
1304                         /* add resulting surface to displist */
1305                         if (process.curindex) {
1306                                 dl = MEM_callocN(sizeof(DispList), "mballdisp");
1307                                 BLI_addtail(dispbase, dl);
1308                                 dl->type = DL_INDEX4;
1309                                 dl->nr = (int)process.curvertex;
1310                                 dl->parts = (int)process.curindex;
1311
1312                                 dl->index = (int *)process.indices;
1313
1314                                 for (a = 0; a < process.curvertex; a++) {
1315                                         normalize_v3(process.no[a]);
1316                                 }
1317
1318                                 dl->verts = (float *)process.co;
1319                                 dl->nors = (float *)process.no;
1320                         }
1321                 }
1322         }
1323
1324         freepolygonize(&process);
1325 }