04024943310de6648b251c32fce15eb2530f8952
[blender-staging.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_path_util.h"
45 #include "BLI_math.h"
46 #include "BLI_utildefines.h"
47
48 #include "BKE_global.h"
49
50 #include "BKE_depsgraph.h"
51 #include "BKE_scene.h"
52 #include "BKE_displist.h"
53 #include "BKE_mball.h"
54 #include "BKE_mball_tessellate.h"  /* own include */
55
56 /* Data types */
57
58 typedef struct vertex {         /* surface vertex */
59         float co[3];  /* position and surface normal */
60         float no[3];
61 } VERTEX;
62
63 typedef struct vertices {       /* list of vertices in polygonization */
64         int count, max;             /* # vertices, max # allowed */
65         VERTEX *ptr;                /* dynamically allocated */
66 } VERTICES;
67
68 typedef struct corner {         /* corner of a cube */
69         int i, j, k;                /* (i, j, k) is index within lattice */
70         float co[3], value;       /* location and function value */
71         struct corner *next;
72 } CORNER;
73
74 typedef struct cube {           /* partitioning cell (cube) */
75         int i, j, k;                /* lattice location of cube */
76         CORNER *corners[8];         /* eight corners */
77 } CUBE;
78
79 typedef struct cubes {          /* linked list of cubes acting as stack */
80         CUBE cube;                  /* a single cube */
81         struct cubes *next;         /* remaining elements */
82 } CUBES;
83
84 typedef struct centerlist {     /* list of cube locations */
85         int i, j, k;                /* cube location */
86         struct centerlist *next;    /* remaining elements */
87 } CENTERLIST;
88
89 typedef struct edgelist {       /* list of edges */
90         int i1, j1, k1, i2, j2, k2; /* edge corner ids */
91         int vid;                    /* vertex id */
92         struct edgelist *next;      /* remaining elements */
93 } EDGELIST;
94
95 typedef struct intlist {        /* list of integers */
96         int i;                      /* an integer */
97         struct intlist *next;       /* remaining elements */
98 } INTLIST;
99
100 typedef struct intlists {       /* list of list of integers */
101         INTLIST *list;              /* a list of integers */
102         struct intlists *next;      /* remaining elements */
103 } INTLISTS;
104
105 /* dividing scene using octal tree makes polygonisation faster */
106 typedef struct ml_pointer {
107         struct ml_pointer *next, *prev;
108         struct MetaElem *ml;
109 } ml_pointer;
110
111 typedef struct octal_node {
112         struct octal_node *nodes[8];/* children of current node */
113         struct octal_node *parent;  /* parent of current node */
114         struct ListBase elems;      /* ListBase of MetaElem pointers (ml_pointer) */
115         float x_min, y_min, z_min;  /* 1st border point */
116         float x_max, y_max, z_max;  /* 7th border point */
117         float x, y, z;              /* center of node */
118         int pos, neg;               /* number of positive and negative MetaElements in the node */
119         int count;                  /* number of MetaElems, which belongs to the node */
120 } octal_node;
121
122 typedef struct octal_tree {
123         struct octal_node *first;   /* first node */
124         int pos, neg;               /* number of positive and negative MetaElements in the scene */
125         short depth;                /* number of scene subdivision */
126 } octal_tree;
127
128 struct pgn_elements {
129         struct pgn_elements *next, *prev;
130         char *data;
131 };
132
133 typedef struct process {        /* parameters, function, storage */
134         /* ** old G_mb contents ** */
135         float thresh;
136         int totelem;
137         MetaElem **mainb;
138         octal_tree *metaball_tree;
139
140         /* ** old process contents ** */
141
142         /* what happens here? floats, I think. */
143         /*  float (*function)(void);     */     /* implicit surface function */
144         float (*function)(struct process *, float, float, float);
145         float size, delta;          /* cube size, normal delta */
146         int bounds;                 /* cube range within lattice */
147         CUBES *cubes;               /* active cubes */
148         VERTICES vertices;          /* surface vertices */
149         CENTERLIST **centers;       /* cube center hash table */
150         CORNER **corners;           /* corner value hash table */
151         EDGELIST **edges;           /* edge and vertex id hash table */
152
153         /* Runtime things */
154         int *indices;
155         int totindex, curindex;
156
157         int pgn_offset;
158         struct pgn_elements *pgn_current;
159         ListBase pgn_list;
160 } PROCESS;
161
162 /* Forward declarations */
163 static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2, MetaBall *mb);
164 static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const int j, const int k);
165 static CORNER *setcorner(PROCESS *process, int i, int j, int k);
166 static void converge(PROCESS *process, const float p1[3], const float p2[3], float v1, float v2,
167                      float p[3], MetaBall *mb, int f);
168
169
170 /* ******************** ARITH ************************* */
171
172 /* BASED AT CODE (but mostly rewritten) :
173  * C code from the article
174  * "An Implicit Surface Polygonizer"
175  * by Jules Bloomenthal, jbloom@beauty.gmu.edu
176  * in "Graphics Gems IV", Academic Press, 1994
177  *
178  * Authored by Jules Bloomenthal, Xerox PARC.
179  * Copyright (c) Xerox Corporation, 1991.  All rights reserved.
180  * Permission is granted to reproduce, use and distribute this code for
181  * any and all purposes, provided that this notice appears in all copies. */
182
183 #define RES 12 /* # converge iterations    */
184
185 #define L   0  /* left direction:       -x, -i */
186 #define R   1  /* right direction:      +x, +i */
187 #define B   2  /* bottom direction: -y, -j */
188 #define T   3  /* top direction:        +y, +j */
189 #define N   4  /* near direction:       -z, -k */
190 #define F   5  /* far direction:        +z, +k */
191 #define LBN 0  /* left bottom near corner  */
192 #define LBF 1  /* left bottom far corner   */
193 #define LTN 2  /* left top near corner     */
194 #define LTF 3  /* left top far corner      */
195 #define RBN 4  /* right bottom near corner */
196 #define RBF 5  /* right bottom far corner  */
197 #define RTN 6  /* right top near corner    */
198 #define RTF 7  /* right top far corner     */
199
200 /* the LBN corner of cube (i, j, k), corresponds with location
201  * (i-0.5)*size, (j-0.5)*size, (k-0.5)*size) */
202
203 #define HASHBIT     (5)
204 #define HASHSIZE    (size_t)(1 << (3 * HASHBIT))   /*! < hash table size (32768) */
205
206 #define HASH(i, j, k) ((((( (i) & 31) << 5) | ( (j) & 31)) << 5) | ( (k) & 31) )
207
208 #define MB_BIT(i, bit) (((i) >> (bit)) & 1)
209 #define FLIP(i, bit) ((i) ^ 1 << (bit)) /* flip the given bit of i */
210
211
212 /* **************** POLYGONIZATION ************************ */
213
214 static void calc_mballco(MetaElem *ml, float vec[3])
215 {
216         if (ml->mat) {
217                 mul_m4_v3((float (*)[4])ml->mat, vec);
218         }
219 }
220
221 static float densfunc(MetaElem *ball, float x, float y, float z)
222 {
223         float dist2;
224         float dvec[3] = {x, y, z};
225
226         mul_m4_v3((float (*)[4])ball->imat, dvec);
227
228         switch (ball->type) {
229                 case MB_BALL:
230                         /* do nothing */
231                         break;
232                 case MB_TUBE:
233                         if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
234                         else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
235                         else                            dvec[0] = 0.0;
236                         break;
237                 case MB_PLANE:
238                         if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
239                         else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
240                         else                            dvec[0] = 0.0;
241                         if      (dvec[1] >  ball->expy) dvec[1] -= ball->expy;
242                         else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
243                         else                            dvec[1] = 0.0;
244                         break;
245                 case MB_ELIPSOID:
246                         dvec[0] /= ball->expx;
247                         dvec[1] /= ball->expy;
248                         dvec[2] /= ball->expz;
249                         break;
250                 case MB_CUBE:
251                         if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
252                         else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
253                         else                            dvec[0] = 0.0;
254
255                         if      (dvec[1] >  ball->expy) dvec[1] -= ball->expy;
256                         else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
257                         else                            dvec[1] = 0.0;
258
259                         if      (dvec[2] >  ball->expz) dvec[2] -= ball->expz;
260                         else if (dvec[2] < -ball->expz) dvec[2] += ball->expz;
261                         else                            dvec[2] = 0.0;
262                         break;
263
264                 /* *** deprecated, could be removed?, do-versioned at least *** */
265                 case MB_TUBEX:
266                         if      (dvec[0] >  ball->len) dvec[0] -= ball->len;
267                         else if (dvec[0] < -ball->len) dvec[0] += ball->len;
268                         else                           dvec[0] = 0.0;
269                         break;
270                 case MB_TUBEY:
271                         if      (dvec[1] >  ball->len) dvec[1] -= ball->len;
272                         else if (dvec[1] < -ball->len) dvec[1] += ball->len;
273                         else                           dvec[1] = 0.0;
274                         break;
275                 case MB_TUBEZ:
276                         if      (dvec[2] >  ball->len) dvec[2] -= ball->len;
277                         else if (dvec[2] < -ball->len) dvec[2] += ball->len;
278                         else                           dvec[2] = 0.0;
279                         break;
280                 /* *** end deprecated *** */
281         }
282
283         dist2 = 1.0f - (len_squared_v3(dvec) / ball->rad2);
284
285         if ((ball->flag & MB_NEGATIVE) == 0) {
286                 return (dist2 < 0.0f) ? -0.5f : (ball->s * dist2 * dist2 * dist2) - 0.5f;
287         }
288         else {
289                 return (dist2 < 0.0f) ? 0.5f : 0.5f - (ball->s * dist2 * dist2 * dist2);
290         }
291 }
292
293 static octal_node *find_metaball_octal_node(octal_node *node, float x, float y, float z, short depth)
294 {
295         if (!depth) return node;
296         
297         if (z < node->z) {
298                 if (y < node->y) {
299                         if (x < node->x) {
300                                 if (node->nodes[0])
301                                         return find_metaball_octal_node(node->nodes[0], x, y, z, depth--);
302                                 else
303                                         return node;
304                         }
305                         else {
306                                 if (node->nodes[1])
307                                         return find_metaball_octal_node(node->nodes[1], x, y, z, depth--);
308                                 else
309                                         return node;
310                         }
311                 }
312                 else {
313                         if (x < node->x) {
314                                 if (node->nodes[3])
315                                         return find_metaball_octal_node(node->nodes[3], x, y, z, depth--);
316                                 else
317                                         return node;
318                         }
319                         else {
320                                 if (node->nodes[2])
321                                         return find_metaball_octal_node(node->nodes[2], x, y, z, depth--);
322                                 else
323                                         return node;
324                         }
325                 }
326         }
327         else {
328                 if (y < node->y) {
329                         if (x < node->x) {
330                                 if (node->nodes[4])
331                                         return find_metaball_octal_node(node->nodes[4], x, y, z, depth--);
332                                 else
333                                         return node;
334                         }
335                         else {
336                                 if (node->nodes[5])
337                                         return find_metaball_octal_node(node->nodes[5], x, y, z, depth--);
338                                 else
339                                         return node;
340                         }
341                 }
342                 else {
343                         if (x < node->x) {
344                                 if (node->nodes[7])
345                                         return find_metaball_octal_node(node->nodes[7], x, y, z, depth--);
346                                 else
347                                         return node;
348                         }
349                         else {
350                                 if (node->nodes[6])
351                                         return find_metaball_octal_node(node->nodes[6], x, y, z, depth--);
352                                 else
353                                         return node;
354                         }
355                 }
356         }
357         
358         /* all cases accounted for */
359         BLI_assert(0);
360 }
361
362 static float metaball(PROCESS *process, float x, float y, float z)
363 /*  float x, y, z; */
364 {
365         octal_tree *metaball_tree = process->metaball_tree;
366         struct octal_node *node;
367         struct ml_pointer *ml_p;
368         float dens = 0;
369         int a;
370         
371         if (process->totelem > 1) {
372                 node = find_metaball_octal_node(metaball_tree->first, x, y, z, metaball_tree->depth);
373                 if (node) {
374                         for (ml_p = node->elems.first; ml_p; ml_p = ml_p->next) {
375                                 dens += densfunc(ml_p->ml, x, y, z);
376                         }
377
378                         dens += -0.5f * (metaball_tree->pos - node->pos);
379                         dens +=  0.5f * (metaball_tree->neg - node->neg);
380                 }
381                 else {
382                         for (a = 0; a < process->totelem; a++) {
383                                 dens += densfunc(process->mainb[a], x, y, z);
384                         }
385                 }
386         }
387         else {
388                 dens += densfunc(process->mainb[0], x, y, z);
389         }
390
391         return process->thresh - dens;
392 }
393
394 /* ******************************************** */
395
396 static void accum_mballfaces(PROCESS *process, int i1, int i2, int i3, int i4)
397 {
398         int *newi, *cur;
399         /* static int i = 0; I would like to delete altogether, but I don't dare to, yet */
400
401         if (process->totindex == process->curindex) {
402                 process->totindex += 256;
403                 newi = MEM_mallocN(4 * sizeof(int) * process->totindex, "vertindex");
404                 
405                 if (process->indices) {
406                         memcpy(newi, process->indices, 4 * sizeof(int) * (process->totindex - 256));
407                         MEM_freeN(process->indices);
408                 }
409                 process->indices = newi;
410         }
411         
412         cur = process->indices + 4 * process->curindex;
413
414         /* displists now support array drawing, we treat tri's as fake quad */
415         
416         cur[0] = i1;
417         cur[1] = i2;
418         cur[2] = i3;
419         if (i4 == 0)
420                 cur[3] = i3;
421         else 
422                 cur[3] = i4;
423         
424         process->curindex++;
425
426 }
427
428 /* ******************* MEMORY MANAGEMENT *********************** */
429 static void *new_pgn_element(PROCESS *process, int size)
430 {
431         /* during polygonize 1000s of elements are allocated
432          * and never freed in between. Freeing only done at the end.
433          */
434         int blocksize = 16384;
435         void *adr;
436         
437         if (size > 10000 || size == 0) {
438                 printf("incorrect use of new_pgn_element\n");
439         }
440         else if (size == -1) {
441                 struct pgn_elements *cur = process->pgn_list.first;
442                 while (cur) {
443                         MEM_freeN(cur->data);
444                         cur = cur->next;
445                 }
446                 BLI_freelistN(&process->pgn_list);
447                 
448                 return NULL;
449         }
450         
451         size = 4 * ( (size + 3) / 4);
452         
453         if (process->pgn_current) {
454                 if (size + process->pgn_offset < blocksize) {
455                         adr = (void *) (process->pgn_current->data + process->pgn_offset);
456                         process->pgn_offset += size;
457                         return adr;
458                 }
459         }
460         
461         process->pgn_current = MEM_callocN(sizeof(struct pgn_elements), "newpgn");
462         process->pgn_current->data = MEM_callocN(blocksize, "newpgn");
463         BLI_addtail(&process->pgn_list, process->pgn_current);
464         
465         process->pgn_offset = size;
466         return process->pgn_current->data;
467 }
468
469 static void freepolygonize(PROCESS *process)
470 {
471         MEM_freeN(process->corners);
472         MEM_freeN(process->edges);
473         MEM_freeN(process->centers);
474
475         new_pgn_element(process, -1);
476
477         if (process->vertices.ptr) {
478                 MEM_freeN(process->vertices.ptr);
479         }
480 }
481
482 /**** Cubical Polygonization (optional) ****/
483
484 #define LB  0  /* left bottom edge      */
485 #define LT  1  /* left top edge */
486 #define LN  2  /* left near edge        */
487 #define LF  3  /* left far edge */
488 #define RB  4  /* right bottom edge */
489 #define RT  5  /* right top edge        */
490 #define RN  6  /* right near edge       */
491 #define RF  7  /* right far edge        */
492 #define BN  8  /* bottom near edge      */
493 #define BF  9  /* bottom far edge       */
494 #define TN  10 /* top near edge */
495 #define TF  11 /* top far edge  */
496
497 static INTLISTS *cubetable[256];
498
499 /* edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
500 static int corner1[12] = {
501         LBN, LTN, LBN, LBF, RBN, RTN, RBN, RBF, LBN, LBF, LTN, LTF
502 };
503 static int corner2[12] = {
504         LBF, LTF, LTN, LTF, RBF, RTF, RTN, RTF, RBN, RBF, RTN, RTF
505 };
506 static int leftface[12] = {
507         B,  L,  L,  F,  R,  T,  N,  R,  N,  B,  T,  F
508 };
509 /* face on left when going corner1 to corner2 */
510 static int rightface[12] = {
511         L,  T,  N,  L,  B,  R,  R,  F,  B,  F,  N,  T
512 };
513 /* face on right when going corner1 to corner2 */
514
515
516 /* docube: triangulate the cube directly, without decomposition */
517
518 static void docube(PROCESS *process, CUBE *cube, MetaBall *mb)
519 {
520         INTLISTS *polys;
521         CORNER *c1, *c2;
522         int i, index = 0, count, indexar[8];
523         
524         for (i = 0; i < 8; i++) {
525                 if (cube->corners[i]->value > 0.0f) {
526                         index += (1 << i);
527                 }
528         }
529         
530         for (polys = cubetable[index]; polys; polys = polys->next) {
531                 INTLIST *edges;
532                 
533                 count = 0;
534                 
535                 for (edges = polys->list; edges; edges = edges->next) {
536                         c1 = cube->corners[corner1[edges->i]];
537                         c2 = cube->corners[corner2[edges->i]];
538                         
539                         indexar[count] = vertid(process, c1, c2, mb);
540                         count++;
541                 }
542                 if (count > 2) {
543                         switch (count) {
544                                 case 3:
545                                         accum_mballfaces(process, indexar[2], indexar[1], indexar[0], 0);
546                                         break;
547                                 case 4:
548                                         if (indexar[0] == 0) accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
549                                         else accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
550                                         break;
551                                 case 5:
552                                         if (indexar[0] == 0) accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
553                                         else accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
554                                 
555                                         accum_mballfaces(process, indexar[4], indexar[3], indexar[0], 0);
556                                         break;
557                                 case 6:
558                                         if (indexar[0] == 0) {
559                                                 accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
560                                                 accum_mballfaces(process, indexar[0], indexar[5], indexar[4], indexar[3]);
561                                         }
562                                         else {
563                                                 accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
564                                                 accum_mballfaces(process, indexar[5], indexar[4], indexar[3], indexar[0]);
565                                         }
566                                         break;
567                                 case 7:
568                                         if (indexar[0] == 0) {
569                                                 accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
570                                                 accum_mballfaces(process, indexar[0], indexar[5], indexar[4], indexar[3]);
571                                         }
572                                         else {
573                                                 accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
574                                                 accum_mballfaces(process, indexar[5], indexar[4], indexar[3], indexar[0]);
575                                         }
576                                 
577                                         accum_mballfaces(process, indexar[6], indexar[5], indexar[0], 0);
578
579                                         break;
580                         }
581                 }
582         }
583 }
584
585
586 /* testface: given cube at lattice (i, j, k), and four corners of face,
587  * if surface crosses face, compute other four corners of adjacent cube
588  * and add new cube to cube stack */
589
590 static void testface(PROCESS *process, int i, int j, int k, CUBE *old, int bit, int c1, int c2, int c3, int c4)
591 {
592         CUBE newc;
593         CUBES *oldcubes = process->cubes;
594         CORNER *corn1, *corn2, *corn3, *corn4;
595         int n, pos;
596
597         corn1 = old->corners[c1];
598         corn2 = old->corners[c2];
599         corn3 = old->corners[c3];
600         corn4 = old->corners[c4];
601         
602         pos = corn1->value > 0.0f ? 1 : 0;
603
604         /* test if no surface crossing */
605         if ( (corn2->value > 0) == pos && (corn3->value > 0) == pos && (corn4->value > 0) == pos) return;
606         /* test if cube out of bounds */
607         /*if ( abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) return;*/
608         /* test if already visited (always as last) */
609         if (setcenter(process, process->centers, i, j, k)) {
610                 return;
611         }
612
613         /* create new cube and add cube to top of stack: */
614         process->cubes = (CUBES *) new_pgn_element(process, sizeof(CUBES));
615         process->cubes->next = oldcubes;
616         
617         newc.i = i;
618         newc.j = j;
619         newc.k = k;
620         for (n = 0; n < 8; n++) newc.corners[n] = NULL;
621         
622         newc.corners[FLIP(c1, bit)] = corn1;
623         newc.corners[FLIP(c2, bit)] = corn2;
624         newc.corners[FLIP(c3, bit)] = corn3;
625         newc.corners[FLIP(c4, bit)] = corn4;
626
627         if (newc.corners[0] == NULL) newc.corners[0] = setcorner(process, i, j, k);
628         if (newc.corners[1] == NULL) newc.corners[1] = setcorner(process, i, j, k + 1);
629         if (newc.corners[2] == NULL) newc.corners[2] = setcorner(process, i, j + 1, k);
630         if (newc.corners[3] == NULL) newc.corners[3] = setcorner(process, i, j + 1, k + 1);
631         if (newc.corners[4] == NULL) newc.corners[4] = setcorner(process, i + 1, j, k);
632         if (newc.corners[5] == NULL) newc.corners[5] = setcorner(process, i + 1, j, k + 1);
633         if (newc.corners[6] == NULL) newc.corners[6] = setcorner(process, i + 1, j + 1, k);
634         if (newc.corners[7] == NULL) newc.corners[7] = setcorner(process, i + 1, j + 1, k + 1);
635
636         process->cubes->cube = newc;
637 }
638
639 /* setcorner: return corner with the given lattice location
640  * set (and cache) its function value */
641
642 static CORNER *setcorner(PROCESS *process, int i, int j, int k)
643 {
644         /* for speed, do corner value caching here */
645         CORNER *c;
646         int index;
647
648         /* does corner exist? */
649         index = HASH(i, j, k);
650         c = process->corners[index];
651         
652         for (; c != NULL; c = c->next) {
653                 if (c->i == i && c->j == j && c->k == k) {
654                         return c;
655                 }
656         }
657
658         c = (CORNER *) new_pgn_element(process, sizeof(CORNER));
659
660         c->i = i; 
661         c->co[0] = ((float)i - 0.5f) * process->size;
662         c->j = j; 
663         c->co[1] = ((float)j - 0.5f) * process->size;
664         c->k = k; 
665         c->co[2] = ((float)k - 0.5f) * process->size;
666         c->value = process->function(process, c->co[0], c->co[1], c->co[2]);
667         
668         c->next = process->corners[index];
669         process->corners[index] = c;
670         
671         return c;
672 }
673
674
675 /* nextcwedge: return next clockwise edge from given edge around given face */
676
677 static int nextcwedge(int edge, int face)
678 {
679         switch (edge) {
680                 case LB:
681                         return (face == L) ? LF : BN;
682                 case LT:
683                         return (face == L) ? LN : TF;
684                 case LN:
685                         return (face == L) ? LB : TN;
686                 case LF:
687                         return (face == L) ? LT : BF;
688                 case RB:
689                         return (face == R) ? RN : BF;
690                 case RT:
691                         return (face == R) ? RF : TN;
692                 case RN:
693                         return (face == R) ? RT : BN;
694                 case RF:
695                         return (face == R) ? RB : TF;
696                 case BN:
697                         return (face == B) ? RB : LN;
698                 case BF:
699                         return (face == B) ? LB : RF;
700                 case TN:
701                         return (face == T) ? LT : RN;
702                 case TF:
703                         return (face == T) ? RT : LF;
704         }
705         return 0;
706 }
707
708
709 /* otherface: return face adjoining edge that is not the given face */
710
711 static int otherface(int edge, int face)
712 {
713         int other = leftface[edge];
714         return face == other ? rightface[edge] : other;
715 }
716
717
718 /* makecubetable: create the 256 entry table for cubical polygonization */
719
720 static void makecubetable(void)
721 {
722         static bool is_done = false;
723         int i, e, c, done[12], pos[8];
724
725         if (is_done) return;
726         is_done = true;
727
728         for (i = 0; i < 256; i++) {
729                 for (e = 0; e < 12; e++) done[e] = 0;
730                 for (c = 0; c < 8; c++) pos[c] = MB_BIT(i, c);
731                 for (e = 0; e < 12; e++)
732                         if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
733                                 INTLIST *ints = NULL;
734                                 INTLISTS *lists = (INTLISTS *) MEM_callocN(sizeof(INTLISTS), "mball_intlist");
735                                 int start = e, edge = e;
736                                 
737                                 /* get face that is to right of edge from pos to neg corner: */
738                                 int face = pos[corner1[e]] ? rightface[e] : leftface[e];
739                                 
740                                 while (1) {
741                                         edge = nextcwedge(edge, face);
742                                         done[edge] = 1;
743                                         if (pos[corner1[edge]] != pos[corner2[edge]]) {
744                                                 INTLIST *tmp = ints;
745                                                 
746                                                 ints = (INTLIST *) MEM_callocN(sizeof(INTLIST), "mball_intlist");
747                                                 ints->i = edge;
748                                                 ints->next = tmp; /* add edge to head of list */
749                                                 
750                                                 if (edge == start) break;
751                                                 face = otherface(edge, face);
752                                         }
753                                 }
754                                 lists->list = ints; /* add ints to head of table entry */
755                                 lists->next = cubetable[i];
756                                 cubetable[i] = lists;
757                         }
758         }
759 }
760
761 void BKE_mball_cubeTable_free(void)
762 {
763         int i;
764         INTLISTS *lists, *nlists;
765         INTLIST *ints, *nints;
766
767         for (i = 0; i < 256; i++) {
768                 lists = cubetable[i];
769                 while (lists) {
770                         nlists = lists->next;
771                         
772                         ints = lists->list;
773                         while (ints) {
774                                 nints = ints->next;
775                                 MEM_freeN(ints);
776                                 ints = nints;
777                         }
778                         
779                         MEM_freeN(lists);
780                         lists = nlists;
781                 }
782                 cubetable[i] = NULL;
783         }
784 }
785
786 /**** Storage ****/
787
788 /* setcenter: set (i, j, k) entry of table[]
789  * return 1 if already set; otherwise, set and return 0 */
790
791 static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const int j, const int k)
792 {
793         int index;
794         CENTERLIST *newc, *l, *q;
795
796         index = HASH(i, j, k);
797         q = table[index];
798
799         for (l = q; l != NULL; l = l->next) {
800                 if (l->i == i && l->j == j && l->k == k) return 1;
801         }
802         
803         newc = (CENTERLIST *) new_pgn_element(process, sizeof(CENTERLIST));
804         newc->i = i; 
805         newc->j = j; 
806         newc->k = k; 
807         newc->next = q;
808         table[index] = newc;
809         
810         return 0;
811 }
812
813
814 /* setedge: set vertex id for edge */
815
816 static void setedge(PROCESS *process,
817                     EDGELIST *table[],
818                     int i1, int j1,
819                     int k1, int i2,
820                     int j2, int k2,
821                     int vid)
822 {
823         unsigned int index;
824         EDGELIST *newe;
825         
826         if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
827                 int t = i1;
828                 i1 = i2;
829                 i2 = t;
830                 t = j1;
831                 j1 = j2;
832                 j2 = t;
833                 t = k1;
834                 k1 = k2;
835                 k2 = t;
836         }
837         index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
838         newe = (EDGELIST *) new_pgn_element(process, sizeof(EDGELIST));
839         newe->i1 = i1; 
840         newe->j1 = j1; 
841         newe->k1 = k1;
842         newe->i2 = i2; 
843         newe->j2 = j2; 
844         newe->k2 = k2;
845         newe->vid = vid;
846         newe->next = table[index];
847         table[index] = newe;
848 }
849
850
851 /* getedge: return vertex id for edge; return -1 if not set */
852
853 static int getedge(EDGELIST *table[],
854                    int i1, int j1, int k1,
855                    int i2, int j2, int k2)
856 {
857         EDGELIST *q;
858         
859         if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
860                 int t = i1;
861                 i1 = i2;
862                 i2 = t;
863                 t = j1;
864                 j1 = j2;
865                 j2 = t;
866                 t = k1;
867                 k1 = k2;
868                 k2 = t;
869         }
870         q = table[HASH(i1, j1, k1) + HASH(i2, j2, k2)];
871         for (; q != NULL; q = q->next) {
872                 if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 &&
873                     q->i2 == i2 && q->j2 == j2 && q->k2 == k2)
874                 {
875                         return q->vid;
876                 }
877         }
878         return -1;
879 }
880
881
882 /**** Vertices ****/
883
884 #undef R
885
886
887
888 /* vertid: return index for vertex on edge:
889  * c1->value and c2->value are presumed of different sign
890  * return saved index if any; else compute vertex and save */
891
892 /* addtovertices: add v to sequence of vertices */
893
894 static void addtovertices(VERTICES *vertices, VERTEX v)
895 {
896         if (vertices->count == vertices->max) {
897                 int i;
898                 VERTEX *newv;
899                 vertices->max = vertices->count == 0 ? 10 : 2 * vertices->count;
900                 newv = (VERTEX *) MEM_callocN(vertices->max * sizeof(VERTEX), "addtovertices");
901                 
902                 for (i = 0; i < vertices->count; i++) newv[i] = vertices->ptr[i];
903                 
904                 if (vertices->ptr != NULL) MEM_freeN(vertices->ptr);
905                 vertices->ptr = newv;
906         }
907         vertices->ptr[vertices->count++] = v;
908 }
909
910 /* vnormal: compute unit length surface normal at point */
911
912 static void vnormal(PROCESS *process, const float point[3], float r_no[3])
913 {
914         const float delta = 0.2f * process->delta;
915         const float f = process->function(process, point[0], point[1], point[2]);
916
917         r_no[0] = process->function(process, point[0] + delta, point[1], point[2]) - f;
918         r_no[1] = process->function(process, point[0], point[1] + delta, point[2]) - f;
919         r_no[2] = process->function(process, point[0], point[1], point[2] + delta) - f;
920
921 #if 1
922         normalize_v3(r_no);
923 #else
924         f = normalize_v3(r_no);
925         
926         if (0) {
927                 float tvec[3];
928                 
929                 delta *= 2.0f;
930                 
931                 f = process->function(process, point[0], point[1], point[2]);
932         
933                 tvec[0] = process->function(process, point[0] + delta, point[1], point[2]) - f;
934                 tvec[1] = process->function(process, point[0], point[1] + delta, point[2]) - f;
935                 tvec[2] = process->function(process, point[0], point[1], point[2] + delta) - f;
936         
937                 if (normalize_v3(tvec) != 0.0f) {
938                         add_v3_v3(r_no, tvec);
939                         normalize_v3(r_no);
940                 }
941         }
942 #endif
943 }
944
945
946 static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2, MetaBall *mb)
947 {
948         VERTEX v;
949         int vid = getedge(process->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
950
951         if (vid != -1) {
952                 return vid;  /* previously computed */
953         }
954
955         converge(process, c1->co, c2->co, c1->value, c2->value, v.co, mb, 1); /* position */
956         vnormal(process, v.co, v.no);
957
958         addtovertices(&process->vertices, v);            /* save vertex */
959         vid = process->vertices.count - 1;
960         setedge(process, process->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
961         
962         return vid;
963 }
964
965
966 /* converge: from two points of differing sign, converge to zero crossing */
967 /* watch it: p1 and p2 are used to calculate */
968 static void converge(PROCESS *process, const float p1[3], const float p2[3], float v1, float v2,
969                      float p[3], MetaBall *mb, int f)
970 {
971         int i = 0;
972         float pos[3], neg[3];
973         float positive = 0.0f, negative = 0.0f;
974         float dvec[3];
975         
976         if (v1 < 0) {
977                 copy_v3_v3(pos, p2);
978                 copy_v3_v3(neg, p1);
979                 positive = v2;
980                 negative = v1;
981         }
982         else {
983                 copy_v3_v3(pos, p1);
984                 copy_v3_v3(neg, p2);
985                 positive = v1;
986                 negative = v2;
987         }
988
989         sub_v3_v3v3(dvec, pos, neg);
990
991 /* Approximation by linear interpolation is faster then binary subdivision,
992  * but it results sometimes (mb->thresh < 0.2) into the strange results */
993         if ((mb->thresh > 0.2f) && (f == 1)) {
994                 if ((dvec[1] == 0.0f) && (dvec[2] == 0.0f)) {
995                         p[0] = neg[0] - negative * dvec[0] / (positive - negative);
996                         p[1] = neg[1];
997                         p[2] = neg[2];
998                         return;
999                 }
1000                 if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
1001                         p[0] = neg[0];
1002                         p[1] = neg[1] - negative * dvec[1] / (positive - negative);
1003                         p[2] = neg[2];
1004                         return;
1005                 }
1006                 if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
1007                         p[0] = neg[0];
1008                         p[1] = neg[1];
1009                         p[2] = neg[2] - negative * dvec[2] / (positive - negative);
1010                         return;
1011                 }
1012         }
1013
1014         if ((dvec[1] == 0.0f) && (dvec[2] == 0.0f)) {
1015                 p[1] = neg[1];
1016                 p[2] = neg[2];
1017                 while (1) {
1018                         if (i++ == RES) return;
1019                         p[0] = 0.5f * (pos[0] + neg[0]);
1020                         if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[0] = p[0];
1021                         else                                                       neg[0] = p[0];
1022                 }
1023         }
1024
1025         if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
1026                 p[0] = neg[0];
1027                 p[2] = neg[2];
1028                 while (1) {
1029                         if (i++ == RES) return;
1030                         p[1] = 0.5f * (pos[1] + neg[1]);
1031                         if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[1] = p[1];
1032                         else                                                       neg[1] = p[1];
1033                 }
1034         }
1035
1036         if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
1037                 p[0] = neg[0];
1038                 p[1] = neg[1];
1039                 while (1) {
1040                         if (i++ == RES) return;
1041                         p[2] = 0.5f * (pos[2] + neg[2]);
1042                         if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[2] = p[2];
1043                         else                                                       neg[2] = p[2];
1044                 }
1045         }
1046
1047         /* This is necessary to find start point */
1048         while (1) {
1049                 mid_v3_v3v3(&p[0], pos, neg);
1050
1051                 if (i++ == RES) {
1052                         return;
1053                 }
1054
1055                 if ((process->function(process, p[0], p[1], p[2])) > 0.0f) {
1056                         copy_v3_v3(pos, &p[0]);
1057                 }
1058                 else {
1059                         copy_v3_v3(neg, &p[0]);
1060                 }
1061         }
1062 }
1063
1064 /* ************************************** */
1065 static void add_cube(PROCESS *process, int i, int j, int k, int count)
1066 {
1067         CUBES *ncube;
1068         int n;
1069         int a, b, c;
1070
1071         /* hmmm, not only one, but eight cube will be added on the stack 
1072          * ... */
1073         for (a = i - 1; a < i + count; a++)
1074                 for (b = j - 1; b < j + count; b++)
1075                         for (c = k - 1; c < k + count; c++) {
1076                                 /* test if cube has been found before */
1077                                 if (setcenter(process, process->centers, a, b, c) == 0) {
1078                                         /* push cube on stack: */
1079                                         ncube = (CUBES *) new_pgn_element(process, sizeof(CUBES));
1080                                         ncube->next = process->cubes;
1081                                         process->cubes = ncube;
1082
1083                                         ncube->cube.i = a;
1084                                         ncube->cube.j = b;
1085                                         ncube->cube.k = c;
1086
1087                                         /* set corners of initial cube: */
1088                                         for (n = 0; n < 8; n++)
1089                                                 ncube->cube.corners[n] = setcorner(process, a + MB_BIT(n, 2), b + MB_BIT(n, 1), c + MB_BIT(n, 0));
1090                                 }
1091                         }
1092 }
1093
1094
1095 static void find_first_points(PROCESS *process, MetaBall *mb, int a)
1096 {
1097         MetaElem *ml;
1098         float f;
1099
1100         ml = process->mainb[a];
1101         f = 1.0f - (mb->thresh / ml->s);
1102
1103         /* Skip, when Stiffness of MetaElement is too small ... MetaElement can't be
1104          * visible alone ... but still can influence others MetaElements :-) */
1105         if (f > 0.0f) {
1106                 float IN[3] = {0.0f}, OUT[3] = {0.0f}, in[3] = {0.0f}, out[3];
1107                 int i, j, k, c_i, c_j, c_k;
1108                 int index[3] = {1, 0, -1};
1109                 float in_v /*, out_v*/;
1110                 float workp[3];
1111                 float dvec[3];
1112                 float tmp_v, workp_v, max_len_sq, nx, ny, nz, max_dim;
1113
1114                 calc_mballco(ml, in);
1115                 in_v = process->function(process, in[0], in[1], in[2]);
1116
1117                 for (i = 0; i < 3; i++) {
1118                         switch (ml->type) {
1119                                 case MB_BALL:
1120                                         OUT[0] = out[0] = IN[0] + index[i] * ml->rad;
1121                                         break;
1122                                 case MB_TUBE:
1123                                 case MB_PLANE:
1124                                 case MB_ELIPSOID:
1125                                 case MB_CUBE:
1126                                         OUT[0] = out[0] = IN[0] + index[i] * (ml->expx + ml->rad);
1127                                         break;
1128                         }
1129
1130                         for (j = 0; j < 3; j++) {
1131                                 switch (ml->type) {
1132                                         case MB_BALL:
1133                                                 OUT[1] = out[1] = IN[1] + index[j] * ml->rad;
1134                                                 break;
1135                                         case MB_TUBE:
1136                                         case MB_PLANE:
1137                                         case MB_ELIPSOID:
1138                                         case MB_CUBE:
1139                                                 OUT[1] = out[1] = IN[1] + index[j] * (ml->expy + ml->rad);
1140                                                 break;
1141                                 }
1142                         
1143                                 for (k = 0; k < 3; k++) {
1144                                         out[0] = OUT[0];
1145                                         out[1] = OUT[1];
1146                                         switch (ml->type) {
1147                                                 case MB_BALL:
1148                                                 case MB_TUBE:
1149                                                 case MB_PLANE:
1150                                                         out[2] = IN[2] + index[k] * ml->rad;
1151                                                         break;
1152                                                 case MB_ELIPSOID:
1153                                                 case MB_CUBE:
1154                                                         out[2] = IN[2] + index[k] * (ml->expz + ml->rad);
1155                                                         break;
1156                                         }
1157
1158                                         calc_mballco(ml, out);
1159
1160                                         /*out_v = process->function(out[0], out[1], out[2]);*/ /*UNUSED*/
1161
1162                                         /* find "first points" on Implicit Surface of MetaElemnt ml */
1163                                         copy_v3_v3(workp, in);
1164                                         workp_v = in_v;
1165                                         max_len_sq = len_squared_v3v3(out, in);
1166
1167                                         nx = fabsf((out[0] - in[0]) / process->size);
1168                                         ny = fabsf((out[1] - in[1]) / process->size);
1169                                         nz = fabsf((out[2] - in[2]) / process->size);
1170                                         
1171                                         max_dim = max_fff(nx, ny, nz);
1172                                         if (max_dim != 0.0f) {
1173                                                 float len_sq = 0.0f;
1174
1175                                                 dvec[0] = (out[0] - in[0]) / max_dim;
1176                                                 dvec[1] = (out[1] - in[1]) / max_dim;
1177                                                 dvec[2] = (out[2] - in[2]) / max_dim;
1178
1179                                                 while (len_sq <= max_len_sq) {
1180                                                         add_v3_v3(workp, dvec);
1181
1182                                                         /* compute value of implicite function */
1183                                                         tmp_v = process->function(process, workp[0], workp[1], workp[2]);
1184                                                         /* add cube to the stack, when value of implicite function crosses zero value */
1185                                                         if ((tmp_v < 0.0f && workp_v >= 0.0f) || (tmp_v > 0.0f && workp_v <= 0.0f)) {
1186
1187                                                                 /* indexes of CUBE, which includes "first point" */
1188                                                                 c_i = (int)floor(workp[0] / process->size);
1189                                                                 c_j = (int)floor(workp[1] / process->size);
1190                                                                 c_k = (int)floor(workp[2] / process->size);
1191                                                                 
1192                                                                 /* add CUBE (with indexes c_i, c_j, c_k) to the stack,
1193                                                                  * this cube includes found point of Implicit Surface */
1194                                                                 if ((ml->flag & MB_NEGATIVE) == 0) {
1195                                                                         add_cube(process, c_i, c_j, c_k, 1);
1196                                                                 }
1197                                                                 else {
1198                                                                         add_cube(process, c_i, c_j, c_k, 2);
1199                                                                 }
1200                                                         }
1201                                                         len_sq = len_squared_v3v3(workp, in);
1202                                                         workp_v = tmp_v;
1203
1204                                                 }
1205                                         }
1206                                 }
1207                         }
1208                 }
1209         }
1210 }
1211
1212 static void polygonize(PROCESS *process, MetaBall *mb)
1213 {
1214         CUBE c;
1215         int a;
1216
1217         process->vertices.count = process->vertices.max = 0;
1218         process->vertices.ptr = NULL;
1219
1220         /* allocate hash tables and build cube polygon table: */
1221         process->centers = MEM_callocN(HASHSIZE * sizeof(CENTERLIST *), "mbproc->centers");
1222         process->corners = MEM_callocN(HASHSIZE * sizeof(CORNER *), "mbproc->corners");
1223         process->edges = MEM_callocN(2 * HASHSIZE * sizeof(EDGELIST *), "mbproc->edges");
1224         makecubetable();
1225
1226         for (a = 0; a < process->totelem; a++) {
1227
1228                 /* try to find 8 points on the surface for each MetaElem */
1229                 find_first_points(process, mb, a);
1230         }
1231
1232         /* polygonize all MetaElems of current MetaBall */
1233         while (process->cubes != NULL) { /* process active cubes till none left */
1234                 c = process->cubes->cube;
1235
1236                 /* polygonize the cube directly: */
1237                 docube(process, &c, mb);
1238                 
1239                 /* pop current cube from stack */
1240                 process->cubes = process->cubes->next;
1241                 
1242                 /* test six face directions, maybe add to stack: */
1243                 testface(process, c.i - 1, c.j, c.k, &c, 2, LBN, LBF, LTN, LTF);
1244                 testface(process, c.i + 1, c.j, c.k, &c, 2, RBN, RBF, RTN, RTF);
1245                 testface(process, c.i, c.j - 1, c.k, &c, 1, LBN, LBF, RBN, RBF);
1246                 testface(process, c.i, c.j + 1, c.k, &c, 1, LTN, LTF, RTN, RTF);
1247                 testface(process, c.i, c.j, c.k - 1, &c, 0, LBN, LTN, RBN, RTN);
1248                 testface(process, c.i, c.j, c.k + 1, &c, 0, LBF, LTF, RBF, RTF);
1249         }
1250 }
1251
1252 static float init_meta(EvaluationContext *eval_ctx, PROCESS *process, Scene *scene, Object *ob)    /* return totsize */
1253 {
1254         Scene *sce_iter = scene;
1255         Base *base;
1256         Object *bob;
1257         MetaBall *mb;
1258         MetaElem *ml;
1259         float size, totsize, obinv[4][4], obmat[4][4], vec[3];
1260         //float max = 0.0f;
1261         int a, obnr, zero_size = 0;
1262         char obname[MAX_ID_NAME];
1263         SceneBaseIter iter;
1264
1265         copy_m4_m4(obmat, ob->obmat);   /* to cope with duplicators from BKE_scene_base_iter_next */
1266         invert_m4_m4(obinv, ob->obmat);
1267         a = 0;
1268         
1269         BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
1270         
1271         /* make main array */
1272         BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 0, NULL, NULL);
1273         while (BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 1, &base, &bob)) {
1274
1275                 if (bob->type == OB_MBALL) {
1276                         zero_size = 0;
1277                         ml = NULL;
1278
1279                         if (bob == ob && (base->flag & OB_FROMDUPLI) == 0) {
1280                                 mb = ob->data;
1281         
1282                                 if (mb->editelems) ml = mb->editelems->first;
1283                                 else ml = mb->elems.first;
1284                         }
1285                         else {
1286                                 char name[MAX_ID_NAME];
1287                                 int nr;
1288                                 
1289                                 BLI_split_name_num(name, &nr, bob->id.name + 2, '.');
1290                                 if (STREQ(obname, name)) {
1291                                         mb = bob->data;
1292                                         
1293                                         if (mb->editelems) ml = mb->editelems->first;
1294                                         else ml = mb->elems.first;
1295                                 }
1296                         }
1297
1298                         /* when metaball object has zero scale, then MetaElem to this MetaBall
1299                          * will not be put to mainb array */
1300                         if (has_zero_axis_m4(bob->obmat)) {
1301                                 zero_size = 1;
1302                         }
1303                         else if (bob->parent) {
1304                                 struct Object *pob = bob->parent;
1305                                 while (pob) {
1306                                         if (has_zero_axis_m4(pob->obmat)) {
1307                                                 zero_size = 1;
1308                                                 break;
1309                                         }
1310                                         pob = pob->parent;
1311                                 }
1312                         }
1313
1314                         if (zero_size) {
1315                                 unsigned int ml_count = 0;
1316                                 while (ml) {
1317                                         ml_count++;
1318                                         ml = ml->next;
1319                                 }
1320                                 process->totelem -= ml_count;
1321                         }
1322                         else {
1323                                 while (ml) {
1324                                         if (!(ml->flag & MB_HIDE)) {
1325                                                 int i;
1326                                                 float temp1[4][4], temp2[4][4], temp3[4][4];
1327                                                 float (*mat)[4] = NULL, (*imat)[4] = NULL;
1328                                                 float max_x, max_y, max_z, min_x, min_y, min_z;
1329                                                 float expx, expy, expz;
1330
1331                                                 max_x = max_y = max_z = -3.4e38;
1332                                                 min_x = min_y = min_z =  3.4e38;
1333
1334                                                 /* too big stiffness seems only ugly due to linear interpolation
1335                                                  * no need to have possibility for too big stiffness */
1336                                                 if (ml->s > 10.0f) ml->s = 10.0f;
1337
1338                                                 /* Rotation of MetaElem is stored in quat */
1339                                                 quat_to_mat4(temp3, ml->quat);
1340
1341                                                 /* Translation of MetaElem */
1342                                                 unit_m4(temp2);
1343                                                 temp2[3][0] = ml->x;
1344                                                 temp2[3][1] = ml->y;
1345                                                 temp2[3][2] = ml->z;
1346
1347                                                 mul_m4_m4m4(temp1, temp2, temp3);
1348
1349                                                 /* make a copy because of duplicates */
1350                                                 process->mainb[a] = new_pgn_element(process, sizeof(MetaElem));
1351                                                 *(process->mainb[a]) = *ml;
1352                                                 process->mainb[a]->bb = new_pgn_element(process, sizeof(BoundBox));
1353
1354                                                 mat = new_pgn_element(process, 4 * 4 * sizeof(float));
1355                                                 imat = new_pgn_element(process, 4 * 4 * sizeof(float));
1356
1357                                                 /* mat is the matrix to transform from mball into the basis-mball */
1358                                                 invert_m4_m4(obinv, obmat);
1359                                                 mul_m4_m4m4(temp2, obinv, bob->obmat);
1360                                                 /* MetaBall transformation */
1361                                                 mul_m4_m4m4(mat, temp2, temp1);
1362
1363                                                 invert_m4_m4(imat, mat);
1364
1365                                                 process->mainb[a]->rad2 = ml->rad * ml->rad;
1366
1367                                                 process->mainb[a]->mat = (float *) mat;
1368                                                 process->mainb[a]->imat = (float *) imat;
1369
1370                                                 if (!MB_TYPE_SIZE_SQUARED(ml->type)) {
1371                                                         expx = ml->expx;
1372                                                         expy = ml->expy;
1373                                                         expz = ml->expz;
1374                                                 }
1375                                                 else {
1376                                                         expx = ml->expx * ml->expx;
1377                                                         expy = ml->expy * ml->expy;
1378                                                         expz = ml->expz * ml->expz;
1379                                                 }
1380
1381                                                 /* untransformed Bounding Box of MetaElem */
1382                                                 /* TODO, its possible the elem type has been changed and the exp* values can use a fallback */
1383                                                 copy_v3_fl3(process->mainb[a]->bb->vec[0], -expx, -expy, -expz);  /* 0 */
1384                                                 copy_v3_fl3(process->mainb[a]->bb->vec[1], +expx, -expy, -expz);  /* 1 */
1385                                                 copy_v3_fl3(process->mainb[a]->bb->vec[2], +expx, +expy, -expz);  /* 2 */
1386                                                 copy_v3_fl3(process->mainb[a]->bb->vec[3], -expx, +expy, -expz);  /* 3 */
1387                                                 copy_v3_fl3(process->mainb[a]->bb->vec[4], -expx, -expy, +expz);  /* 4 */
1388                                                 copy_v3_fl3(process->mainb[a]->bb->vec[5], +expx, -expy, +expz);  /* 5 */
1389                                                 copy_v3_fl3(process->mainb[a]->bb->vec[6], +expx, +expy, +expz);  /* 6 */
1390                                                 copy_v3_fl3(process->mainb[a]->bb->vec[7], -expx, +expy, +expz);  /* 7 */
1391
1392                                                 /* transformation of Metalem bb */
1393                                                 for (i = 0; i < 8; i++)
1394                                                         mul_m4_v3((float (*)[4])mat, process->mainb[a]->bb->vec[i]);
1395
1396                                                 /* find max and min of transformed bb */
1397                                                 for (i = 0; i < 8; i++) {
1398                                                         /* find maximums */
1399                                                         if (process->mainb[a]->bb->vec[i][0] > max_x) max_x = process->mainb[a]->bb->vec[i][0];
1400                                                         if (process->mainb[a]->bb->vec[i][1] > max_y) max_y = process->mainb[a]->bb->vec[i][1];
1401                                                         if (process->mainb[a]->bb->vec[i][2] > max_z) max_z = process->mainb[a]->bb->vec[i][2];
1402                                                         /* find  minimums */
1403                                                         if (process->mainb[a]->bb->vec[i][0] < min_x) min_x = process->mainb[a]->bb->vec[i][0];
1404                                                         if (process->mainb[a]->bb->vec[i][1] < min_y) min_y = process->mainb[a]->bb->vec[i][1];
1405                                                         if (process->mainb[a]->bb->vec[i][2] < min_z) min_z = process->mainb[a]->bb->vec[i][2];
1406                                                 }
1407                                         
1408                                                 /* create "new" bb, only point 0 and 6, which are
1409                                                  * necessary for octal tree filling */
1410                                                 process->mainb[a]->bb->vec[0][0] = min_x - ml->rad;
1411                                                 process->mainb[a]->bb->vec[0][1] = min_y - ml->rad;
1412                                                 process->mainb[a]->bb->vec[0][2] = min_z - ml->rad;
1413
1414                                                 process->mainb[a]->bb->vec[6][0] = max_x + ml->rad;
1415                                                 process->mainb[a]->bb->vec[6][1] = max_y + ml->rad;
1416                                                 process->mainb[a]->bb->vec[6][2] = max_z + ml->rad;
1417
1418                                                 a++;
1419                                         }
1420                                         ml = ml->next;
1421                                 }
1422                         }
1423                 }
1424         }
1425
1426         
1427         /* totsize (= 'manhattan' radius) */
1428         totsize = 0.0;
1429         for (a = 0; a < process->totelem; a++) {
1430                 
1431                 vec[0] = process->mainb[a]->x + process->mainb[a]->rad + process->mainb[a]->expx;
1432                 vec[1] = process->mainb[a]->y + process->mainb[a]->rad + process->mainb[a]->expy;
1433                 vec[2] = process->mainb[a]->z + process->mainb[a]->rad + process->mainb[a]->expz;
1434
1435                 calc_mballco(process->mainb[a], vec);
1436         
1437                 size = fabsf(vec[0]);
1438                 if (size > totsize) totsize = size;
1439                 size = fabsf(vec[1]);
1440                 if (size > totsize) totsize = size;
1441                 size = fabsf(vec[2]);
1442                 if (size > totsize) totsize = size;
1443
1444                 vec[0] = process->mainb[a]->x - process->mainb[a]->rad;
1445                 vec[1] = process->mainb[a]->y - process->mainb[a]->rad;
1446                 vec[2] = process->mainb[a]->z - process->mainb[a]->rad;
1447                                 
1448                 calc_mballco(process->mainb[a], vec);
1449         
1450                 size = fabsf(vec[0]);
1451                 if (size > totsize) totsize = size;
1452                 size = fabsf(vec[1]);
1453                 if (size > totsize) totsize = size;
1454                 size = fabsf(vec[2]);
1455                 if (size > totsize) totsize = size;
1456         }
1457
1458         for (a = 0; a < process->totelem; a++) {
1459                 process->thresh += densfunc(process->mainb[a], 2.0f * totsize, 2.0f * totsize, 2.0f * totsize);
1460         }
1461
1462         return totsize;
1463 }
1464
1465 /* if MetaElem lies in node, then node includes MetaElem pointer (ml_p)
1466  * pointing at MetaElem (ml)
1467  */
1468 static void fill_metaball_octal_node(octal_node *node, MetaElem *ml, short i)
1469 {
1470         ml_pointer *ml_p;
1471
1472         ml_p = MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
1473         ml_p->ml = ml;
1474         BLI_addtail(&(node->nodes[i]->elems), ml_p);
1475         node->count++;
1476         
1477         if ((ml->flag & MB_NEGATIVE) == 0) {
1478                 node->nodes[i]->pos++;
1479         }
1480         else {
1481                 node->nodes[i]->neg++;
1482         }
1483 }
1484
1485 /* Node is subdivided as is illustrated on the following figure:
1486  * 
1487  *      +------+------+
1488  *     /      /      /|
1489  *    +------+------+ |
1490  *   /      /      /| +
1491  *  +------+------+ |/|
1492  *  |      |      | + |
1493  *  |      |      |/| +
1494  *  +------+------+ |/
1495  *  |      |      | +
1496  *  |      |      |/
1497  *  +------+------+
1498  *  
1499  */
1500 static void subdivide_metaball_octal_node(octal_node *node, float size_x, float size_y, float size_z, short depth)
1501 {
1502         MetaElem *ml;
1503         ml_pointer *ml_p;
1504         float x, y, z;
1505         int a, i;
1506
1507         /* create new nodes */
1508         for (a = 0; a < 8; a++) {
1509                 node->nodes[a] = MEM_mallocN(sizeof(octal_node), "octal_node");
1510                 for (i = 0; i < 8; i++)
1511                         node->nodes[a]->nodes[i] = NULL;
1512                 node->nodes[a]->parent = node;
1513                 BLI_listbase_clear(&node->nodes[a]->elems);
1514                 node->nodes[a]->count = 0;
1515                 node->nodes[a]->neg = 0;
1516                 node->nodes[a]->pos = 0;
1517         }
1518
1519         size_x /= 2;
1520         size_y /= 2;
1521         size_z /= 2;
1522         
1523         /* center of node */
1524         node->x = x = node->x_min + size_x;
1525         node->y = y = node->y_min + size_y;
1526         node->z = z = node->z_min + size_z;
1527
1528         /* setting up of border points of new nodes */
1529         node->nodes[0]->x_min = node->x_min;
1530         node->nodes[0]->y_min = node->y_min;
1531         node->nodes[0]->z_min = node->z_min;
1532         node->nodes[0]->x = node->nodes[0]->x_min + size_x / 2;
1533         node->nodes[0]->y = node->nodes[0]->y_min + size_y / 2;
1534         node->nodes[0]->z = node->nodes[0]->z_min + size_z / 2;
1535         
1536         node->nodes[1]->x_min = x;
1537         node->nodes[1]->y_min = node->y_min;
1538         node->nodes[1]->z_min = node->z_min;
1539         node->nodes[1]->x = node->nodes[1]->x_min + size_x / 2;
1540         node->nodes[1]->y = node->nodes[1]->y_min + size_y / 2;
1541         node->nodes[1]->z = node->nodes[1]->z_min + size_z / 2;
1542
1543         node->nodes[2]->x_min = x;
1544         node->nodes[2]->y_min = y;
1545         node->nodes[2]->z_min = node->z_min;
1546         node->nodes[2]->x = node->nodes[2]->x_min + size_x / 2;
1547         node->nodes[2]->y = node->nodes[2]->y_min + size_y / 2;
1548         node->nodes[2]->z = node->nodes[2]->z_min + size_z / 2;
1549
1550         node->nodes[3]->x_min = node->x_min;
1551         node->nodes[3]->y_min = y;
1552         node->nodes[3]->z_min = node->z_min;
1553         node->nodes[3]->x = node->nodes[3]->x_min + size_x / 2;
1554         node->nodes[3]->y = node->nodes[3]->y_min + size_y / 2;
1555         node->nodes[3]->z = node->nodes[3]->z_min + size_z / 2;
1556
1557         node->nodes[4]->x_min = node->x_min;
1558         node->nodes[4]->y_min = node->y_min;
1559         node->nodes[4]->z_min = z;
1560         node->nodes[4]->x = node->nodes[4]->x_min + size_x / 2;
1561         node->nodes[4]->y = node->nodes[4]->y_min + size_y / 2;
1562         node->nodes[4]->z = node->nodes[4]->z_min + size_z / 2;
1563         
1564         node->nodes[5]->x_min = x;
1565         node->nodes[5]->y_min = node->y_min;
1566         node->nodes[5]->z_min = z;
1567         node->nodes[5]->x = node->nodes[5]->x_min + size_x / 2;
1568         node->nodes[5]->y = node->nodes[5]->y_min + size_y / 2;
1569         node->nodes[5]->z = node->nodes[5]->z_min + size_z / 2;
1570
1571         node->nodes[6]->x_min = x;
1572         node->nodes[6]->y_min = y;
1573         node->nodes[6]->z_min = z;
1574         node->nodes[6]->x = node->nodes[6]->x_min + size_x / 2;
1575         node->nodes[6]->y = node->nodes[6]->y_min + size_y / 2;
1576         node->nodes[6]->z = node->nodes[6]->z_min + size_z / 2;
1577
1578         node->nodes[7]->x_min = node->x_min;
1579         node->nodes[7]->y_min = y;
1580         node->nodes[7]->z_min = z;
1581         node->nodes[7]->x = node->nodes[7]->x_min + size_x / 2;
1582         node->nodes[7]->y = node->nodes[7]->y_min + size_y / 2;
1583         node->nodes[7]->z = node->nodes[7]->z_min + size_z / 2;
1584
1585         ml_p = node->elems.first;
1586         
1587         /* setting up references of MetaElems for new nodes */
1588         while (ml_p) {
1589                 ml = ml_p->ml;
1590                 if (ml->bb->vec[0][2] < z) {
1591                         if (ml->bb->vec[0][1] < y) {
1592                                 /* vec[0][0] lies in first octant */
1593                                 if (ml->bb->vec[0][0] < x) {
1594                                         /* ml belongs to the (0)1st node */
1595                                         fill_metaball_octal_node(node, ml, 0);
1596
1597                                         /* ml belongs to the (3)4th node */
1598                                         if (ml->bb->vec[6][1] >= y) {
1599                                                 fill_metaball_octal_node(node, ml, 3);
1600
1601                                                 /* ml belongs to the (7)8th node */
1602                                                 if (ml->bb->vec[6][2] >= z) {
1603                                                         fill_metaball_octal_node(node, ml, 7);
1604                                                 }
1605                                         }
1606         
1607                                         /* ml belongs to the (1)2nd node */
1608                                         if (ml->bb->vec[6][0] >= x) {
1609                                                 fill_metaball_octal_node(node, ml, 1);
1610
1611                                                 /* ml belongs to the (5)6th node */
1612                                                 if (ml->bb->vec[6][2] >= z) {
1613                                                         fill_metaball_octal_node(node, ml, 5);
1614                                                 }
1615                                         }
1616
1617                                         /* ml belongs to the (2)3th node */
1618                                         if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
1619                                                 fill_metaball_octal_node(node, ml, 2);
1620                                                 
1621                                                 /* ml belong to the (6)7th node */
1622                                                 if (ml->bb->vec[6][2] >= z) {
1623                                                         fill_metaball_octal_node(node, ml, 6);
1624                                                 }
1625                                                 
1626                                         }
1627                         
1628                                         /* ml belongs to the (4)5th node too */
1629                                         if (ml->bb->vec[6][2] >= z) {
1630                                                 fill_metaball_octal_node(node, ml, 4);
1631                                         }
1632
1633                                         
1634                                         
1635                                 }
1636                                 /* vec[0][0] is in the (1)second octant */
1637                                 else {
1638                                         /* ml belong to the (1)2nd node */
1639                                         fill_metaball_octal_node(node, ml, 1);
1640
1641                                         /* ml belongs to the (2)3th node */
1642                                         if (ml->bb->vec[6][1] >= y) {
1643                                                 fill_metaball_octal_node(node, ml, 2);
1644
1645                                                 /* ml belongs to the (6)7th node */
1646                                                 if (ml->bb->vec[6][2] >= z) {
1647                                                         fill_metaball_octal_node(node, ml, 6);
1648                                                 }
1649                                                 
1650                                         }
1651                                         
1652                                         /* ml belongs to the (5)6th node */
1653                                         if (ml->bb->vec[6][2] >= z) {
1654                                                 fill_metaball_octal_node(node, ml, 5);
1655                                         }
1656                                 }
1657                         }
1658                         else {
1659                                 /* vec[0][0] is in the (3)4th octant */
1660                                 if (ml->bb->vec[0][0] < x) {
1661                                         /* ml belongs to the (3)4nd node */
1662                                         fill_metaball_octal_node(node, ml, 3);
1663                                         
1664                                         /* ml belongs to the (7)8th node */
1665                                         if (ml->bb->vec[6][2] >= z) {
1666                                                 fill_metaball_octal_node(node, ml, 7);
1667                                         }
1668                                 
1669
1670                                         /* ml belongs to the (2)3th node */
1671                                         if (ml->bb->vec[6][0] >= x) {
1672                                                 fill_metaball_octal_node(node, ml, 2);
1673                                         
1674                                                 /* ml belongs to the (6)7th node */
1675                                                 if (ml->bb->vec[6][2] >= z) {
1676                                                         fill_metaball_octal_node(node, ml, 6);
1677                                                 }
1678                                         }
1679                                 }
1680
1681                         }
1682
1683                         /* vec[0][0] is in the (2)3th octant */
1684                         if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
1685                                 /* ml belongs to the (2)3th node */
1686                                 fill_metaball_octal_node(node, ml, 2);
1687                                 
1688                                 /* ml belongs to the (6)7th node */
1689                                 if (ml->bb->vec[6][2] >= z) {
1690                                         fill_metaball_octal_node(node, ml, 6);
1691                                 }
1692                         }
1693                 }
1694                 else {
1695                         if (ml->bb->vec[0][1] < y) {
1696                                 /* vec[0][0] lies in (4)5th octant */
1697                                 if (ml->bb->vec[0][0] < x) {
1698                                         /* ml belongs to the (4)5th node */
1699                                         fill_metaball_octal_node(node, ml, 4);
1700
1701                                         if (ml->bb->vec[6][0] >= x) {
1702                                                 fill_metaball_octal_node(node, ml, 5);
1703                                         }
1704
1705                                         if (ml->bb->vec[6][1] >= y) {
1706                                                 fill_metaball_octal_node(node, ml, 7);
1707                                         }
1708                                         
1709                                         if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
1710                                                 fill_metaball_octal_node(node, ml, 6);
1711                                         }
1712                                 }
1713                                 /* vec[0][0] lies in (5)6th octant */
1714                                 else {
1715                                         fill_metaball_octal_node(node, ml, 5);
1716
1717                                         if (ml->bb->vec[6][1] >= y) {
1718                                                 fill_metaball_octal_node(node, ml, 6);
1719                                         }
1720                                 }
1721                         }
1722                         else {
1723                                 /* vec[0][0] lies in (7)8th octant */
1724                                 if (ml->bb->vec[0][0] < x) {
1725                                         fill_metaball_octal_node(node, ml, 7);
1726
1727                                         if (ml->bb->vec[6][0] >= x) {
1728                                                 fill_metaball_octal_node(node, ml, 6);
1729                                         }
1730                                 }
1731
1732                         }
1733                         
1734                         /* vec[0][0] lies in (6)7th octant */
1735                         if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
1736                                 fill_metaball_octal_node(node, ml, 6);
1737                         }
1738                 }
1739                 ml_p = ml_p->next;
1740         }
1741
1742         /* free references of MetaElems for curent node (it is not needed anymore) */
1743         BLI_freelistN(&node->elems);
1744
1745         depth--;
1746         
1747         if (depth > 0) {
1748                 for (a = 0; a < 8; a++) {
1749                         if (node->nodes[a]->count > 0) /* if node is not empty, then it is subdivided */
1750                                 subdivide_metaball_octal_node(node->nodes[a], size_x, size_y, size_z, depth);
1751                 }
1752         }
1753 }
1754
1755 /* free all octal nodes recursively */
1756 static void free_metaball_octal_node(octal_node *node)
1757 {
1758         int a;
1759         for (a = 0; a < 8; a++) {
1760                 if (node->nodes[a] != NULL) free_metaball_octal_node(node->nodes[a]);
1761         }
1762         BLI_freelistN(&node->elems);
1763         MEM_freeN(node);
1764 }
1765
1766 /* If scene include more than one MetaElem, then octree is used */
1767 static void init_metaball_octal_tree(PROCESS *process, int depth)
1768 {
1769         struct octal_node *node;
1770         ml_pointer *ml_p;
1771         float size[3];
1772         int a;
1773         
1774         process->metaball_tree = MEM_mallocN(sizeof(octal_tree), "metaball_octal_tree");
1775         process->metaball_tree->first = node = MEM_mallocN(sizeof(octal_node), "metaball_octal_node");
1776         /* maximal depth of octree */
1777         process->metaball_tree->depth = depth;
1778
1779         process->metaball_tree->neg = node->neg = 0;
1780         process->metaball_tree->pos = node->pos = 0;
1781         
1782         BLI_listbase_clear(&node->elems);
1783         node->count = 0;
1784
1785         for (a = 0; a < 8; a++)
1786                 node->nodes[a] = NULL;
1787
1788         node->x_min = node->y_min = node->z_min = FLT_MAX;
1789         node->x_max = node->y_max = node->z_max = -FLT_MAX;
1790
1791         /* size of octal tree scene */
1792         for (a = 0; a < process->totelem; a++) {
1793                 if (process->mainb[a]->bb->vec[0][0] < node->x_min) node->x_min = process->mainb[a]->bb->vec[0][0];
1794                 if (process->mainb[a]->bb->vec[0][1] < node->y_min) node->y_min = process->mainb[a]->bb->vec[0][1];
1795                 if (process->mainb[a]->bb->vec[0][2] < node->z_min) node->z_min = process->mainb[a]->bb->vec[0][2];
1796
1797                 if (process->mainb[a]->bb->vec[6][0] > node->x_max) node->x_max = process->mainb[a]->bb->vec[6][0];
1798                 if (process->mainb[a]->bb->vec[6][1] > node->y_max) node->y_max = process->mainb[a]->bb->vec[6][1];
1799                 if (process->mainb[a]->bb->vec[6][2] > node->z_max) node->z_max = process->mainb[a]->bb->vec[6][2];
1800
1801                 ml_p = MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
1802                 ml_p->ml = process->mainb[a];
1803                 BLI_addtail(&node->elems, ml_p);
1804
1805                 if ((process->mainb[a]->flag & MB_NEGATIVE) == 0) {
1806                         /* number of positive MetaElem in scene */
1807                         process->metaball_tree->pos++;
1808                 }
1809                 else {
1810                         /* number of negative MetaElem in scene */
1811                         process->metaball_tree->neg++;
1812                 }
1813         }
1814
1815         /* size of first node */
1816         size[0] = node->x_max - node->x_min;
1817         size[1] = node->y_max - node->y_min;
1818         size[2] = node->z_max - node->z_min;
1819
1820         /* first node is subdivided recursively */
1821         subdivide_metaball_octal_node(node, size[0], size[1], size[2], process->metaball_tree->depth);
1822 }
1823
1824 static void mball_count(EvaluationContext *eval_ctx, PROCESS *process, Scene *scene, Object *basis)
1825 {
1826         Scene *sce_iter = scene;
1827         Base *base;
1828         Object *ob, *bob = basis;
1829         MetaElem *ml = NULL;
1830         int basisnr, obnr;
1831         char basisname[MAX_ID_NAME], obname[MAX_ID_NAME];
1832         SceneBaseIter iter;
1833
1834         BLI_split_name_num(basisname, &basisnr, basis->id.name + 2, '.');
1835         process->totelem = 0;
1836
1837         BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 0, NULL, NULL);
1838         while (BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 1, &base, &ob)) {
1839                 if (ob->type == OB_MBALL) {
1840                         if (ob == bob) {
1841                                 MetaBall *mb = ob->data;
1842
1843                                 /* if bob object is in edit mode, then dynamic list of all MetaElems
1844                                  * is stored in editelems */
1845                                 if (mb->editelems) ml = mb->editelems->first;
1846                                 /* if bob object is in object mode */
1847                                 else ml = mb->elems.first;
1848                         }
1849                         else {
1850                                 BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
1851
1852                                 /* object ob has to be in same "group" ... it means, that it has to have
1853                                  * same base of its name */
1854                                 if (STREQ(obname, basisname)) {
1855                                         MetaBall *mb = ob->data;
1856
1857                                         /* if object is in edit mode, then dynamic list of all MetaElems
1858                                          * is stored in editelems */
1859                                         if (mb->editelems) ml = mb->editelems->first;
1860                                         /* if bob object is in object mode */
1861                                         else ml = mb->elems.first;
1862                                 }
1863                         }
1864
1865                         for ( ; ml; ml = ml->next) {
1866                                 if (!(ml->flag & MB_HIDE)) {
1867                                         process->totelem++;
1868                                 }
1869                         }
1870                 }
1871         }
1872 }
1873
1874 void BKE_mball_polygonize(EvaluationContext *eval_ctx, Scene *scene, Object *ob, ListBase *dispbase)
1875 {
1876         MetaBall *mb;
1877         DispList *dl;
1878         int a, nr_cubes;
1879         float *co, *no, totsize, width;
1880         PROCESS process = {0};
1881
1882         mb = ob->data;
1883
1884         mball_count(eval_ctx, &process, scene, ob);
1885
1886         if (process.totelem == 0) return;
1887         if ((eval_ctx->mode != DAG_EVAL_RENDER) && (mb->flag == MB_UPDATE_NEVER)) return;
1888         if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_FAST) return;
1889
1890         process.thresh = mb->thresh;
1891
1892         /* total number of MetaElems (totelem) is precomputed in find_basis_mball() function */
1893         process.mainb = MEM_mallocN(sizeof(void *) * process.totelem, "mainb");
1894         
1895         /* initialize all mainb (MetaElems) */
1896         totsize = init_meta(eval_ctx, &process, scene, ob);
1897
1898         /* if scene includes more than one MetaElem, then octal tree optimization is used */
1899         if ((process.totelem >    1) && (process.totelem <=   64)) init_metaball_octal_tree(&process, 1);
1900         if ((process.totelem >   64) && (process.totelem <=  128)) init_metaball_octal_tree(&process, 2);
1901         if ((process.totelem >  128) && (process.totelem <=  512)) init_metaball_octal_tree(&process, 3);
1902         if ((process.totelem >  512) && (process.totelem <= 1024)) init_metaball_octal_tree(&process, 4);
1903         if (process.totelem  > 1024)                               init_metaball_octal_tree(&process, 5);
1904
1905         /* don't polygonize metaballs with too high resolution (base mball to small)
1906          * note: Eps was 0.0001f but this was giving problems for blood animation for durian, using 0.00001f */
1907         if (process.metaball_tree) {
1908                 if (ob->size[0] <= 0.00001f * (process.metaball_tree->first->x_max - process.metaball_tree->first->x_min) ||
1909                     ob->size[1] <= 0.00001f * (process.metaball_tree->first->y_max - process.metaball_tree->first->y_min) ||
1910                     ob->size[2] <= 0.00001f * (process.metaball_tree->first->z_max - process.metaball_tree->first->z_min))
1911                 {
1912                         new_pgn_element(&process, -1); /* free values created by init_meta */
1913
1914                         MEM_freeN(process.mainb);
1915
1916                         /* free tree */
1917                         free_metaball_octal_node(process.metaball_tree->first);
1918                         MEM_freeN(process.metaball_tree);
1919
1920                         return;
1921                 }
1922         }
1923
1924         /* width is size per polygonize cube */
1925         if (eval_ctx->mode == DAG_EVAL_RENDER) {
1926                 width = mb->rendersize;
1927         }
1928         else {
1929                 width = mb->wiresize;
1930                 if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_HALFRES) {
1931                         width *= 2;
1932                 }
1933         }
1934         /* nr_cubes is just for safety, minimum is totsize */
1935         nr_cubes = (int)(0.5f + totsize / width);
1936
1937         /* init process */
1938         process.function = metaball;
1939         process.size = width;
1940         process.bounds = nr_cubes;
1941         process.cubes = NULL;
1942         process.delta = width / (float)(RES * RES);
1943
1944         polygonize(&process, mb);
1945         
1946         MEM_freeN(process.mainb);
1947
1948         /* free octal tree */
1949         if (process.totelem > 1) {
1950                 free_metaball_octal_node(process.metaball_tree->first);
1951                 MEM_freeN(process.metaball_tree);
1952                 process.metaball_tree = NULL;
1953         }
1954
1955         if (process.curindex) {
1956                 VERTEX *ptr = process.vertices.ptr;
1957
1958                 dl = MEM_callocN(sizeof(DispList), "mbaldisp");
1959                 BLI_addtail(dispbase, dl);
1960                 dl->type = DL_INDEX4;
1961                 dl->nr = process.vertices.count;
1962                 dl->parts = process.curindex;
1963
1964                 dl->index = process.indices;
1965                 process.indices = NULL;
1966
1967                 a = process.vertices.count;
1968                 dl->verts = co = MEM_mallocN(sizeof(float[3]) * a, "mballverts");
1969                 dl->nors = no = MEM_mallocN(sizeof(float[3]) * a, "mballnors");
1970
1971                 for (a = 0; a < process.vertices.count; ptr++, a++, no += 3, co += 3) {
1972                         copy_v3_v3(co, ptr->co);
1973                         copy_v3_v3(no, ptr->no);
1974                 }
1975         }
1976
1977         freepolygonize(&process);
1978 }