BKE_mball: split tessellation into its own file
authorCampbell Barton <ideasman42@gmail.com>
Thu, 26 Feb 2015 03:39:57 +0000 (14:39 +1100)
committerCampbell Barton <ideasman42@gmail.com>
Wed, 25 Mar 2015 11:40:09 +0000 (22:40 +1100)
this has a lot of its own local structs, functions,
better to keep isolated from general metaball selection/library logic.

source/blender/blenkernel/BKE_mball_tessellate.h [new file with mode: 0644]
source/blender/blenkernel/CMakeLists.txt
source/blender/blenkernel/intern/mball.c
source/blender/blenkernel/intern/mball_tessellate.c [new file with mode: 0644]

diff --git a/source/blender/blenkernel/BKE_mball_tessellate.h b/source/blender/blenkernel/BKE_mball_tessellate.h
new file mode 100644 (file)
index 0000000..e69de29
index a826fe26e68c7612d2fc674cfc15a4a1589f518f..f19c1066766ebd23b56db5f28d0be176b780d70a 100644 (file)
@@ -119,6 +119,7 @@ set(SRC
        intern/mask_rasterize.c
        intern/material.c
        intern/mball.c
+       intern/mball_tessellate.c
        intern/mesh.c
        intern/mesh_evaluate.c
        intern/mesh_mapping.c
@@ -233,6 +234,7 @@ set(SRC
        BKE_mask.h
        BKE_material.h
        BKE_mball.h
+       BKE_mball_tessellate.h
        BKE_mesh.h
        BKE_mesh_mapping.h
        BKE_mesh_remap.h
index ce20636233a6bbfe3f63cd986795a0b3ab8eadbc..efc5b23a9d02e47d0594c3cb30b99af026d7dbe6 100644 (file)
 #include "BKE_object.h"
 #include "BKE_material.h"
 
-/* Data types */
-
-typedef struct vertex {         /* surface vertex */
-       float co[3];  /* position and surface normal */
-       float no[3];
-} VERTEX;
-
-typedef struct vertices {       /* list of vertices in polygonization */
-       int count, max;             /* # vertices, max # allowed */
-       VERTEX *ptr;                /* dynamically allocated */
-} VERTICES;
-
-typedef struct corner {         /* corner of a cube */
-       int i, j, k;                /* (i, j, k) is index within lattice */
-       float co[3], value;       /* location and function value */
-       struct corner *next;
-} CORNER;
-
-typedef struct cube {           /* partitioning cell (cube) */
-       int i, j, k;                /* lattice location of cube */
-       CORNER *corners[8];         /* eight corners */
-} CUBE;
-
-typedef struct cubes {          /* linked list of cubes acting as stack */
-       CUBE cube;                  /* a single cube */
-       struct cubes *next;         /* remaining elements */
-} CUBES;
-
-typedef struct centerlist {     /* list of cube locations */
-       int i, j, k;                /* cube location */
-       struct centerlist *next;    /* remaining elements */
-} CENTERLIST;
-
-typedef struct edgelist {       /* list of edges */
-       int i1, j1, k1, i2, j2, k2; /* edge corner ids */
-       int vid;                    /* vertex id */
-       struct edgelist *next;      /* remaining elements */
-} EDGELIST;
-
-typedef struct intlist {        /* list of integers */
-       int i;                      /* an integer */
-       struct intlist *next;       /* remaining elements */
-} INTLIST;
-
-typedef struct intlists {       /* list of list of integers */
-       INTLIST *list;              /* a list of integers */
-       struct intlists *next;      /* remaining elements */
-} INTLISTS;
-
-/* dividing scene using octal tree makes polygonisation faster */
-typedef struct ml_pointer {
-       struct ml_pointer *next, *prev;
-       struct MetaElem *ml;
-} ml_pointer;
-
-typedef struct octal_node {
-       struct octal_node *nodes[8];/* children of current node */
-       struct octal_node *parent;  /* parent of current node */
-       struct ListBase elems;      /* ListBase of MetaElem pointers (ml_pointer) */
-       float x_min, y_min, z_min;  /* 1st border point */
-       float x_max, y_max, z_max;  /* 7th border point */
-       float x, y, z;              /* center of node */
-       int pos, neg;               /* number of positive and negative MetaElements in the node */
-       int count;                  /* number of MetaElems, which belongs to the node */
-} octal_node;
-
-typedef struct octal_tree {
-       struct octal_node *first;   /* first node */
-       int pos, neg;               /* number of positive and negative MetaElements in the scene */
-       short depth;                /* number of scene subdivision */
-} octal_tree;
-
-struct pgn_elements {
-       struct pgn_elements *next, *prev;
-       char *data;
-};
-
-typedef struct process {        /* parameters, function, storage */
-       /* ** old G_mb contents ** */
-       float thresh;
-       int totelem;
-       MetaElem **mainb;
-       octal_tree *metaball_tree;
-
-       /* ** old process contents ** */
-
-       /* what happens here? floats, I think. */
-       /*  float (*function)(void);     */     /* implicit surface function */
-       float (*function)(struct process *, float, float, float);
-       float size, delta;          /* cube size, normal delta */
-       int bounds;                 /* cube range within lattice */
-       CUBES *cubes;               /* active cubes */
-       VERTICES vertices;          /* surface vertices */
-       CENTERLIST **centers;       /* cube center hash table */
-       CORNER **corners;           /* corner value hash table */
-       EDGELIST **edges;           /* edge and vertex id hash table */
-
-       /* Runtime things */
-       int *indices;
-       int totindex, curindex;
-
-       int pgn_offset;
-       struct pgn_elements *pgn_current;
-       ListBase pgn_list;
-} PROCESS;
-
-/* Forward declarations */
-static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2, MetaBall *mb);
-static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const int j, const int k);
-static CORNER *setcorner(PROCESS *process, int i, int j, int k);
-static void converge(PROCESS *process, const float p1[3], const float p2[3], float v1, float v2,
-                     float p[3], MetaBall *mb, int f);
-
 /* Functions */
 
 void BKE_mball_unlink(MetaBall *mb)
@@ -558,1817 +445,6 @@ Object *BKE_mball_basis_find(Scene *scene, Object *basis)
        return basis;
 }
 
-
-/* ******************** ARITH ************************* */
-
-/* BASED AT CODE (but mostly rewritten) :
- * C code from the article
- * "An Implicit Surface Polygonizer"
- * by Jules Bloomenthal, jbloom@beauty.gmu.edu
- * in "Graphics Gems IV", Academic Press, 1994
- *
- * Authored by Jules Bloomenthal, Xerox PARC.
- * Copyright (c) Xerox Corporation, 1991.  All rights reserved.
- * Permission is granted to reproduce, use and distribute this code for
- * any and all purposes, provided that this notice appears in all copies. */
-
-#define RES 12 /* # converge iterations    */
-
-#define L   0  /* left direction:      -x, -i */
-#define R   1  /* right direction:     +x, +i */
-#define B   2  /* bottom direction: -y, -j */
-#define T   3  /* top direction:       +y, +j */
-#define N   4  /* near direction:      -z, -k */
-#define F   5  /* far direction:       +z, +k */
-#define LBN 0  /* left bottom near corner  */
-#define LBF 1  /* left bottom far corner   */
-#define LTN 2  /* left top near corner     */
-#define LTF 3  /* left top far corner      */
-#define RBN 4  /* right bottom near corner */
-#define RBF 5  /* right bottom far corner  */
-#define RTN 6  /* right top near corner    */
-#define RTF 7  /* right top far corner     */
-
-/* the LBN corner of cube (i, j, k), corresponds with location
- * (i-0.5)*size, (j-0.5)*size, (k-0.5)*size) */
-
-#define HASHBIT     (5)
-#define HASHSIZE    (size_t)(1 << (3 * HASHBIT))   /*! < hash table size (32768) */
-
-#define HASH(i, j, k) ((((( (i) & 31) << 5) | ( (j) & 31)) << 5) | ( (k) & 31) )
-
-#define MB_BIT(i, bit) (((i) >> (bit)) & 1)
-#define FLIP(i, bit) ((i) ^ 1 << (bit)) /* flip the given bit of i */
-
-
-/* **************** POLYGONIZATION ************************ */
-
-static void calc_mballco(MetaElem *ml, float vec[3])
-{
-       if (ml->mat) {
-               mul_m4_v3((float (*)[4])ml->mat, vec);
-       }
-}
-
-static float densfunc(MetaElem *ball, float x, float y, float z)
-{
-       float dist2;
-       float dvec[3] = {x, y, z};
-
-       mul_m4_v3((float (*)[4])ball->imat, dvec);
-
-       switch (ball->type) {
-               case MB_BALL:
-                       /* do nothing */
-                       break;
-               case MB_TUBE:
-                       if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
-                       else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
-                       else                            dvec[0] = 0.0;
-                       break;
-               case MB_PLANE:
-                       if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
-                       else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
-                       else                            dvec[0] = 0.0;
-                       if      (dvec[1] >  ball->expy) dvec[1] -= ball->expy;
-                       else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
-                       else                            dvec[1] = 0.0;
-                       break;
-               case MB_ELIPSOID:
-                       dvec[0] /= ball->expx;
-                       dvec[1] /= ball->expy;
-                       dvec[2] /= ball->expz;
-                       break;
-               case MB_CUBE:
-                       if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
-                       else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
-                       else                            dvec[0] = 0.0;
-
-                       if      (dvec[1] >  ball->expy) dvec[1] -= ball->expy;
-                       else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
-                       else                            dvec[1] = 0.0;
-
-                       if      (dvec[2] >  ball->expz) dvec[2] -= ball->expz;
-                       else if (dvec[2] < -ball->expz) dvec[2] += ball->expz;
-                       else                            dvec[2] = 0.0;
-                       break;
-
-               /* *** deprecated, could be removed?, do-versioned at least *** */
-               case MB_TUBEX:
-                       if      (dvec[0] >  ball->len) dvec[0] -= ball->len;
-                       else if (dvec[0] < -ball->len) dvec[0] += ball->len;
-                       else                           dvec[0] = 0.0;
-                       break;
-               case MB_TUBEY:
-                       if      (dvec[1] >  ball->len) dvec[1] -= ball->len;
-                       else if (dvec[1] < -ball->len) dvec[1] += ball->len;
-                       else                           dvec[1] = 0.0;
-                       break;
-               case MB_TUBEZ:
-                       if      (dvec[2] >  ball->len) dvec[2] -= ball->len;
-                       else if (dvec[2] < -ball->len) dvec[2] += ball->len;
-                       else                           dvec[2] = 0.0;
-                       break;
-               /* *** end deprecated *** */
-       }
-
-       dist2 = 1.0f - (len_squared_v3(dvec) / ball->rad2);
-
-       if ((ball->flag & MB_NEGATIVE) == 0) {
-               return (dist2 < 0.0f) ? -0.5f : (ball->s * dist2 * dist2 * dist2) - 0.5f;
-       }
-       else {
-               return (dist2 < 0.0f) ? 0.5f : 0.5f - (ball->s * dist2 * dist2 * dist2);
-       }
-}
-
-static octal_node *find_metaball_octal_node(octal_node *node, float x, float y, float z, short depth)
-{
-       if (!depth) return node;
-       
-       if (z < node->z) {
-               if (y < node->y) {
-                       if (x < node->x) {
-                               if (node->nodes[0])
-                                       return find_metaball_octal_node(node->nodes[0], x, y, z, depth--);
-                               else
-                                       return node;
-                       }
-                       else {
-                               if (node->nodes[1])
-                                       return find_metaball_octal_node(node->nodes[1], x, y, z, depth--);
-                               else
-                                       return node;
-                       }
-               }
-               else {
-                       if (x < node->x) {
-                               if (node->nodes[3])
-                                       return find_metaball_octal_node(node->nodes[3], x, y, z, depth--);
-                               else
-                                       return node;
-                       }
-                       else {
-                               if (node->nodes[2])
-                                       return find_metaball_octal_node(node->nodes[2], x, y, z, depth--);
-                               else
-                                       return node;
-                       }
-               }
-       }
-       else {
-               if (y < node->y) {
-                       if (x < node->x) {
-                               if (node->nodes[4])
-                                       return find_metaball_octal_node(node->nodes[4], x, y, z, depth--);
-                               else
-                                       return node;
-                       }
-                       else {
-                               if (node->nodes[5])
-                                       return find_metaball_octal_node(node->nodes[5], x, y, z, depth--);
-                               else
-                                       return node;
-                       }
-               }
-               else {
-                       if (x < node->x) {
-                               if (node->nodes[7])
-                                       return find_metaball_octal_node(node->nodes[7], x, y, z, depth--);
-                               else
-                                       return node;
-                       }
-                       else {
-                               if (node->nodes[6])
-                                       return find_metaball_octal_node(node->nodes[6], x, y, z, depth--);
-                               else
-                                       return node;
-                       }
-               }
-       }
-       
-       /* all cases accounted for */
-       BLI_assert(0);
-}
-
-static float metaball(PROCESS *process, float x, float y, float z)
-/*  float x, y, z; */
-{
-       octal_tree *metaball_tree = process->metaball_tree;
-       struct octal_node *node;
-       struct ml_pointer *ml_p;
-       float dens = 0;
-       int a;
-       
-       if (process->totelem > 1) {
-               node = find_metaball_octal_node(metaball_tree->first, x, y, z, metaball_tree->depth);
-               if (node) {
-                       for (ml_p = node->elems.first; ml_p; ml_p = ml_p->next) {
-                               dens += densfunc(ml_p->ml, x, y, z);
-                       }
-
-                       dens += -0.5f * (metaball_tree->pos - node->pos);
-                       dens +=  0.5f * (metaball_tree->neg - node->neg);
-               }
-               else {
-                       for (a = 0; a < process->totelem; a++) {
-                               dens += densfunc(process->mainb[a], x, y, z);
-                       }
-               }
-       }
-       else {
-               dens += densfunc(process->mainb[0], x, y, z);
-       }
-
-       return process->thresh - dens;
-}
-
-/* ******************************************** */
-
-static void accum_mballfaces(PROCESS *process, int i1, int i2, int i3, int i4)
-{
-       int *newi, *cur;
-       /* static int i = 0; I would like to delete altogether, but I don't dare to, yet */
-
-       if (process->totindex == process->curindex) {
-               process->totindex += 256;
-               newi = MEM_mallocN(4 * sizeof(int) * process->totindex, "vertindex");
-               
-               if (process->indices) {
-                       memcpy(newi, process->indices, 4 * sizeof(int) * (process->totindex - 256));
-                       MEM_freeN(process->indices);
-               }
-               process->indices = newi;
-       }
-       
-       cur = process->indices + 4 * process->curindex;
-
-       /* displists now support array drawing, we treat tri's as fake quad */
-       
-       cur[0] = i1;
-       cur[1] = i2;
-       cur[2] = i3;
-       if (i4 == 0)
-               cur[3] = i3;
-       else 
-               cur[3] = i4;
-       
-       process->curindex++;
-
-}
-
-/* ******************* MEMORY MANAGEMENT *********************** */
-static void *new_pgn_element(PROCESS *process, int size)
-{
-       /* during polygonize 1000s of elements are allocated
-        * and never freed in between. Freeing only done at the end.
-        */
-       int blocksize = 16384;
-       void *adr;
-       
-       if (size > 10000 || size == 0) {
-               printf("incorrect use of new_pgn_element\n");
-       }
-       else if (size == -1) {
-               struct pgn_elements *cur = process->pgn_list.first;
-               while (cur) {
-                       MEM_freeN(cur->data);
-                       cur = cur->next;
-               }
-               BLI_freelistN(&process->pgn_list);
-               
-               return NULL;
-       }
-       
-       size = 4 * ( (size + 3) / 4);
-       
-       if (process->pgn_current) {
-               if (size + process->pgn_offset < blocksize) {
-                       adr = (void *) (process->pgn_current->data + process->pgn_offset);
-                       process->pgn_offset += size;
-                       return adr;
-               }
-       }
-       
-       process->pgn_current = MEM_callocN(sizeof(struct pgn_elements), "newpgn");
-       process->pgn_current->data = MEM_callocN(blocksize, "newpgn");
-       BLI_addtail(&process->pgn_list, process->pgn_current);
-       
-       process->pgn_offset = size;
-       return process->pgn_current->data;
-}
-
-static void freepolygonize(PROCESS *process)
-{
-       MEM_freeN(process->corners);
-       MEM_freeN(process->edges);
-       MEM_freeN(process->centers);
-
-       new_pgn_element(process, -1);
-
-       if (process->vertices.ptr) {
-               MEM_freeN(process->vertices.ptr);
-       }
-}
-
-/**** Cubical Polygonization (optional) ****/
-
-#define LB  0  /* left bottom edge     */
-#define LT  1  /* left top edge        */
-#define LN  2  /* left near edge       */
-#define LF  3  /* left far edge        */
-#define RB  4  /* right bottom edge */
-#define RT  5  /* right top edge       */
-#define RN  6  /* right near edge      */
-#define RF  7  /* right far edge       */
-#define BN  8  /* bottom near edge     */
-#define BF  9  /* bottom far edge      */
-#define TN  10 /* top near edge        */
-#define TF  11 /* top far edge */
-
-static INTLISTS *cubetable[256];
-
-/* edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
-static int corner1[12] = {
-       LBN, LTN, LBN, LBF, RBN, RTN, RBN, RBF, LBN, LBF, LTN, LTF
-};
-static int corner2[12] = {
-       LBF, LTF, LTN, LTF, RBF, RTF, RTN, RTF, RBN, RBF, RTN, RTF
-};
-static int leftface[12] = {
-       B,  L,  L,  F,  R,  T,  N,  R,  N,  B,  T,  F
-};
-/* face on left when going corner1 to corner2 */
-static int rightface[12] = {
-       L,  T,  N,  L,  B,  R,  R,  F,  B,  F,  N,  T
-};
-/* face on right when going corner1 to corner2 */
-
-
-/* docube: triangulate the cube directly, without decomposition */
-
-static void docube(PROCESS *process, CUBE *cube, MetaBall *mb)
-{
-       INTLISTS *polys;
-       CORNER *c1, *c2;
-       int i, index = 0, count, indexar[8];
-       
-       for (i = 0; i < 8; i++) {
-               if (cube->corners[i]->value > 0.0f) {
-                       index += (1 << i);
-               }
-       }
-       
-       for (polys = cubetable[index]; polys; polys = polys->next) {
-               INTLIST *edges;
-               
-               count = 0;
-               
-               for (edges = polys->list; edges; edges = edges->next) {
-                       c1 = cube->corners[corner1[edges->i]];
-                       c2 = cube->corners[corner2[edges->i]];
-                       
-                       indexar[count] = vertid(process, c1, c2, mb);
-                       count++;
-               }
-               if (count > 2) {
-                       switch (count) {
-                               case 3:
-                                       accum_mballfaces(process, indexar[2], indexar[1], indexar[0], 0);
-                                       break;
-                               case 4:
-                                       if (indexar[0] == 0) accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
-                                       else accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
-                                       break;
-                               case 5:
-                                       if (indexar[0] == 0) accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
-                                       else accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
-                               
-                                       accum_mballfaces(process, indexar[4], indexar[3], indexar[0], 0);
-                                       break;
-                               case 6:
-                                       if (indexar[0] == 0) {
-                                               accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
-                                               accum_mballfaces(process, indexar[0], indexar[5], indexar[4], indexar[3]);
-                                       }
-                                       else {
-                                               accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
-                                               accum_mballfaces(process, indexar[5], indexar[4], indexar[3], indexar[0]);
-                                       }
-                                       break;
-                               case 7:
-                                       if (indexar[0] == 0) {
-                                               accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
-                                               accum_mballfaces(process, indexar[0], indexar[5], indexar[4], indexar[3]);
-                                       }
-                                       else {
-                                               accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
-                                               accum_mballfaces(process, indexar[5], indexar[4], indexar[3], indexar[0]);
-                                       }
-                               
-                                       accum_mballfaces(process, indexar[6], indexar[5], indexar[0], 0);
-
-                                       break;
-                       }
-               }
-       }
-}
-
-
-/* testface: given cube at lattice (i, j, k), and four corners of face,
- * if surface crosses face, compute other four corners of adjacent cube
- * and add new cube to cube stack */
-
-static void testface(PROCESS *process, int i, int j, int k, CUBE *old, int bit, int c1, int c2, int c3, int c4)
-{
-       CUBE newc;
-       CUBES *oldcubes = process->cubes;
-       CORNER *corn1, *corn2, *corn3, *corn4;
-       int n, pos;
-
-       corn1 = old->corners[c1];
-       corn2 = old->corners[c2];
-       corn3 = old->corners[c3];
-       corn4 = old->corners[c4];
-       
-       pos = corn1->value > 0.0f ? 1 : 0;
-
-       /* test if no surface crossing */
-       if ( (corn2->value > 0) == pos && (corn3->value > 0) == pos && (corn4->value > 0) == pos) return;
-       /* test if cube out of bounds */
-       /*if ( abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) return;*/
-       /* test if already visited (always as last) */
-       if (setcenter(process, process->centers, i, j, k)) {
-               return;
-       }
-
-       /* create new cube and add cube to top of stack: */
-       process->cubes = (CUBES *) new_pgn_element(process, sizeof(CUBES));
-       process->cubes->next = oldcubes;
-       
-       newc.i = i;
-       newc.j = j;
-       newc.k = k;
-       for (n = 0; n < 8; n++) newc.corners[n] = NULL;
-       
-       newc.corners[FLIP(c1, bit)] = corn1;
-       newc.corners[FLIP(c2, bit)] = corn2;
-       newc.corners[FLIP(c3, bit)] = corn3;
-       newc.corners[FLIP(c4, bit)] = corn4;
-
-       if (newc.corners[0] == NULL) newc.corners[0] = setcorner(process, i, j, k);
-       if (newc.corners[1] == NULL) newc.corners[1] = setcorner(process, i, j, k + 1);
-       if (newc.corners[2] == NULL) newc.corners[2] = setcorner(process, i, j + 1, k);
-       if (newc.corners[3] == NULL) newc.corners[3] = setcorner(process, i, j + 1, k + 1);
-       if (newc.corners[4] == NULL) newc.corners[4] = setcorner(process, i + 1, j, k);
-       if (newc.corners[5] == NULL) newc.corners[5] = setcorner(process, i + 1, j, k + 1);
-       if (newc.corners[6] == NULL) newc.corners[6] = setcorner(process, i + 1, j + 1, k);
-       if (newc.corners[7] == NULL) newc.corners[7] = setcorner(process, i + 1, j + 1, k + 1);
-
-       process->cubes->cube = newc;
-}
-
-/* setcorner: return corner with the given lattice location
- * set (and cache) its function value */
-
-static CORNER *setcorner(PROCESS *process, int i, int j, int k)
-{
-       /* for speed, do corner value caching here */
-       CORNER *c;
-       int index;
-
-       /* does corner exist? */
-       index = HASH(i, j, k);
-       c = process->corners[index];
-       
-       for (; c != NULL; c = c->next) {
-               if (c->i == i && c->j == j && c->k == k) {
-                       return c;
-               }
-       }
-
-       c = (CORNER *) new_pgn_element(process, sizeof(CORNER));
-
-       c->i = i; 
-       c->co[0] = ((float)i - 0.5f) * process->size;
-       c->j = j; 
-       c->co[1] = ((float)j - 0.5f) * process->size;
-       c->k = k; 
-       c->co[2] = ((float)k - 0.5f) * process->size;
-       c->value = process->function(process, c->co[0], c->co[1], c->co[2]);
-       
-       c->next = process->corners[index];
-       process->corners[index] = c;
-       
-       return c;
-}
-
-
-/* nextcwedge: return next clockwise edge from given edge around given face */
-
-static int nextcwedge(int edge, int face)
-{
-       switch (edge) {
-               case LB:
-                       return (face == L) ? LF : BN;
-               case LT:
-                       return (face == L) ? LN : TF;
-               case LN:
-                       return (face == L) ? LB : TN;
-               case LF:
-                       return (face == L) ? LT : BF;
-               case RB:
-                       return (face == R) ? RN : BF;
-               case RT:
-                       return (face == R) ? RF : TN;
-               case RN:
-                       return (face == R) ? RT : BN;
-               case RF:
-                       return (face == R) ? RB : TF;
-               case BN:
-                       return (face == B) ? RB : LN;
-               case BF:
-                       return (face == B) ? LB : RF;
-               case TN:
-                       return (face == T) ? LT : RN;
-               case TF:
-                       return (face == T) ? RT : LF;
-       }
-       return 0;
-}
-
-
-/* otherface: return face adjoining edge that is not the given face */
-
-static int otherface(int edge, int face)
-{
-       int other = leftface[edge];
-       return face == other ? rightface[edge] : other;
-}
-
-
-/* makecubetable: create the 256 entry table for cubical polygonization */
-
-static void makecubetable(void)
-{
-       static bool is_done = false;
-       int i, e, c, done[12], pos[8];
-
-       if (is_done) return;
-       is_done = true;
-
-       for (i = 0; i < 256; i++) {
-               for (e = 0; e < 12; e++) done[e] = 0;
-               for (c = 0; c < 8; c++) pos[c] = MB_BIT(i, c);
-               for (e = 0; e < 12; e++)
-                       if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
-                               INTLIST *ints = NULL;
-                               INTLISTS *lists = (INTLISTS *) MEM_callocN(sizeof(INTLISTS), "mball_intlist");
-                               int start = e, edge = e;
-                               
-                               /* get face that is to right of edge from pos to neg corner: */
-                               int face = pos[corner1[e]] ? rightface[e] : leftface[e];
-                               
-                               while (1) {
-                                       edge = nextcwedge(edge, face);
-                                       done[edge] = 1;
-                                       if (pos[corner1[edge]] != pos[corner2[edge]]) {
-                                               INTLIST *tmp = ints;
-                                               
-                                               ints = (INTLIST *) MEM_callocN(sizeof(INTLIST), "mball_intlist");
-                                               ints->i = edge;
-                                               ints->next = tmp; /* add edge to head of list */
-                                               
-                                               if (edge == start) break;
-                                               face = otherface(edge, face);
-                                       }
-                               }
-                               lists->list = ints; /* add ints to head of table entry */
-                               lists->next = cubetable[i];
-                               cubetable[i] = lists;
-                       }
-       }
-}
-
-void BKE_mball_cubeTable_free(void)
-{
-       int i;
-       INTLISTS *lists, *nlists;
-       INTLIST *ints, *nints;
-
-       for (i = 0; i < 256; i++) {
-               lists = cubetable[i];
-               while (lists) {
-                       nlists = lists->next;
-                       
-                       ints = lists->list;
-                       while (ints) {
-                               nints = ints->next;
-                               MEM_freeN(ints);
-                               ints = nints;
-                       }
-                       
-                       MEM_freeN(lists);
-                       lists = nlists;
-               }
-               cubetable[i] = NULL;
-       }
-}
-
-/**** Storage ****/
-
-/* setcenter: set (i, j, k) entry of table[]
- * return 1 if already set; otherwise, set and return 0 */
-
-static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const int j, const int k)
-{
-       int index;
-       CENTERLIST *newc, *l, *q;
-
-       index = HASH(i, j, k);
-       q = table[index];
-
-       for (l = q; l != NULL; l = l->next) {
-               if (l->i == i && l->j == j && l->k == k) return 1;
-       }
-       
-       newc = (CENTERLIST *) new_pgn_element(process, sizeof(CENTERLIST));
-       newc->i = i; 
-       newc->j = j; 
-       newc->k = k; 
-       newc->next = q;
-       table[index] = newc;
-       
-       return 0;
-}
-
-
-/* setedge: set vertex id for edge */
-
-static void setedge(PROCESS *process,
-                    EDGELIST *table[],
-                    int i1, int j1,
-                    int k1, int i2,
-                    int j2, int k2,
-                    int vid)
-{
-       unsigned int index;
-       EDGELIST *newe;
-       
-       if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
-               int t = i1;
-               i1 = i2;
-               i2 = t;
-               t = j1;
-               j1 = j2;
-               j2 = t;
-               t = k1;
-               k1 = k2;
-               k2 = t;
-       }
-       index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
-       newe = (EDGELIST *) new_pgn_element(process, sizeof(EDGELIST));
-       newe->i1 = i1; 
-       newe->j1 = j1; 
-       newe->k1 = k1;
-       newe->i2 = i2; 
-       newe->j2 = j2; 
-       newe->k2 = k2;
-       newe->vid = vid;
-       newe->next = table[index];
-       table[index] = newe;
-}
-
-
-/* getedge: return vertex id for edge; return -1 if not set */
-
-static int getedge(EDGELIST *table[],
-                   int i1, int j1, int k1,
-                   int i2, int j2, int k2)
-{
-       EDGELIST *q;
-       
-       if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
-               int t = i1;
-               i1 = i2;
-               i2 = t;
-               t = j1;
-               j1 = j2;
-               j2 = t;
-               t = k1;
-               k1 = k2;
-               k2 = t;
-       }
-       q = table[HASH(i1, j1, k1) + HASH(i2, j2, k2)];
-       for (; q != NULL; q = q->next) {
-               if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 &&
-                   q->i2 == i2 && q->j2 == j2 && q->k2 == k2)
-               {
-                       return q->vid;
-               }
-       }
-       return -1;
-}
-
-
-/**** Vertices ****/
-
-#undef R
-
-
-
-/* vertid: return index for vertex on edge:
- * c1->value and c2->value are presumed of different sign
- * return saved index if any; else compute vertex and save */
-
-/* addtovertices: add v to sequence of vertices */
-
-static void addtovertices(VERTICES *vertices, VERTEX v)
-{
-       if (vertices->count == vertices->max) {
-               int i;
-               VERTEX *newv;
-               vertices->max = vertices->count == 0 ? 10 : 2 * vertices->count;
-               newv = (VERTEX *) MEM_callocN(vertices->max * sizeof(VERTEX), "addtovertices");
-               
-               for (i = 0; i < vertices->count; i++) newv[i] = vertices->ptr[i];
-               
-               if (vertices->ptr != NULL) MEM_freeN(vertices->ptr);
-               vertices->ptr = newv;
-       }
-       vertices->ptr[vertices->count++] = v;
-}
-
-/* vnormal: compute unit length surface normal at point */
-
-static void vnormal(PROCESS *process, const float point[3], float r_no[3])
-{
-       const float delta = 0.2f * process->delta;
-       const float f = process->function(process, point[0], point[1], point[2]);
-
-       r_no[0] = process->function(process, point[0] + delta, point[1], point[2]) - f;
-       r_no[1] = process->function(process, point[0], point[1] + delta, point[2]) - f;
-       r_no[2] = process->function(process, point[0], point[1], point[2] + delta) - f;
-
-#if 1
-       normalize_v3(r_no);
-#else
-       f = normalize_v3(r_no);
-       
-       if (0) {
-               float tvec[3];
-               
-               delta *= 2.0f;
-               
-               f = process->function(process, point[0], point[1], point[2]);
-       
-               tvec[0] = process->function(process, point[0] + delta, point[1], point[2]) - f;
-               tvec[1] = process->function(process, point[0], point[1] + delta, point[2]) - f;
-               tvec[2] = process->function(process, point[0], point[1], point[2] + delta) - f;
-       
-               if (normalize_v3(tvec) != 0.0f) {
-                       add_v3_v3(r_no, tvec);
-                       normalize_v3(r_no);
-               }
-       }
-#endif
-}
-
-
-static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2, MetaBall *mb)
-{
-       VERTEX v;
-       int vid = getedge(process->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
-
-       if (vid != -1) {
-               return vid;  /* previously computed */
-       }
-
-       converge(process, c1->co, c2->co, c1->value, c2->value, v.co, mb, 1); /* position */
-       vnormal(process, v.co, v.no);
-
-       addtovertices(&process->vertices, v);            /* save vertex */
-       vid = process->vertices.count - 1;
-       setedge(process, process->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
-       
-       return vid;
-}
-
-
-/* converge: from two points of differing sign, converge to zero crossing */
-/* watch it: p1 and p2 are used to calculate */
-static void converge(PROCESS *process, const float p1[3], const float p2[3], float v1, float v2,
-                     float p[3], MetaBall *mb, int f)
-{
-       int i = 0;
-       float pos[3], neg[3];
-       float positive = 0.0f, negative = 0.0f;
-       float dvec[3];
-       
-       if (v1 < 0) {
-               copy_v3_v3(pos, p2);
-               copy_v3_v3(neg, p1);
-               positive = v2;
-               negative = v1;
-       }
-       else {
-               copy_v3_v3(pos, p1);
-               copy_v3_v3(neg, p2);
-               positive = v1;
-               negative = v2;
-       }
-
-       sub_v3_v3v3(dvec, pos, neg);
-
-/* Approximation by linear interpolation is faster then binary subdivision,
- * but it results sometimes (mb->thresh < 0.2) into the strange results */
-       if ((mb->thresh > 0.2f) && (f == 1)) {
-               if ((dvec[1] == 0.0f) && (dvec[2] == 0.0f)) {
-                       p[0] = neg[0] - negative * dvec[0] / (positive - negative);
-                       p[1] = neg[1];
-                       p[2] = neg[2];
-                       return;
-               }
-               if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
-                       p[0] = neg[0];
-                       p[1] = neg[1] - negative * dvec[1] / (positive - negative);
-                       p[2] = neg[2];
-                       return;
-               }
-               if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
-                       p[0] = neg[0];
-                       p[1] = neg[1];
-                       p[2] = neg[2] - negative * dvec[2] / (positive - negative);
-                       return;
-               }
-       }
-
-       if ((dvec[1] == 0.0f) && (dvec[2] == 0.0f)) {
-               p[1] = neg[1];
-               p[2] = neg[2];
-               while (1) {
-                       if (i++ == RES) return;
-                       p[0] = 0.5f * (pos[0] + neg[0]);
-                       if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[0] = p[0];
-                       else                                                       neg[0] = p[0];
-               }
-       }
-
-       if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
-               p[0] = neg[0];
-               p[2] = neg[2];
-               while (1) {
-                       if (i++ == RES) return;
-                       p[1] = 0.5f * (pos[1] + neg[1]);
-                       if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[1] = p[1];
-                       else                                                       neg[1] = p[1];
-               }
-       }
-
-       if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
-               p[0] = neg[0];
-               p[1] = neg[1];
-               while (1) {
-                       if (i++ == RES) return;
-                       p[2] = 0.5f * (pos[2] + neg[2]);
-                       if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[2] = p[2];
-                       else                                                       neg[2] = p[2];
-               }
-       }
-
-       /* This is necessary to find start point */
-       while (1) {
-               mid_v3_v3v3(&p[0], pos, neg);
-
-               if (i++ == RES) {
-                       return;
-               }
-
-               if ((process->function(process, p[0], p[1], p[2])) > 0.0f) {
-                       copy_v3_v3(pos, &p[0]);
-               }
-               else {
-                       copy_v3_v3(neg, &p[0]);
-               }
-       }
-}
-
-/* ************************************** */
-static void add_cube(PROCESS *process, int i, int j, int k, int count)
-{
-       CUBES *ncube;
-       int n;
-       int a, b, c;
-
-       /* hmmm, not only one, but eight cube will be added on the stack 
-        * ... */
-       for (a = i - 1; a < i + count; a++)
-               for (b = j - 1; b < j + count; b++)
-                       for (c = k - 1; c < k + count; c++) {
-                               /* test if cube has been found before */
-                               if (setcenter(process, process->centers, a, b, c) == 0) {
-                                       /* push cube on stack: */
-                                       ncube = (CUBES *) new_pgn_element(process, sizeof(CUBES));
-                                       ncube->next = process->cubes;
-                                       process->cubes = ncube;
-
-                                       ncube->cube.i = a;
-                                       ncube->cube.j = b;
-                                       ncube->cube.k = c;
-
-                                       /* set corners of initial cube: */
-                                       for (n = 0; n < 8; n++)
-                                               ncube->cube.corners[n] = setcorner(process, a + MB_BIT(n, 2), b + MB_BIT(n, 1), c + MB_BIT(n, 0));
-                               }
-                       }
-}
-
-
-static void find_first_points(PROCESS *process, MetaBall *mb, int a)
-{
-       MetaElem *ml;
-       float f;
-
-       ml = process->mainb[a];
-       f = 1.0f - (mb->thresh / ml->s);
-
-       /* Skip, when Stiffness of MetaElement is too small ... MetaElement can't be
-        * visible alone ... but still can influence others MetaElements :-) */
-       if (f > 0.0f) {
-               float IN[3] = {0.0f}, OUT[3] = {0.0f}, in[3] = {0.0f}, out[3];
-               int i, j, k, c_i, c_j, c_k;
-               int index[3] = {1, 0, -1};
-               float in_v /*, out_v*/;
-               float workp[3];
-               float dvec[3];
-               float tmp_v, workp_v, max_len_sq, nx, ny, nz, max_dim;
-
-               calc_mballco(ml, in);
-               in_v = process->function(process, in[0], in[1], in[2]);
-
-               for (i = 0; i < 3; i++) {
-                       switch (ml->type) {
-                               case MB_BALL:
-                                       OUT[0] = out[0] = IN[0] + index[i] * ml->rad;
-                                       break;
-                               case MB_TUBE:
-                               case MB_PLANE:
-                               case MB_ELIPSOID:
-                               case MB_CUBE:
-                                       OUT[0] = out[0] = IN[0] + index[i] * (ml->expx + ml->rad);
-                                       break;
-                       }
-
-                       for (j = 0; j < 3; j++) {
-                               switch (ml->type) {
-                                       case MB_BALL:
-                                               OUT[1] = out[1] = IN[1] + index[j] * ml->rad;
-                                               break;
-                                       case MB_TUBE:
-                                       case MB_PLANE:
-                                       case MB_ELIPSOID:
-                                       case MB_CUBE:
-                                               OUT[1] = out[1] = IN[1] + index[j] * (ml->expy + ml->rad);
-                                               break;
-                               }
-                       
-                               for (k = 0; k < 3; k++) {
-                                       out[0] = OUT[0];
-                                       out[1] = OUT[1];
-                                       switch (ml->type) {
-                                               case MB_BALL:
-                                               case MB_TUBE:
-                                               case MB_PLANE:
-                                                       out[2] = IN[2] + index[k] * ml->rad;
-                                                       break;
-                                               case MB_ELIPSOID:
-                                               case MB_CUBE:
-                                                       out[2] = IN[2] + index[k] * (ml->expz + ml->rad);
-                                                       break;
-                                       }
-
-                                       calc_mballco(ml, out);
-
-                                       /*out_v = process->function(out[0], out[1], out[2]);*/ /*UNUSED*/
-
-                                       /* find "first points" on Implicit Surface of MetaElemnt ml */
-                                       copy_v3_v3(workp, in);
-                                       workp_v = in_v;
-                                       max_len_sq = len_squared_v3v3(out, in);
-
-                                       nx = fabsf((out[0] - in[0]) / process->size);
-                                       ny = fabsf((out[1] - in[1]) / process->size);
-                                       nz = fabsf((out[2] - in[2]) / process->size);
-                                       
-                                       max_dim = max_fff(nx, ny, nz);
-                                       if (max_dim != 0.0f) {
-                                               float len_sq = 0.0f;
-
-                                               dvec[0] = (out[0] - in[0]) / max_dim;
-                                               dvec[1] = (out[1] - in[1]) / max_dim;
-                                               dvec[2] = (out[2] - in[2]) / max_dim;
-
-                                               while (len_sq <= max_len_sq) {
-                                                       add_v3_v3(workp, dvec);
-
-                                                       /* compute value of implicite function */
-                                                       tmp_v = process->function(process, workp[0], workp[1], workp[2]);
-                                                       /* add cube to the stack, when value of implicite function crosses zero value */
-                                                       if ((tmp_v < 0.0f && workp_v >= 0.0f) || (tmp_v > 0.0f && workp_v <= 0.0f)) {
-
-                                                               /* indexes of CUBE, which includes "first point" */
-                                                               c_i = (int)floor(workp[0] / process->size);
-                                                               c_j = (int)floor(workp[1] / process->size);
-                                                               c_k = (int)floor(workp[2] / process->size);
-                                                               
-                                                               /* add CUBE (with indexes c_i, c_j, c_k) to the stack,
-                                                                * this cube includes found point of Implicit Surface */
-                                                               if ((ml->flag & MB_NEGATIVE) == 0) {
-                                                                       add_cube(process, c_i, c_j, c_k, 1);
-                                                               }
-                                                               else {
-                                                                       add_cube(process, c_i, c_j, c_k, 2);
-                                                               }
-                                                       }
-                                                       len_sq = len_squared_v3v3(workp, in);
-                                                       workp_v = tmp_v;
-
-                                               }
-                                       }
-                               }
-                       }
-               }
-       }
-}
-
-static void polygonize(PROCESS *process, MetaBall *mb)
-{
-       CUBE c;
-       int a;
-
-       process->vertices.count = process->vertices.max = 0;
-       process->vertices.ptr = NULL;
-
-       /* allocate hash tables and build cube polygon table: */
-       process->centers = MEM_callocN(HASHSIZE * sizeof(CENTERLIST *), "mbproc->centers");
-       process->corners = MEM_callocN(HASHSIZE * sizeof(CORNER *), "mbproc->corners");
-       process->edges = MEM_callocN(2 * HASHSIZE * sizeof(EDGELIST *), "mbproc->edges");
-       makecubetable();
-
-       for (a = 0; a < process->totelem; a++) {
-
-               /* try to find 8 points on the surface for each MetaElem */
-               find_first_points(process, mb, a);
-       }
-
-       /* polygonize all MetaElems of current MetaBall */
-       while (process->cubes != NULL) { /* process active cubes till none left */
-               c = process->cubes->cube;
-
-               /* polygonize the cube directly: */
-               docube(process, &c, mb);
-               
-               /* pop current cube from stack */
-               process->cubes = process->cubes->next;
-               
-               /* test six face directions, maybe add to stack: */
-               testface(process, c.i - 1, c.j, c.k, &c, 2, LBN, LBF, LTN, LTF);
-               testface(process, c.i + 1, c.j, c.k, &c, 2, RBN, RBF, RTN, RTF);
-               testface(process, c.i, c.j - 1, c.k, &c, 1, LBN, LBF, RBN, RBF);
-               testface(process, c.i, c.j + 1, c.k, &c, 1, LTN, LTF, RTN, RTF);
-               testface(process, c.i, c.j, c.k - 1, &c, 0, LBN, LTN, RBN, RTN);
-               testface(process, c.i, c.j, c.k + 1, &c, 0, LBF, LTF, RBF, RTF);
-       }
-}
-
-static float init_meta(EvaluationContext *eval_ctx, PROCESS *process, Scene *scene, Object *ob)    /* return totsize */
-{
-       Scene *sce_iter = scene;
-       Base *base;
-       Object *bob;
-       MetaBall *mb;
-       MetaElem *ml;
-       float size, totsize, obinv[4][4], obmat[4][4], vec[3];
-       //float max = 0.0f;
-       int a, obnr, zero_size = 0;
-       char obname[MAX_ID_NAME];
-       SceneBaseIter iter;
-
-       copy_m4_m4(obmat, ob->obmat);   /* to cope with duplicators from BKE_scene_base_iter_next */
-       invert_m4_m4(obinv, ob->obmat);
-       a = 0;
-       
-       BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
-       
-       /* make main array */
-       BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 0, NULL, NULL);
-       while (BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 1, &base, &bob)) {
-
-               if (bob->type == OB_MBALL) {
-                       zero_size = 0;
-                       ml = NULL;
-
-                       if (bob == ob && (base->flag & OB_FROMDUPLI) == 0) {
-                               mb = ob->data;
-       
-                               if (mb->editelems) ml = mb->editelems->first;
-                               else ml = mb->elems.first;
-                       }
-                       else {
-                               char name[MAX_ID_NAME];
-                               int nr;
-                               
-                               BLI_split_name_num(name, &nr, bob->id.name + 2, '.');
-                               if (STREQ(obname, name)) {
-                                       mb = bob->data;
-                                       
-                                       if (mb->editelems) ml = mb->editelems->first;
-                                       else ml = mb->elems.first;
-                               }
-                       }
-
-                       /* when metaball object has zero scale, then MetaElem to this MetaBall
-                        * will not be put to mainb array */
-                       if (has_zero_axis_m4(bob->obmat)) {
-                               zero_size = 1;
-                       }
-                       else if (bob->parent) {
-                               struct Object *pob = bob->parent;
-                               while (pob) {
-                                       if (has_zero_axis_m4(pob->obmat)) {
-                                               zero_size = 1;
-                                               break;
-                                       }
-                                       pob = pob->parent;
-                               }
-                       }
-
-                       if (zero_size) {
-                               unsigned int ml_count = 0;
-                               while (ml) {
-                                       ml_count++;
-                                       ml = ml->next;
-                               }
-                               process->totelem -= ml_count;
-                       }
-                       else {
-                               while (ml) {
-                                       if (!(ml->flag & MB_HIDE)) {
-                                               int i;
-                                               float temp1[4][4], temp2[4][4], temp3[4][4];
-                                               float (*mat)[4] = NULL, (*imat)[4] = NULL;
-                                               float max_x, max_y, max_z, min_x, min_y, min_z;
-                                               float expx, expy, expz;
-
-                                               max_x = max_y = max_z = -3.4e38;
-                                               min_x = min_y = min_z =  3.4e38;
-
-                                               /* too big stiffness seems only ugly due to linear interpolation
-                                                * no need to have possibility for too big stiffness */
-                                               if (ml->s > 10.0f) ml->s = 10.0f;
-
-                                               /* Rotation of MetaElem is stored in quat */
-                                               quat_to_mat4(temp3, ml->quat);
-
-                                               /* Translation of MetaElem */
-                                               unit_m4(temp2);
-                                               temp2[3][0] = ml->x;
-                                               temp2[3][1] = ml->y;
-                                               temp2[3][2] = ml->z;
-
-                                               mul_m4_m4m4(temp1, temp2, temp3);
-
-                                               /* make a copy because of duplicates */
-                                               process->mainb[a] = new_pgn_element(process, sizeof(MetaElem));
-                                               *(process->mainb[a]) = *ml;
-                                               process->mainb[a]->bb = new_pgn_element(process, sizeof(BoundBox));
-
-                                               mat = new_pgn_element(process, 4 * 4 * sizeof(float));
-                                               imat = new_pgn_element(process, 4 * 4 * sizeof(float));
-
-                                               /* mat is the matrix to transform from mball into the basis-mball */
-                                               invert_m4_m4(obinv, obmat);
-                                               mul_m4_m4m4(temp2, obinv, bob->obmat);
-                                               /* MetaBall transformation */
-                                               mul_m4_m4m4(mat, temp2, temp1);
-
-                                               invert_m4_m4(imat, mat);
-
-                                               process->mainb[a]->rad2 = ml->rad * ml->rad;
-
-                                               process->mainb[a]->mat = (float *) mat;
-                                               process->mainb[a]->imat = (float *) imat;
-
-                                               if (!MB_TYPE_SIZE_SQUARED(ml->type)) {
-                                                       expx = ml->expx;
-                                                       expy = ml->expy;
-                                                       expz = ml->expz;
-                                               }
-                                               else {
-                                                       expx = ml->expx * ml->expx;
-                                                       expy = ml->expy * ml->expy;
-                                                       expz = ml->expz * ml->expz;
-                                               }
-
-                                               /* untransformed Bounding Box of MetaElem */
-                                               /* TODO, its possible the elem type has been changed and the exp* values can use a fallback */
-                                               copy_v3_fl3(process->mainb[a]->bb->vec[0], -expx, -expy, -expz);  /* 0 */
-                                               copy_v3_fl3(process->mainb[a]->bb->vec[1], +expx, -expy, -expz);  /* 1 */
-                                               copy_v3_fl3(process->mainb[a]->bb->vec[2], +expx, +expy, -expz);  /* 2 */
-                                               copy_v3_fl3(process->mainb[a]->bb->vec[3], -expx, +expy, -expz);  /* 3 */
-                                               copy_v3_fl3(process->mainb[a]->bb->vec[4], -expx, -expy, +expz);  /* 4 */
-                                               copy_v3_fl3(process->mainb[a]->bb->vec[5], +expx, -expy, +expz);  /* 5 */
-                                               copy_v3_fl3(process->mainb[a]->bb->vec[6], +expx, +expy, +expz);  /* 6 */
-                                               copy_v3_fl3(process->mainb[a]->bb->vec[7], -expx, +expy, +expz);  /* 7 */
-
-                                               /* transformation of Metalem bb */
-                                               for (i = 0; i < 8; i++)
-                                                       mul_m4_v3((float (*)[4])mat, process->mainb[a]->bb->vec[i]);
-
-                                               /* find max and min of transformed bb */
-                                               for (i = 0; i < 8; i++) {
-                                                       /* find maximums */
-                                                       if (process->mainb[a]->bb->vec[i][0] > max_x) max_x = process->mainb[a]->bb->vec[i][0];
-                                                       if (process->mainb[a]->bb->vec[i][1] > max_y) max_y = process->mainb[a]->bb->vec[i][1];
-                                                       if (process->mainb[a]->bb->vec[i][2] > max_z) max_z = process->mainb[a]->bb->vec[i][2];
-                                                       /* find  minimums */
-                                                       if (process->mainb[a]->bb->vec[i][0] < min_x) min_x = process->mainb[a]->bb->vec[i][0];
-                                                       if (process->mainb[a]->bb->vec[i][1] < min_y) min_y = process->mainb[a]->bb->vec[i][1];
-                                                       if (process->mainb[a]->bb->vec[i][2] < min_z) min_z = process->mainb[a]->bb->vec[i][2];
-                                               }
-                                       
-                                               /* create "new" bb, only point 0 and 6, which are
-                                                * necessary for octal tree filling */
-                                               process->mainb[a]->bb->vec[0][0] = min_x - ml->rad;
-                                               process->mainb[a]->bb->vec[0][1] = min_y - ml->rad;
-                                               process->mainb[a]->bb->vec[0][2] = min_z - ml->rad;
-
-                                               process->mainb[a]->bb->vec[6][0] = max_x + ml->rad;
-                                               process->mainb[a]->bb->vec[6][1] = max_y + ml->rad;
-                                               process->mainb[a]->bb->vec[6][2] = max_z + ml->rad;
-
-                                               a++;
-                                       }
-                                       ml = ml->next;
-                               }
-                       }
-               }
-       }
-
-       
-       /* totsize (= 'manhattan' radius) */
-       totsize = 0.0;
-       for (a = 0; a < process->totelem; a++) {
-               
-               vec[0] = process->mainb[a]->x + process->mainb[a]->rad + process->mainb[a]->expx;
-               vec[1] = process->mainb[a]->y + process->mainb[a]->rad + process->mainb[a]->expy;
-               vec[2] = process->mainb[a]->z + process->mainb[a]->rad + process->mainb[a]->expz;
-
-               calc_mballco(process->mainb[a], vec);
-       
-               size = fabsf(vec[0]);
-               if (size > totsize) totsize = size;
-               size = fabsf(vec[1]);
-               if (size > totsize) totsize = size;
-               size = fabsf(vec[2]);
-               if (size > totsize) totsize = size;
-
-               vec[0] = process->mainb[a]->x - process->mainb[a]->rad;
-               vec[1] = process->mainb[a]->y - process->mainb[a]->rad;
-               vec[2] = process->mainb[a]->z - process->mainb[a]->rad;
-                               
-               calc_mballco(process->mainb[a], vec);
-       
-               size = fabsf(vec[0]);
-               if (size > totsize) totsize = size;
-               size = fabsf(vec[1]);
-               if (size > totsize) totsize = size;
-               size = fabsf(vec[2]);
-               if (size > totsize) totsize = size;
-       }
-
-       for (a = 0; a < process->totelem; a++) {
-               process->thresh += densfunc(process->mainb[a], 2.0f * totsize, 2.0f * totsize, 2.0f * totsize);
-       }
-
-       return totsize;
-}
-
-/* if MetaElem lies in node, then node includes MetaElem pointer (ml_p)
- * pointing at MetaElem (ml)
- */
-static void fill_metaball_octal_node(octal_node *node, MetaElem *ml, short i)
-{
-       ml_pointer *ml_p;
-
-       ml_p = MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
-       ml_p->ml = ml;
-       BLI_addtail(&(node->nodes[i]->elems), ml_p);
-       node->count++;
-       
-       if ((ml->flag & MB_NEGATIVE) == 0) {
-               node->nodes[i]->pos++;
-       }
-       else {
-               node->nodes[i]->neg++;
-       }
-}
-
-/* Node is subdivided as is illustrated on the following figure:
- * 
- *      +------+------+
- *     /      /      /|
- *    +------+------+ |
- *   /      /      /| +
- *  +------+------+ |/|
- *  |      |      | + |
- *  |      |      |/| +
- *  +------+------+ |/
- *  |      |      | +
- *  |      |      |/
- *  +------+------+
- *  
- */
-static void subdivide_metaball_octal_node(octal_node *node, float size_x, float size_y, float size_z, short depth)
-{
-       MetaElem *ml;
-       ml_pointer *ml_p;
-       float x, y, z;
-       int a, i;
-
-       /* create new nodes */
-       for (a = 0; a < 8; a++) {
-               node->nodes[a] = MEM_mallocN(sizeof(octal_node), "octal_node");
-               for (i = 0; i < 8; i++)
-                       node->nodes[a]->nodes[i] = NULL;
-               node->nodes[a]->parent = node;
-               BLI_listbase_clear(&node->nodes[a]->elems);
-               node->nodes[a]->count = 0;
-               node->nodes[a]->neg = 0;
-               node->nodes[a]->pos = 0;
-       }
-
-       size_x /= 2;
-       size_y /= 2;
-       size_z /= 2;
-       
-       /* center of node */
-       node->x = x = node->x_min + size_x;
-       node->y = y = node->y_min + size_y;
-       node->z = z = node->z_min + size_z;
-
-       /* setting up of border points of new nodes */
-       node->nodes[0]->x_min = node->x_min;
-       node->nodes[0]->y_min = node->y_min;
-       node->nodes[0]->z_min = node->z_min;
-       node->nodes[0]->x = node->nodes[0]->x_min + size_x / 2;
-       node->nodes[0]->y = node->nodes[0]->y_min + size_y / 2;
-       node->nodes[0]->z = node->nodes[0]->z_min + size_z / 2;
-       
-       node->nodes[1]->x_min = x;
-       node->nodes[1]->y_min = node->y_min;
-       node->nodes[1]->z_min = node->z_min;
-       node->nodes[1]->x = node->nodes[1]->x_min + size_x / 2;
-       node->nodes[1]->y = node->nodes[1]->y_min + size_y / 2;
-       node->nodes[1]->z = node->nodes[1]->z_min + size_z / 2;
-
-       node->nodes[2]->x_min = x;
-       node->nodes[2]->y_min = y;
-       node->nodes[2]->z_min = node->z_min;
-       node->nodes[2]->x = node->nodes[2]->x_min + size_x / 2;
-       node->nodes[2]->y = node->nodes[2]->y_min + size_y / 2;
-       node->nodes[2]->z = node->nodes[2]->z_min + size_z / 2;
-
-       node->nodes[3]->x_min = node->x_min;
-       node->nodes[3]->y_min = y;
-       node->nodes[3]->z_min = node->z_min;
-       node->nodes[3]->x = node->nodes[3]->x_min + size_x / 2;
-       node->nodes[3]->y = node->nodes[3]->y_min + size_y / 2;
-       node->nodes[3]->z = node->nodes[3]->z_min + size_z / 2;
-
-       node->nodes[4]->x_min = node->x_min;
-       node->nodes[4]->y_min = node->y_min;
-       node->nodes[4]->z_min = z;
-       node->nodes[4]->x = node->nodes[4]->x_min + size_x / 2;
-       node->nodes[4]->y = node->nodes[4]->y_min + size_y / 2;
-       node->nodes[4]->z = node->nodes[4]->z_min + size_z / 2;
-       
-       node->nodes[5]->x_min = x;
-       node->nodes[5]->y_min = node->y_min;
-       node->nodes[5]->z_min = z;
-       node->nodes[5]->x = node->nodes[5]->x_min + size_x / 2;
-       node->nodes[5]->y = node->nodes[5]->y_min + size_y / 2;
-       node->nodes[5]->z = node->nodes[5]->z_min + size_z / 2;
-
-       node->nodes[6]->x_min = x;
-       node->nodes[6]->y_min = y;
-       node->nodes[6]->z_min = z;
-       node->nodes[6]->x = node->nodes[6]->x_min + size_x / 2;
-       node->nodes[6]->y = node->nodes[6]->y_min + size_y / 2;
-       node->nodes[6]->z = node->nodes[6]->z_min + size_z / 2;
-
-       node->nodes[7]->x_min = node->x_min;
-       node->nodes[7]->y_min = y;
-       node->nodes[7]->z_min = z;
-       node->nodes[7]->x = node->nodes[7]->x_min + size_x / 2;
-       node->nodes[7]->y = node->nodes[7]->y_min + size_y / 2;
-       node->nodes[7]->z = node->nodes[7]->z_min + size_z / 2;
-
-       ml_p = node->elems.first;
-       
-       /* setting up references of MetaElems for new nodes */
-       while (ml_p) {
-               ml = ml_p->ml;
-               if (ml->bb->vec[0][2] < z) {
-                       if (ml->bb->vec[0][1] < y) {
-                               /* vec[0][0] lies in first octant */
-                               if (ml->bb->vec[0][0] < x) {
-                                       /* ml belongs to the (0)1st node */
-                                       fill_metaball_octal_node(node, ml, 0);
-
-                                       /* ml belongs to the (3)4th node */
-                                       if (ml->bb->vec[6][1] >= y) {
-                                               fill_metaball_octal_node(node, ml, 3);
-
-                                               /* ml belongs to the (7)8th node */
-                                               if (ml->bb->vec[6][2] >= z) {
-                                                       fill_metaball_octal_node(node, ml, 7);
-                                               }
-                                       }
-       
-                                       /* ml belongs to the (1)2nd node */
-                                       if (ml->bb->vec[6][0] >= x) {
-                                               fill_metaball_octal_node(node, ml, 1);
-
-                                               /* ml belongs to the (5)6th node */
-                                               if (ml->bb->vec[6][2] >= z) {
-                                                       fill_metaball_octal_node(node, ml, 5);
-                                               }
-                                       }
-
-                                       /* ml belongs to the (2)3th node */
-                                       if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
-                                               fill_metaball_octal_node(node, ml, 2);
-                                               
-                                               /* ml belong to the (6)7th node */
-                                               if (ml->bb->vec[6][2] >= z) {
-                                                       fill_metaball_octal_node(node, ml, 6);
-                                               }
-                                               
-                                       }
-                       
-                                       /* ml belongs to the (4)5th node too */
-                                       if (ml->bb->vec[6][2] >= z) {
-                                               fill_metaball_octal_node(node, ml, 4);
-                                       }
-
-                                       
-                                       
-                               }
-                               /* vec[0][0] is in the (1)second octant */
-                               else {
-                                       /* ml belong to the (1)2nd node */
-                                       fill_metaball_octal_node(node, ml, 1);
-
-                                       /* ml belongs to the (2)3th node */
-                                       if (ml->bb->vec[6][1] >= y) {
-                                               fill_metaball_octal_node(node, ml, 2);
-
-                                               /* ml belongs to the (6)7th node */
-                                               if (ml->bb->vec[6][2] >= z) {
-                                                       fill_metaball_octal_node(node, ml, 6);
-                                               }
-                                               
-                                       }
-                                       
-                                       /* ml belongs to the (5)6th node */
-                                       if (ml->bb->vec[6][2] >= z) {
-                                               fill_metaball_octal_node(node, ml, 5);
-                                       }
-                               }
-                       }
-                       else {
-                               /* vec[0][0] is in the (3)4th octant */
-                               if (ml->bb->vec[0][0] < x) {
-                                       /* ml belongs to the (3)4nd node */
-                                       fill_metaball_octal_node(node, ml, 3);
-                                       
-                                       /* ml belongs to the (7)8th node */
-                                       if (ml->bb->vec[6][2] >= z) {
-                                               fill_metaball_octal_node(node, ml, 7);
-                                       }
-                               
-
-                                       /* ml belongs to the (2)3th node */
-                                       if (ml->bb->vec[6][0] >= x) {
-                                               fill_metaball_octal_node(node, ml, 2);
-                                       
-                                               /* ml belongs to the (6)7th node */
-                                               if (ml->bb->vec[6][2] >= z) {
-                                                       fill_metaball_octal_node(node, ml, 6);
-                                               }
-                                       }
-                               }
-
-                       }
-
-                       /* vec[0][0] is in the (2)3th octant */
-                       if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
-                               /* ml belongs to the (2)3th node */
-                               fill_metaball_octal_node(node, ml, 2);
-                               
-                               /* ml belongs to the (6)7th node */
-                               if (ml->bb->vec[6][2] >= z) {
-                                       fill_metaball_octal_node(node, ml, 6);
-                               }
-                       }
-               }
-               else {
-                       if (ml->bb->vec[0][1] < y) {
-                               /* vec[0][0] lies in (4)5th octant */
-                               if (ml->bb->vec[0][0] < x) {
-                                       /* ml belongs to the (4)5th node */
-                                       fill_metaball_octal_node(node, ml, 4);
-
-                                       if (ml->bb->vec[6][0] >= x) {
-                                               fill_metaball_octal_node(node, ml, 5);
-                                       }
-
-                                       if (ml->bb->vec[6][1] >= y) {
-                                               fill_metaball_octal_node(node, ml, 7);
-                                       }
-                                       
-                                       if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
-                                               fill_metaball_octal_node(node, ml, 6);
-                                       }
-                               }
-                               /* vec[0][0] lies in (5)6th octant */
-                               else {
-                                       fill_metaball_octal_node(node, ml, 5);
-
-                                       if (ml->bb->vec[6][1] >= y) {
-                                               fill_metaball_octal_node(node, ml, 6);
-                                       }
-                               }
-                       }
-                       else {
-                               /* vec[0][0] lies in (7)8th octant */
-                               if (ml->bb->vec[0][0] < x) {
-                                       fill_metaball_octal_node(node, ml, 7);
-
-                                       if (ml->bb->vec[6][0] >= x) {
-                                               fill_metaball_octal_node(node, ml, 6);
-                                       }
-                               }
-
-                       }
-                       
-                       /* vec[0][0] lies in (6)7th octant */
-                       if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
-                               fill_metaball_octal_node(node, ml, 6);
-                       }
-               }
-               ml_p = ml_p->next;
-       }
-
-       /* free references of MetaElems for curent node (it is not needed anymore) */
-       BLI_freelistN(&node->elems);
-
-       depth--;
-       
-       if (depth > 0) {
-               for (a = 0; a < 8; a++) {
-                       if (node->nodes[a]->count > 0) /* if node is not empty, then it is subdivided */
-                               subdivide_metaball_octal_node(node->nodes[a], size_x, size_y, size_z, depth);
-               }
-       }
-}
-
-/* free all octal nodes recursively */
-static void free_metaball_octal_node(octal_node *node)
-{
-       int a;
-       for (a = 0; a < 8; a++) {
-               if (node->nodes[a] != NULL) free_metaball_octal_node(node->nodes[a]);
-       }
-       BLI_freelistN(&node->elems);
-       MEM_freeN(node);
-}
-
-/* If scene include more than one MetaElem, then octree is used */
-static void init_metaball_octal_tree(PROCESS *process, int depth)
-{
-       struct octal_node *node;
-       ml_pointer *ml_p;
-       float size[3];
-       int a;
-       
-       process->metaball_tree = MEM_mallocN(sizeof(octal_tree), "metaball_octal_tree");
-       process->metaball_tree->first = node = MEM_mallocN(sizeof(octal_node), "metaball_octal_node");
-       /* maximal depth of octree */
-       process->metaball_tree->depth = depth;
-
-       process->metaball_tree->neg = node->neg = 0;
-       process->metaball_tree->pos = node->pos = 0;
-       
-       BLI_listbase_clear(&node->elems);
-       node->count = 0;
-
-       for (a = 0; a < 8; a++)
-               node->nodes[a] = NULL;
-
-       node->x_min = node->y_min = node->z_min = FLT_MAX;
-       node->x_max = node->y_max = node->z_max = -FLT_MAX;
-
-       /* size of octal tree scene */
-       for (a = 0; a < process->totelem; a++) {
-               if (process->mainb[a]->bb->vec[0][0] < node->x_min) node->x_min = process->mainb[a]->bb->vec[0][0];
-               if (process->mainb[a]->bb->vec[0][1] < node->y_min) node->y_min = process->mainb[a]->bb->vec[0][1];
-               if (process->mainb[a]->bb->vec[0][2] < node->z_min) node->z_min = process->mainb[a]->bb->vec[0][2];
-
-               if (process->mainb[a]->bb->vec[6][0] > node->x_max) node->x_max = process->mainb[a]->bb->vec[6][0];
-               if (process->mainb[a]->bb->vec[6][1] > node->y_max) node->y_max = process->mainb[a]->bb->vec[6][1];
-               if (process->mainb[a]->bb->vec[6][2] > node->z_max) node->z_max = process->mainb[a]->bb->vec[6][2];
-
-               ml_p = MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
-               ml_p->ml = process->mainb[a];
-               BLI_addtail(&node->elems, ml_p);
-
-               if ((process->mainb[a]->flag & MB_NEGATIVE) == 0) {
-                       /* number of positive MetaElem in scene */
-                       process->metaball_tree->pos++;
-               }
-               else {
-                       /* number of negative MetaElem in scene */
-                       process->metaball_tree->neg++;
-               }
-       }
-
-       /* size of first node */
-       size[0] = node->x_max - node->x_min;
-       size[1] = node->y_max - node->y_min;
-       size[2] = node->z_max - node->z_min;
-
-       /* first node is subdivided recursively */
-       subdivide_metaball_octal_node(node, size[0], size[1], size[2], process->metaball_tree->depth);
-}
-
-static void mball_count(EvaluationContext *eval_ctx, PROCESS *process, Scene *scene, Object *basis)
-{
-       Scene *sce_iter = scene;
-       Base *base;
-       Object *ob, *bob = basis;
-       MetaElem *ml = NULL;
-       int basisnr, obnr;
-       char basisname[MAX_ID_NAME], obname[MAX_ID_NAME];
-       SceneBaseIter iter;
-
-       BLI_split_name_num(basisname, &basisnr, basis->id.name + 2, '.');
-       process->totelem = 0;
-
-       BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 0, NULL, NULL);
-       while (BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 1, &base, &ob)) {
-               if (ob->type == OB_MBALL) {
-                       if (ob == bob) {
-                               MetaBall *mb = ob->data;
-
-                               /* if bob object is in edit mode, then dynamic list of all MetaElems
-                                * is stored in editelems */
-                               if (mb->editelems) ml = mb->editelems->first;
-                               /* if bob object is in object mode */
-                               else ml = mb->elems.first;
-                       }
-                       else {
-                               BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
-
-                               /* object ob has to be in same "group" ... it means, that it has to have
-                                * same base of its name */
-                               if (STREQ(obname, basisname)) {
-                                       MetaBall *mb = ob->data;
-
-                                       /* if object is in edit mode, then dynamic list of all MetaElems
-                                        * is stored in editelems */
-                                       if (mb->editelems) ml = mb->editelems->first;
-                                       /* if bob object is in object mode */
-                                       else ml = mb->elems.first;
-                               }
-                       }
-
-                       for ( ; ml; ml = ml->next) {
-                               if (!(ml->flag & MB_HIDE)) {
-                                       process->totelem++;
-                               }
-                       }
-               }
-       }
-}
-
-void BKE_mball_polygonize(EvaluationContext *eval_ctx, Scene *scene, Object *ob, ListBase *dispbase)
-{
-       MetaBall *mb;
-       DispList *dl;
-       int a, nr_cubes;
-       float *co, *no, totsize, width;
-       PROCESS process = {0};
-
-       mb = ob->data;
-
-       mball_count(eval_ctx, &process, scene, ob);
-
-       if (process.totelem == 0) return;
-       if ((eval_ctx->mode != DAG_EVAL_RENDER) && (mb->flag == MB_UPDATE_NEVER)) return;
-       if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_FAST) return;
-
-       process.thresh = mb->thresh;
-
-       /* total number of MetaElems (totelem) is precomputed in find_basis_mball() function */
-       process.mainb = MEM_mallocN(sizeof(void *) * process.totelem, "mainb");
-       
-       /* initialize all mainb (MetaElems) */
-       totsize = init_meta(eval_ctx, &process, scene, ob);
-
-       /* if scene includes more than one MetaElem, then octal tree optimization is used */
-       if ((process.totelem >    1) && (process.totelem <=   64)) init_metaball_octal_tree(&process, 1);
-       if ((process.totelem >   64) && (process.totelem <=  128)) init_metaball_octal_tree(&process, 2);
-       if ((process.totelem >  128) && (process.totelem <=  512)) init_metaball_octal_tree(&process, 3);
-       if ((process.totelem >  512) && (process.totelem <= 1024)) init_metaball_octal_tree(&process, 4);
-       if (process.totelem  > 1024)                               init_metaball_octal_tree(&process, 5);
-
-       /* don't polygonize metaballs with too high resolution (base mball to small)
-        * note: Eps was 0.0001f but this was giving problems for blood animation for durian, using 0.00001f */
-       if (process.metaball_tree) {
-               if (ob->size[0] <= 0.00001f * (process.metaball_tree->first->x_max - process.metaball_tree->first->x_min) ||
-                   ob->size[1] <= 0.00001f * (process.metaball_tree->first->y_max - process.metaball_tree->first->y_min) ||
-                   ob->size[2] <= 0.00001f * (process.metaball_tree->first->z_max - process.metaball_tree->first->z_min))
-               {
-                       new_pgn_element(&process, -1); /* free values created by init_meta */
-
-                       MEM_freeN(process.mainb);
-
-                       /* free tree */
-                       free_metaball_octal_node(process.metaball_tree->first);
-                       MEM_freeN(process.metaball_tree);
-
-                       return;
-               }
-       }
-
-       /* width is size per polygonize cube */
-       if (eval_ctx->mode == DAG_EVAL_RENDER) {
-               width = mb->rendersize;
-       }
-       else {
-               width = mb->wiresize;
-               if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_HALFRES) {
-                       width *= 2;
-               }
-       }
-       /* nr_cubes is just for safety, minimum is totsize */
-       nr_cubes = (int)(0.5f + totsize / width);
-
-       /* init process */
-       process.function = metaball;
-       process.size = width;
-       process.bounds = nr_cubes;
-       process.cubes = NULL;
-       process.delta = width / (float)(RES * RES);
-
-       polygonize(&process, mb);
-       
-       MEM_freeN(process.mainb);
-
-       /* free octal tree */
-       if (process.totelem > 1) {
-               free_metaball_octal_node(process.metaball_tree->first);
-               MEM_freeN(process.metaball_tree);
-               process.metaball_tree = NULL;
-       }
-
-       if (process.curindex) {
-               VERTEX *ptr = process.vertices.ptr;
-
-               dl = MEM_callocN(sizeof(DispList), "mbaldisp");
-               BLI_addtail(dispbase, dl);
-               dl->type = DL_INDEX4;
-               dl->nr = process.vertices.count;
-               dl->parts = process.curindex;
-
-               dl->index = process.indices;
-               process.indices = NULL;
-
-               a = process.vertices.count;
-               dl->verts = co = MEM_mallocN(sizeof(float[3]) * a, "mballverts");
-               dl->nors = no = MEM_mallocN(sizeof(float[3]) * a, "mballnors");
-
-               for (a = 0; a < process.vertices.count; ptr++, a++, no += 3, co += 3) {
-                       copy_v3_v3(co, ptr->co);
-                       copy_v3_v3(no, ptr->no);
-               }
-       }
-
-       freepolygonize(&process);
-}
-
 bool BKE_mball_minmax_ex(MetaBall *mb, float min[3], float max[3],
                          float obmat[4][4], const short flag)
 {
diff --git a/source/blender/blenkernel/intern/mball_tessellate.c b/source/blender/blenkernel/intern/mball_tessellate.c
new file mode 100644 (file)
index 0000000..b1bef13
--- /dev/null
@@ -0,0 +1,1984 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * Contributor(s): Jiri Hnidek <jiri.hnidek@vslib.cz>.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+/** \file blender/blenkernel/intern/mball_tessellate.c
+ *  \ingroup bke
+ */
+
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+#include <stdlib.h>
+#include <ctype.h>
+#include <float.h>
+#include "MEM_guardedalloc.h"
+
+#include "DNA_material_types.h"
+#include "DNA_object_types.h"
+#include "DNA_meta_types.h"
+#include "DNA_scene_types.h"
+
+
+#include "BLI_blenlib.h"
+#include "BLI_math.h"
+#include "BLI_utildefines.h"
+
+#include "BKE_global.h"
+#include "BKE_main.h"
+
+/*  #include "BKE_object.h" */
+#include "BKE_animsys.h"
+#include "BKE_curve.h"
+#include "BKE_depsgraph.h"
+#include "BKE_scene.h"
+#include "BKE_library.h"
+#include "BKE_displist.h"
+#include "BKE_mball.h"
+#include "BKE_object.h"
+#include "BKE_material.h"
+
+/* Data types */
+
+typedef struct vertex {         /* surface vertex */
+       float co[3];  /* position and surface normal */
+       float no[3];
+} VERTEX;
+
+typedef struct vertices {       /* list of vertices in polygonization */
+       int count, max;             /* # vertices, max # allowed */
+       VERTEX *ptr;                /* dynamically allocated */
+} VERTICES;
+
+typedef struct corner {         /* corner of a cube */
+       int i, j, k;                /* (i, j, k) is index within lattice */
+       float co[3], value;       /* location and function value */
+       struct corner *next;
+} CORNER;
+
+typedef struct cube {           /* partitioning cell (cube) */
+       int i, j, k;                /* lattice location of cube */
+       CORNER *corners[8];         /* eight corners */
+} CUBE;
+
+typedef struct cubes {          /* linked list of cubes acting as stack */
+       CUBE cube;                  /* a single cube */
+       struct cubes *next;         /* remaining elements */
+} CUBES;
+
+typedef struct centerlist {     /* list of cube locations */
+       int i, j, k;                /* cube location */
+       struct centerlist *next;    /* remaining elements */
+} CENTERLIST;
+
+typedef struct edgelist {       /* list of edges */
+       int i1, j1, k1, i2, j2, k2; /* edge corner ids */
+       int vid;                    /* vertex id */
+       struct edgelist *next;      /* remaining elements */
+} EDGELIST;
+
+typedef struct intlist {        /* list of integers */
+       int i;                      /* an integer */
+       struct intlist *next;       /* remaining elements */
+} INTLIST;
+
+typedef struct intlists {       /* list of list of integers */
+       INTLIST *list;              /* a list of integers */
+       struct intlists *next;      /* remaining elements */
+} INTLISTS;
+
+/* dividing scene using octal tree makes polygonisation faster */
+typedef struct ml_pointer {
+       struct ml_pointer *next, *prev;
+       struct MetaElem *ml;
+} ml_pointer;
+
+typedef struct octal_node {
+       struct octal_node *nodes[8];/* children of current node */
+       struct octal_node *parent;  /* parent of current node */
+       struct ListBase elems;      /* ListBase of MetaElem pointers (ml_pointer) */
+       float x_min, y_min, z_min;  /* 1st border point */
+       float x_max, y_max, z_max;  /* 7th border point */
+       float x, y, z;              /* center of node */
+       int pos, neg;               /* number of positive and negative MetaElements in the node */
+       int count;                  /* number of MetaElems, which belongs to the node */
+} octal_node;
+
+typedef struct octal_tree {
+       struct octal_node *first;   /* first node */
+       int pos, neg;               /* number of positive and negative MetaElements in the scene */
+       short depth;                /* number of scene subdivision */
+} octal_tree;
+
+struct pgn_elements {
+       struct pgn_elements *next, *prev;
+       char *data;
+};
+
+typedef struct process {        /* parameters, function, storage */
+       /* ** old G_mb contents ** */
+       float thresh;
+       int totelem;
+       MetaElem **mainb;
+       octal_tree *metaball_tree;
+
+       /* ** old process contents ** */
+
+       /* what happens here? floats, I think. */
+       /*  float (*function)(void);     */     /* implicit surface function */
+       float (*function)(struct process *, float, float, float);
+       float size, delta;          /* cube size, normal delta */
+       int bounds;                 /* cube range within lattice */
+       CUBES *cubes;               /* active cubes */
+       VERTICES vertices;          /* surface vertices */
+       CENTERLIST **centers;       /* cube center hash table */
+       CORNER **corners;           /* corner value hash table */
+       EDGELIST **edges;           /* edge and vertex id hash table */
+
+       /* Runtime things */
+       int *indices;
+       int totindex, curindex;
+
+       int pgn_offset;
+       struct pgn_elements *pgn_current;
+       ListBase pgn_list;
+} PROCESS;
+
+/* Forward declarations */
+static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2, MetaBall *mb);
+static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const int j, const int k);
+static CORNER *setcorner(PROCESS *process, int i, int j, int k);
+static void converge(PROCESS *process, const float p1[3], const float p2[3], float v1, float v2,
+                     float p[3], MetaBall *mb, int f);
+
+
+/* ******************** ARITH ************************* */
+
+/* BASED AT CODE (but mostly rewritten) :
+ * C code from the article
+ * "An Implicit Surface Polygonizer"
+ * by Jules Bloomenthal, jbloom@beauty.gmu.edu
+ * in "Graphics Gems IV", Academic Press, 1994
+ *
+ * Authored by Jules Bloomenthal, Xerox PARC.
+ * Copyright (c) Xerox Corporation, 1991.  All rights reserved.
+ * Permission is granted to reproduce, use and distribute this code for
+ * any and all purposes, provided that this notice appears in all copies. */
+
+#define RES 12 /* # converge iterations    */
+
+#define L   0  /* left direction:      -x, -i */
+#define R   1  /* right direction:     +x, +i */
+#define B   2  /* bottom direction: -y, -j */
+#define T   3  /* top direction:       +y, +j */
+#define N   4  /* near direction:      -z, -k */
+#define F   5  /* far direction:       +z, +k */
+#define LBN 0  /* left bottom near corner  */
+#define LBF 1  /* left bottom far corner   */
+#define LTN 2  /* left top near corner     */
+#define LTF 3  /* left top far corner      */
+#define RBN 4  /* right bottom near corner */
+#define RBF 5  /* right bottom far corner  */
+#define RTN 6  /* right top near corner    */
+#define RTF 7  /* right top far corner     */
+
+/* the LBN corner of cube (i, j, k), corresponds with location
+ * (i-0.5)*size, (j-0.5)*size, (k-0.5)*size) */
+
+#define HASHBIT     (5)
+#define HASHSIZE    (size_t)(1 << (3 * HASHBIT))   /*! < hash table size (32768) */
+
+#define HASH(i, j, k) ((((( (i) & 31) << 5) | ( (j) & 31)) << 5) | ( (k) & 31) )
+
+#define MB_BIT(i, bit) (((i) >> (bit)) & 1)
+#define FLIP(i, bit) ((i) ^ 1 << (bit)) /* flip the given bit of i */
+
+
+/* **************** POLYGONIZATION ************************ */
+
+static void calc_mballco(MetaElem *ml, float vec[3])
+{
+       if (ml->mat) {
+               mul_m4_v3((float (*)[4])ml->mat, vec);
+       }
+}
+
+static float densfunc(MetaElem *ball, float x, float y, float z)
+{
+       float dist2;
+       float dvec[3] = {x, y, z};
+
+       mul_m4_v3((float (*)[4])ball->imat, dvec);
+
+       switch (ball->type) {
+               case MB_BALL:
+                       /* do nothing */
+                       break;
+               case MB_TUBE:
+                       if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
+                       else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
+                       else                            dvec[0] = 0.0;
+                       break;
+               case MB_PLANE:
+                       if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
+                       else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
+                       else                            dvec[0] = 0.0;
+                       if      (dvec[1] >  ball->expy) dvec[1] -= ball->expy;
+                       else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
+                       else                            dvec[1] = 0.0;
+                       break;
+               case MB_ELIPSOID:
+                       dvec[0] /= ball->expx;
+                       dvec[1] /= ball->expy;
+                       dvec[2] /= ball->expz;
+                       break;
+               case MB_CUBE:
+                       if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
+                       else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
+                       else                            dvec[0] = 0.0;
+
+                       if      (dvec[1] >  ball->expy) dvec[1] -= ball->expy;
+                       else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
+                       else                            dvec[1] = 0.0;
+
+                       if      (dvec[2] >  ball->expz) dvec[2] -= ball->expz;
+                       else if (dvec[2] < -ball->expz) dvec[2] += ball->expz;
+                       else                            dvec[2] = 0.0;
+                       break;
+
+               /* *** deprecated, could be removed?, do-versioned at least *** */
+               case MB_TUBEX:
+                       if      (dvec[0] >  ball->len) dvec[0] -= ball->len;
+                       else if (dvec[0] < -ball->len) dvec[0] += ball->len;
+                       else                           dvec[0] = 0.0;
+                       break;
+               case MB_TUBEY:
+                       if      (dvec[1] >  ball->len) dvec[1] -= ball->len;
+                       else if (dvec[1] < -ball->len) dvec[1] += ball->len;
+                       else                           dvec[1] = 0.0;
+                       break;
+               case MB_TUBEZ:
+                       if      (dvec[2] >  ball->len) dvec[2] -= ball->len;
+                       else if (dvec[2] < -ball->len) dvec[2] += ball->len;
+                       else                           dvec[2] = 0.0;
+                       break;
+               /* *** end deprecated *** */
+       }
+
+       dist2 = 1.0f - (len_squared_v3(dvec) / ball->rad2);
+
+       if ((ball->flag & MB_NEGATIVE) == 0) {
+               return (dist2 < 0.0f) ? -0.5f : (ball->s * dist2 * dist2 * dist2) - 0.5f;
+       }
+       else {
+               return (dist2 < 0.0f) ? 0.5f : 0.5f - (ball->s * dist2 * dist2 * dist2);
+       }
+}
+
+static octal_node *find_metaball_octal_node(octal_node *node, float x, float y, float z, short depth)
+{
+       if (!depth) return node;
+       
+       if (z < node->z) {
+               if (y < node->y) {
+                       if (x < node->x) {
+                               if (node->nodes[0])
+                                       return find_metaball_octal_node(node->nodes[0], x, y, z, depth--);
+                               else
+                                       return node;
+                       }
+                       else {
+                               if (node->nodes[1])
+                                       return find_metaball_octal_node(node->nodes[1], x, y, z, depth--);
+                               else
+                                       return node;
+                       }
+               }
+               else {
+                       if (x < node->x) {
+                               if (node->nodes[3])
+                                       return find_metaball_octal_node(node->nodes[3], x, y, z, depth--);
+                               else
+                                       return node;
+                       }
+                       else {
+                               if (node->nodes[2])
+                                       return find_metaball_octal_node(node->nodes[2], x, y, z, depth--);
+                               else
+                                       return node;
+                       }
+               }
+       }
+       else {
+               if (y < node->y) {
+                       if (x < node->x) {
+                               if (node->nodes[4])
+                                       return find_metaball_octal_node(node->nodes[4], x, y, z, depth--);
+                               else
+                                       return node;
+                       }
+                       else {
+                               if (node->nodes[5])
+                                       return find_metaball_octal_node(node->nodes[5], x, y, z, depth--);
+                               else
+                                       return node;
+                       }
+               }
+               else {
+                       if (x < node->x) {
+                               if (node->nodes[7])
+                                       return find_metaball_octal_node(node->nodes[7], x, y, z, depth--);
+                               else
+                                       return node;
+                       }
+                       else {
+                               if (node->nodes[6])
+                                       return find_metaball_octal_node(node->nodes[6], x, y, z, depth--);
+                               else
+                                       return node;
+                       }
+               }
+       }
+       
+       /* all cases accounted for */
+       BLI_assert(0);
+}
+
+static float metaball(PROCESS *process, float x, float y, float z)
+/*  float x, y, z; */
+{
+       octal_tree *metaball_tree = process->metaball_tree;
+       struct octal_node *node;
+       struct ml_pointer *ml_p;
+       float dens = 0;
+       int a;
+       
+       if (process->totelem > 1) {
+               node = find_metaball_octal_node(metaball_tree->first, x, y, z, metaball_tree->depth);
+               if (node) {
+                       for (ml_p = node->elems.first; ml_p; ml_p = ml_p->next) {
+                               dens += densfunc(ml_p->ml, x, y, z);
+                       }
+
+                       dens += -0.5f * (metaball_tree->pos - node->pos);
+                       dens +=  0.5f * (metaball_tree->neg - node->neg);
+               }
+               else {
+                       for (a = 0; a < process->totelem; a++) {
+                               dens += densfunc(process->mainb[a], x, y, z);
+                       }
+               }
+       }
+       else {
+               dens += densfunc(process->mainb[0], x, y, z);
+       }
+
+       return process->thresh - dens;
+}
+
+/* ******************************************** */
+
+static void accum_mballfaces(PROCESS *process, int i1, int i2, int i3, int i4)
+{
+       int *newi, *cur;
+       /* static int i = 0; I would like to delete altogether, but I don't dare to, yet */
+
+       if (process->totindex == process->curindex) {
+               process->totindex += 256;
+               newi = MEM_mallocN(4 * sizeof(int) * process->totindex, "vertindex");
+               
+               if (process->indices) {
+                       memcpy(newi, process->indices, 4 * sizeof(int) * (process->totindex - 256));
+                       MEM_freeN(process->indices);
+               }
+               process->indices = newi;
+       }
+       
+       cur = process->indices + 4 * process->curindex;
+
+       /* displists now support array drawing, we treat tri's as fake quad */
+       
+       cur[0] = i1;
+       cur[1] = i2;
+       cur[2] = i3;
+       if (i4 == 0)
+               cur[3] = i3;
+       else 
+               cur[3] = i4;
+       
+       process->curindex++;
+
+}
+
+/* ******************* MEMORY MANAGEMENT *********************** */
+static void *new_pgn_element(PROCESS *process, int size)
+{
+       /* during polygonize 1000s of elements are allocated
+        * and never freed in between. Freeing only done at the end.
+        */
+       int blocksize = 16384;
+       void *adr;
+       
+       if (size > 10000 || size == 0) {
+               printf("incorrect use of new_pgn_element\n");
+       }
+       else if (size == -1) {
+               struct pgn_elements *cur = process->pgn_list.first;
+               while (cur) {
+                       MEM_freeN(cur->data);
+                       cur = cur->next;
+               }
+               BLI_freelistN(&process->pgn_list);
+               
+               return NULL;
+       }
+       
+       size = 4 * ( (size + 3) / 4);
+       
+       if (process->pgn_current) {
+               if (size + process->pgn_offset < blocksize) {
+                       adr = (void *) (process->pgn_current->data + process->pgn_offset);
+                       process->pgn_offset += size;
+                       return adr;
+               }
+       }
+       
+       process->pgn_current = MEM_callocN(sizeof(struct pgn_elements), "newpgn");
+       process->pgn_current->data = MEM_callocN(blocksize, "newpgn");
+       BLI_addtail(&process->pgn_list, process->pgn_current);
+       
+       process->pgn_offset = size;
+       return process->pgn_current->data;
+}
+
+static void freepolygonize(PROCESS *process)
+{
+       MEM_freeN(process->corners);
+       MEM_freeN(process->edges);
+       MEM_freeN(process->centers);
+
+       new_pgn_element(process, -1);
+
+       if (process->vertices.ptr) {
+               MEM_freeN(process->vertices.ptr);
+       }
+}
+
+/**** Cubical Polygonization (optional) ****/
+
+#define LB  0  /* left bottom edge     */
+#define LT  1  /* left top edge        */
+#define LN  2  /* left near edge       */
+#define LF  3  /* left far edge        */
+#define RB  4  /* right bottom edge */
+#define RT  5  /* right top edge       */
+#define RN  6  /* right near edge      */
+#define RF  7  /* right far edge       */
+#define BN  8  /* bottom near edge     */
+#define BF  9  /* bottom far edge      */
+#define TN  10 /* top near edge        */
+#define TF  11 /* top far edge */
+
+static INTLISTS *cubetable[256];
+
+/* edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
+static int corner1[12] = {
+       LBN, LTN, LBN, LBF, RBN, RTN, RBN, RBF, LBN, LBF, LTN, LTF
+};
+static int corner2[12] = {
+       LBF, LTF, LTN, LTF, RBF, RTF, RTN, RTF, RBN, RBF, RTN, RTF
+};
+static int leftface[12] = {
+       B,  L,  L,  F,  R,  T,  N,  R,  N,  B,  T,  F
+};
+/* face on left when going corner1 to corner2 */
+static int rightface[12] = {
+       L,  T,  N,  L,  B,  R,  R,  F,  B,  F,  N,  T
+};
+/* face on right when going corner1 to corner2 */
+
+
+/* docube: triangulate the cube directly, without decomposition */
+
+static void docube(PROCESS *process, CUBE *cube, MetaBall *mb)
+{
+       INTLISTS *polys;
+       CORNER *c1, *c2;
+       int i, index = 0, count, indexar[8];
+       
+       for (i = 0; i < 8; i++) {
+               if (cube->corners[i]->value > 0.0f) {
+                       index += (1 << i);
+               }
+       }
+       
+       for (polys = cubetable[index]; polys; polys = polys->next) {
+               INTLIST *edges;
+               
+               count = 0;
+               
+               for (edges = polys->list; edges; edges = edges->next) {
+                       c1 = cube->corners[corner1[edges->i]];
+                       c2 = cube->corners[corner2[edges->i]];
+                       
+                       indexar[count] = vertid(process, c1, c2, mb);
+                       count++;
+               }
+               if (count > 2) {
+                       switch (count) {
+                               case 3:
+                                       accum_mballfaces(process, indexar[2], indexar[1], indexar[0], 0);
+                                       break;
+                               case 4:
+                                       if (indexar[0] == 0) accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
+                                       else accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
+                                       break;
+                               case 5:
+                                       if (indexar[0] == 0) accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
+                                       else accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
+                               
+                                       accum_mballfaces(process, indexar[4], indexar[3], indexar[0], 0);
+                                       break;
+                               case 6:
+                                       if (indexar[0] == 0) {
+                                               accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
+                                               accum_mballfaces(process, indexar[0], indexar[5], indexar[4], indexar[3]);
+                                       }
+                                       else {
+                                               accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
+                                               accum_mballfaces(process, indexar[5], indexar[4], indexar[3], indexar[0]);
+                                       }
+                                       break;
+                               case 7:
+                                       if (indexar[0] == 0) {
+                                               accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
+                                               accum_mballfaces(process, indexar[0], indexar[5], indexar[4], indexar[3]);
+                                       }
+                                       else {
+                                               accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
+                                               accum_mballfaces(process, indexar[5], indexar[4], indexar[3], indexar[0]);
+                                       }
+                               
+                                       accum_mballfaces(process, indexar[6], indexar[5], indexar[0], 0);
+
+                                       break;
+                       }
+               }
+       }
+}
+
+
+/* testface: given cube at lattice (i, j, k), and four corners of face,
+ * if surface crosses face, compute other four corners of adjacent cube
+ * and add new cube to cube stack */
+
+static void testface(PROCESS *process, int i, int j, int k, CUBE *old, int bit, int c1, int c2, int c3, int c4)
+{
+       CUBE newc;
+       CUBES *oldcubes = process->cubes;
+       CORNER *corn1, *corn2, *corn3, *corn4;
+       int n, pos;
+
+       corn1 = old->corners[c1];
+       corn2 = old->corners[c2];
+       corn3 = old->corners[c3];
+       corn4 = old->corners[c4];
+       
+       pos = corn1->value > 0.0f ? 1 : 0;
+
+       /* test if no surface crossing */
+       if ( (corn2->value > 0) == pos && (corn3->value > 0) == pos && (corn4->value > 0) == pos) return;
+       /* test if cube out of bounds */
+       /*if ( abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) return;*/
+       /* test if already visited (always as last) */
+       if (setcenter(process, process->centers, i, j, k)) {
+               return;
+       }
+
+       /* create new cube and add cube to top of stack: */
+       process->cubes = (CUBES *) new_pgn_element(process, sizeof(CUBES));
+       process->cubes->next = oldcubes;
+       
+       newc.i = i;
+       newc.j = j;
+       newc.k = k;
+       for (n = 0; n < 8; n++) newc.corners[n] = NULL;
+       
+       newc.corners[FLIP(c1, bit)] = corn1;
+       newc.corners[FLIP(c2, bit)] = corn2;
+       newc.corners[FLIP(c3, bit)] = corn3;
+       newc.corners[FLIP(c4, bit)] = corn4;
+
+       if (newc.corners[0] == NULL) newc.corners[0] = setcorner(process, i, j, k);
+       if (newc.corners[1] == NULL) newc.corners[1] = setcorner(process, i, j, k + 1);
+       if (newc.corners[2] == NULL) newc.corners[2] = setcorner(process, i, j + 1, k);
+       if (newc.corners[3] == NULL) newc.corners[3] = setcorner(process, i, j + 1, k + 1);
+       if (newc.corners[4] == NULL) newc.corners[4] = setcorner(process, i + 1, j, k);
+       if (newc.corners[5] == NULL) newc.corners[5] = setcorner(process, i + 1, j, k + 1);
+       if (newc.corners[6] == NULL) newc.corners[6] = setcorner(process, i + 1, j + 1, k);
+       if (newc.corners[7] == NULL) newc.corners[7] = setcorner(process, i + 1, j + 1, k + 1);
+
+       process->cubes->cube = newc;
+}
+
+/* setcorner: return corner with the given lattice location
+ * set (and cache) its function value */
+
+static CORNER *setcorner(PROCESS *process, int i, int j, int k)
+{
+       /* for speed, do corner value caching here */
+       CORNER *c;
+       int index;
+
+       /* does corner exist? */
+       index = HASH(i, j, k);
+       c = process->corners[index];
+       
+       for (; c != NULL; c = c->next) {
+               if (c->i == i && c->j == j && c->k == k) {
+                       return c;
+               }
+       }
+
+       c = (CORNER *) new_pgn_element(process, sizeof(CORNER));
+
+       c->i = i; 
+       c->co[0] = ((float)i - 0.5f) * process->size;
+       c->j = j; 
+       c->co[1] = ((float)j - 0.5f) * process->size;
+       c->k = k; 
+       c->co[2] = ((float)k - 0.5f) * process->size;
+       c->value = process->function(process, c->co[0], c->co[1], c->co[2]);
+       
+       c->next = process->corners[index];
+       process->corners[index] = c;
+       
+       return c;
+}
+
+
+/* nextcwedge: return next clockwise edge from given edge around given face */
+
+static int nextcwedge(int edge, int face)
+{
+       switch (edge) {
+               case LB:
+                       return (face == L) ? LF : BN;
+               case LT:
+                       return (face == L) ? LN : TF;
+               case LN:
+                       return (face == L) ? LB : TN;
+               case LF:
+                       return (face == L) ? LT : BF;
+               case RB:
+                       return (face == R) ? RN : BF;
+               case RT:
+                       return (face == R) ? RF : TN;
+               case RN:
+                       return (face == R) ? RT : BN;
+               case RF:
+                       return (face == R) ? RB : TF;
+               case BN:
+                       return (face == B) ? RB : LN;
+               case BF:
+                       return (face == B) ? LB : RF;
+               case TN:
+                       return (face == T) ? LT : RN;
+               case TF:
+                       return (face == T) ? RT : LF;
+       }
+       return 0;
+}
+
+
+/* otherface: return face adjoining edge that is not the given face */
+
+static int otherface(int edge, int face)
+{
+       int other = leftface[edge];
+       return face == other ? rightface[edge] : other;
+}
+
+
+/* makecubetable: create the 256 entry table for cubical polygonization */
+
+static void makecubetable(void)
+{
+       static bool is_done = false;
+       int i, e, c, done[12], pos[8];
+
+       if (is_done) return;
+       is_done = true;
+
+       for (i = 0; i < 256; i++) {
+               for (e = 0; e < 12; e++) done[e] = 0;
+               for (c = 0; c < 8; c++) pos[c] = MB_BIT(i, c);
+               for (e = 0; e < 12; e++)
+                       if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
+                               INTLIST *ints = NULL;
+                               INTLISTS *lists = (INTLISTS *) MEM_callocN(sizeof(INTLISTS), "mball_intlist");
+                               int start = e, edge = e;
+                               
+                               /* get face that is to right of edge from pos to neg corner: */
+                               int face = pos[corner1[e]] ? rightface[e] : leftface[e];
+                               
+                               while (1) {
+                                       edge = nextcwedge(edge, face);
+                                       done[edge] = 1;
+                                       if (pos[corner1[edge]] != pos[corner2[edge]]) {
+                                               INTLIST *tmp = ints;
+                                               
+                                               ints = (INTLIST *) MEM_callocN(sizeof(INTLIST), "mball_intlist");
+                                               ints->i = edge;
+                                               ints->next = tmp; /* add edge to head of list */
+                                               
+                                               if (edge == start) break;
+                                               face = otherface(edge, face);
+                                       }
+                               }
+                               lists->list = ints; /* add ints to head of table entry */
+                               lists->next = cubetable[i];
+                               cubetable[i] = lists;
+                       }
+       }
+}
+
+void BKE_mball_cubeTable_free(void)
+{
+       int i;
+       INTLISTS *lists, *nlists;
+       INTLIST *ints, *nints;
+
+       for (i = 0; i < 256; i++) {
+               lists = cubetable[i];
+               while (lists) {
+                       nlists = lists->next;
+                       
+                       ints = lists->list;
+                       while (ints) {
+                               nints = ints->next;
+                               MEM_freeN(ints);
+                               ints = nints;
+                       }
+                       
+                       MEM_freeN(lists);
+                       lists = nlists;
+               }
+               cubetable[i] = NULL;
+       }
+}
+
+/**** Storage ****/
+
+/* setcenter: set (i, j, k) entry of table[]
+ * return 1 if already set; otherwise, set and return 0 */
+
+static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const int j, const int k)
+{
+       int index;
+       CENTERLIST *newc, *l, *q;
+
+       index = HASH(i, j, k);
+       q = table[index];
+
+       for (l = q; l != NULL; l = l->next) {
+               if (l->i == i && l->j == j && l->k == k) return 1;
+       }
+       
+       newc = (CENTERLIST *) new_pgn_element(process, sizeof(CENTERLIST));
+       newc->i = i; 
+       newc->j = j; 
+       newc->k = k; 
+       newc->next = q;
+       table[index] = newc;
+       
+       return 0;
+}
+
+
+/* setedge: set vertex id for edge */
+
+static void setedge(PROCESS *process,
+                    EDGELIST *table[],
+                    int i1, int j1,
+                    int k1, int i2,
+                    int j2, int k2,
+                    int vid)
+{
+       unsigned int index;
+       EDGELIST *newe;
+       
+       if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
+               int t = i1;
+               i1 = i2;
+               i2 = t;
+               t = j1;
+               j1 = j2;
+               j2 = t;
+               t = k1;
+               k1 = k2;
+               k2 = t;
+       }
+       index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
+       newe = (EDGELIST *) new_pgn_element(process, sizeof(EDGELIST));
+       newe->i1 = i1; 
+       newe->j1 = j1; 
+       newe->k1 = k1;
+       newe->i2 = i2; 
+       newe->j2 = j2; 
+       newe->k2 = k2;
+       newe->vid = vid;
+       newe->next = table[index];
+       table[index] = newe;
+}
+
+
+/* getedge: return vertex id for edge; return -1 if not set */
+
+static int getedge(EDGELIST *table[],
+                   int i1, int j1, int k1,
+                   int i2, int j2, int k2)
+{
+       EDGELIST *q;
+       
+       if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
+               int t = i1;
+               i1 = i2;
+               i2 = t;
+               t = j1;
+               j1 = j2;
+               j2 = t;
+               t = k1;
+               k1 = k2;
+               k2 = t;
+       }
+       q = table[HASH(i1, j1, k1) + HASH(i2, j2, k2)];
+       for (; q != NULL; q = q->next) {
+               if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 &&
+                   q->i2 == i2 && q->j2 == j2 && q->k2 == k2)
+               {
+                       return q->vid;
+               }
+       }
+       return -1;
+}
+
+
+/**** Vertices ****/
+
+#undef R
+
+
+
+/* vertid: return index for vertex on edge:
+ * c1->value and c2->value are presumed of different sign
+ * return saved index if any; else compute vertex and save */
+
+/* addtovertices: add v to sequence of vertices */
+
+static void addtovertices(VERTICES *vertices, VERTEX v)
+{
+       if (vertices->count == vertices->max) {
+               int i;
+               VERTEX *newv;
+               vertices->max = vertices->count == 0 ? 10 : 2 * vertices->count;
+               newv = (VERTEX *) MEM_callocN(vertices->max * sizeof(VERTEX), "addtovertices");
+               
+               for (i = 0; i < vertices->count; i++) newv[i] = vertices->ptr[i];
+               
+               if (vertices->ptr != NULL) MEM_freeN(vertices->ptr);
+               vertices->ptr = newv;
+       }
+       vertices->ptr[vertices->count++] = v;
+}
+
+/* vnormal: compute unit length surface normal at point */
+
+static void vnormal(PROCESS *process, const float point[3], float r_no[3])
+{
+       const float delta = 0.2f * process->delta;
+       const float f = process->function(process, point[0], point[1], point[2]);
+
+       r_no[0] = process->function(process, point[0] + delta, point[1], point[2]) - f;
+       r_no[1] = process->function(process, point[0], point[1] + delta, point[2]) - f;
+       r_no[2] = process->function(process, point[0], point[1], point[2] + delta) - f;
+
+#if 1
+       normalize_v3(r_no);
+#else
+       f = normalize_v3(r_no);
+       
+       if (0) {
+               float tvec[3];
+               
+               delta *= 2.0f;
+               
+               f = process->function(process, point[0], point[1], point[2]);
+       
+               tvec[0] = process->function(process, point[0] + delta, point[1], point[2]) - f;
+               tvec[1] = process->function(process, point[0], point[1] + delta, point[2]) - f;
+               tvec[2] = process->function(process, point[0], point[1], point[2] + delta) - f;
+       
+               if (normalize_v3(tvec) != 0.0f) {
+                       add_v3_v3(r_no, tvec);
+                       normalize_v3(r_no);
+               }
+       }
+#endif
+}
+
+
+static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2, MetaBall *mb)
+{
+       VERTEX v;
+       int vid = getedge(process->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
+
+       if (vid != -1) {
+               return vid;  /* previously computed */
+       }
+
+       converge(process, c1->co, c2->co, c1->value, c2->value, v.co, mb, 1); /* position */
+       vnormal(process, v.co, v.no);
+
+       addtovertices(&process->vertices, v);            /* save vertex */
+       vid = process->vertices.count - 1;
+       setedge(process, process->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
+       
+       return vid;
+}
+
+
+/* converge: from two points of differing sign, converge to zero crossing */
+/* watch it: p1 and p2 are used to calculate */
+static void converge(PROCESS *process, const float p1[3], const float p2[3], float v1, float v2,
+                     float p[3], MetaBall *mb, int f)
+{
+       int i = 0;
+       float pos[3], neg[3];
+       float positive = 0.0f, negative = 0.0f;
+       float dvec[3];
+       
+       if (v1 < 0) {
+               copy_v3_v3(pos, p2);
+               copy_v3_v3(neg, p1);
+               positive = v2;
+               negative = v1;
+       }
+       else {
+               copy_v3_v3(pos, p1);
+               copy_v3_v3(neg, p2);
+               positive = v1;
+               negative = v2;
+       }
+
+       sub_v3_v3v3(dvec, pos, neg);
+
+/* Approximation by linear interpolation is faster then binary subdivision,
+ * but it results sometimes (mb->thresh < 0.2) into the strange results */
+       if ((mb->thresh > 0.2f) && (f == 1)) {
+               if ((dvec[1] == 0.0f) && (dvec[2] == 0.0f)) {
+                       p[0] = neg[0] - negative * dvec[0] / (positive - negative);
+                       p[1] = neg[1];
+                       p[2] = neg[2];
+                       return;
+               }
+               if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
+                       p[0] = neg[0];
+                       p[1] = neg[1] - negative * dvec[1] / (positive - negative);
+                       p[2] = neg[2];
+                       return;
+               }
+               if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
+                       p[0] = neg[0];
+                       p[1] = neg[1];
+                       p[2] = neg[2] - negative * dvec[2] / (positive - negative);
+                       return;
+               }
+       }
+
+       if ((dvec[1] == 0.0f) && (dvec[2] == 0.0f)) {
+               p[1] = neg[1];
+               p[2] = neg[2];
+               while (1) {
+                       if (i++ == RES) return;
+                       p[0] = 0.5f * (pos[0] + neg[0]);
+                       if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[0] = p[0];
+                       else                                                       neg[0] = p[0];
+               }
+       }
+
+       if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
+               p[0] = neg[0];
+               p[2] = neg[2];
+               while (1) {
+                       if (i++ == RES) return;
+                       p[1] = 0.5f * (pos[1] + neg[1]);
+                       if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[1] = p[1];
+                       else                                                       neg[1] = p[1];
+               }
+       }
+
+       if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
+               p[0] = neg[0];
+               p[1] = neg[1];
+               while (1) {
+                       if (i++ == RES) return;
+                       p[2] = 0.5f * (pos[2] + neg[2]);
+                       if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[2] = p[2];
+                       else                                                       neg[2] = p[2];
+               }
+       }
+
+       /* This is necessary to find start point */
+       while (1) {
+               mid_v3_v3v3(&p[0], pos, neg);
+
+               if (i++ == RES) {
+                       return;
+               }
+
+               if ((process->function(process, p[0], p[1], p[2])) > 0.0f) {
+                       copy_v3_v3(pos, &p[0]);
+               }
+               else {
+                       copy_v3_v3(neg, &p[0]);
+               }
+       }
+}
+
+/* ************************************** */
+static void add_cube(PROCESS *process, int i, int j, int k, int count)
+{
+       CUBES *ncube;
+       int n;
+       int a, b, c;
+
+       /* hmmm, not only one, but eight cube will be added on the stack 
+        * ... */
+       for (a = i - 1; a < i + count; a++)
+               for (b = j - 1; b < j + count; b++)
+                       for (c = k - 1; c < k + count; c++) {
+                               /* test if cube has been found before */
+                               if (setcenter(process, process->centers, a, b, c) == 0) {
+                                       /* push cube on stack: */
+                                       ncube = (CUBES *) new_pgn_element(process, sizeof(CUBES));
+                                       ncube->next = process->cubes;
+                                       process->cubes = ncube;
+
+                                       ncube->cube.i = a;
+                                       ncube->cube.j = b;
+                                       ncube->cube.k = c;
+
+                                       /* set corners of initial cube: */
+                                       for (n = 0; n < 8; n++)
+                                               ncube->cube.corners[n] = setcorner(process, a + MB_BIT(n, 2), b + MB_BIT(n, 1), c + MB_BIT(n, 0));
+                               }
+                       }
+}
+
+
+static void find_first_points(PROCESS *process, MetaBall *mb, int a)
+{
+       MetaElem *ml;
+       float f;
+
+       ml = process->mainb[a];
+       f = 1.0f - (mb->thresh / ml->s);
+
+       /* Skip, when Stiffness of MetaElement is too small ... MetaElement can't be
+        * visible alone ... but still can influence others MetaElements :-) */
+       if (f > 0.0f) {
+               float IN[3] = {0.0f}, OUT[3] = {0.0f}, in[3] = {0.0f}, out[3];
+               int i, j, k, c_i, c_j, c_k;
+               int index[3] = {1, 0, -1};
+               float in_v /*, out_v*/;
+               float workp[3];
+               float dvec[3];
+               float tmp_v, workp_v, max_len_sq, nx, ny, nz, max_dim;
+
+               calc_mballco(ml, in);
+               in_v = process->function(process, in[0], in[1], in[2]);
+
+               for (i = 0; i < 3; i++) {
+                       switch (ml->type) {
+                               case MB_BALL:
+                                       OUT[0] = out[0] = IN[0] + index[i] * ml->rad;
+                                       break;
+                               case MB_TUBE:
+                               case MB_PLANE:
+                               case MB_ELIPSOID:
+                               case MB_CUBE:
+                                       OUT[0] = out[0] = IN[0] + index[i] * (ml->expx + ml->rad);
+                                       break;
+                       }
+
+                       for (j = 0; j < 3; j++) {
+                               switch (ml->type) {
+                                       case MB_BALL:
+                                               OUT[1] = out[1] = IN[1] + index[j] * ml->rad;
+                                               break;
+                                       case MB_TUBE:
+                                       case MB_PLANE:
+                                       case MB_ELIPSOID:
+                                       case MB_CUBE:
+                                               OUT[1] = out[1] = IN[1] + index[j] * (ml->expy + ml->rad);
+                                               break;
+                               }
+                       
+                               for (k = 0; k < 3; k++) {
+                                       out[0] = OUT[0];
+                                       out[1] = OUT[1];
+                                       switch (ml->type) {
+                                               case MB_BALL:
+                                               case MB_TUBE:
+                                               case MB_PLANE:
+                                                       out[2] = IN[2] + index[k] * ml->rad;
+                                                       break;
+                                               case MB_ELIPSOID:
+                                               case MB_CUBE:
+                                                       out[2] = IN[2] + index[k] * (ml->expz + ml->rad);
+                                                       break;
+                                       }
+
+                                       calc_mballco(ml, out);
+
+                                       /*out_v = process->function(out[0], out[1], out[2]);*/ /*UNUSED*/
+
+                                       /* find "first points" on Implicit Surface of MetaElemnt ml */
+                                       copy_v3_v3(workp, in);
+                                       workp_v = in_v;
+                                       max_len_sq = len_squared_v3v3(out, in);
+
+                                       nx = fabsf((out[0] - in[0]) / process->size);
+                                       ny = fabsf((out[1] - in[1]) / process->size);
+                                       nz = fabsf((out[2] - in[2]) / process->size);
+                                       
+                                       max_dim = max_fff(nx, ny, nz);
+                                       if (max_dim != 0.0f) {
+                                               float len_sq = 0.0f;
+
+                                               dvec[0] = (out[0] - in[0]) / max_dim;
+                                               dvec[1] = (out[1] - in[1]) / max_dim;
+                                               dvec[2] = (out[2] - in[2]) / max_dim;
+
+                                               while (len_sq <= max_len_sq) {
+                                                       add_v3_v3(workp, dvec);
+
+                                                       /* compute value of implicite function */
+                                                       tmp_v = process->function(process, workp[0], workp[1], workp[2]);
+                                                       /* add cube to the stack, when value of implicite function crosses zero value */
+                                                       if ((tmp_v < 0.0f && workp_v >= 0.0f) || (tmp_v > 0.0f && workp_v <= 0.0f)) {
+
+                                                               /* indexes of CUBE, which includes "first point" */
+                                                               c_i = (int)floor(workp[0] / process->size);
+                                                               c_j = (int)floor(workp[1] / process->size);
+                                                               c_k = (int)floor(workp[2] / process->size);
+                                                               
+                                                               /* add CUBE (with indexes c_i, c_j, c_k) to the stack,
+                                                                * this cube includes found point of Implicit Surface */
+                                                               if ((ml->flag & MB_NEGATIVE) == 0) {
+                                                                       add_cube(process, c_i, c_j, c_k, 1);
+                                                               }
+                                                               else {
+                                                                       add_cube(process, c_i, c_j, c_k, 2);
+                                                               }
+                                                       }
+                                                       len_sq = len_squared_v3v3(workp, in);
+                                                       workp_v = tmp_v;
+
+                                               }
+                                       }
+                               }
+                       }
+               }
+       }
+}
+
+static void polygonize(PROCESS *process, MetaBall *mb)
+{
+       CUBE c;
+       int a;
+
+       process->vertices.count = process->vertices.max = 0;
+       process->vertices.ptr = NULL;
+
+       /* allocate hash tables and build cube polygon table: */
+       process->centers = MEM_callocN(HASHSIZE * sizeof(CENTERLIST *), "mbproc->centers");
+       process->corners = MEM_callocN(HASHSIZE * sizeof(CORNER *), "mbproc->corners");
+       process->edges = MEM_callocN(2 * HASHSIZE * sizeof(EDGELIST *), "mbproc->edges");
+       makecubetable();
+
+       for (a = 0; a < process->totelem; a++) {
+
+               /* try to find 8 points on the surface for each MetaElem */
+               find_first_points(process, mb, a);
+       }
+
+       /* polygonize all MetaElems of current MetaBall */
+       while (process->cubes != NULL) { /* process active cubes till none left */
+               c = process->cubes->cube;
+
+               /* polygonize the cube directly: */
+               docube(process, &c, mb);
+               
+               /* pop current cube from stack */
+               process->cubes = process->cubes->next;
+               
+               /* test six face directions, maybe add to stack: */
+               testface(process, c.i - 1, c.j, c.k, &c, 2, LBN, LBF, LTN, LTF);
+               testface(process, c.i + 1, c.j, c.k, &c, 2, RBN, RBF, RTN, RTF);
+               testface(process, c.i, c.j - 1, c.k, &c, 1, LBN, LBF, RBN, RBF);
+               testface(process, c.i, c.j + 1, c.k, &c, 1, LTN, LTF, RTN, RTF);
+               testface(process, c.i, c.j, c.k - 1, &c, 0, LBN, LTN, RBN, RTN);
+               testface(process, c.i, c.j, c.k + 1, &c, 0, LBF, LTF, RBF, RTF);
+       }
+}
+
+static float init_meta(EvaluationContext *eval_ctx, PROCESS *process, Scene *scene, Object *ob)    /* return totsize */
+{
+       Scene *sce_iter = scene;
+       Base *base;
+       Object *bob;
+       MetaBall *mb;
+       MetaElem *ml;
+       float size, totsize, obinv[4][4], obmat[4][4], vec[3];
+       //float max = 0.0f;
+       int a, obnr, zero_size = 0;
+       char obname[MAX_ID_NAME];
+       SceneBaseIter iter;
+
+       copy_m4_m4(obmat, ob->obmat);   /* to cope with duplicators from BKE_scene_base_iter_next */
+       invert_m4_m4(obinv, ob->obmat);
+       a = 0;
+       
+       BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
+       
+       /* make main array */
+       BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 0, NULL, NULL);
+       while (BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 1, &base, &bob)) {
+
+               if (bob->type == OB_MBALL) {
+                       zero_size = 0;
+                       ml = NULL;
+
+                       if (bob == ob && (base->flag & OB_FROMDUPLI) == 0) {
+                               mb = ob->data;
+       
+                               if (mb->editelems) ml = mb->editelems->first;
+                               else ml = mb->elems.first;
+                       }
+                       else {
+                               char name[MAX_ID_NAME];
+                               int nr;
+                               
+                               BLI_split_name_num(name, &nr, bob->id.name + 2, '.');
+                               if (STREQ(obname, name)) {
+                                       mb = bob->data;
+                                       
+                                       if (mb->editelems) ml = mb->editelems->first;
+                                       else ml = mb->elems.first;
+                               }
+                       }
+
+                       /* when metaball object has zero scale, then MetaElem to this MetaBall
+                        * will not be put to mainb array */
+                       if (has_zero_axis_m4(bob->obmat)) {
+                               zero_size = 1;
+                       }
+                       else if (bob->parent) {
+                               struct Object *pob = bob->parent;
+                               while (pob) {
+                                       if (has_zero_axis_m4(pob->obmat)) {
+                                               zero_size = 1;
+                                               break;
+                                       }
+                                       pob = pob->parent;
+                               }
+                       }
+
+                       if (zero_size) {
+                               unsigned int ml_count = 0;
+                               while (ml) {
+                                       ml_count++;
+                                       ml = ml->next;
+                               }
+                               process->totelem -= ml_count;
+                       }
+                       else {
+                               while (ml) {
+                                       if (!(ml->flag & MB_HIDE)) {
+                                               int i;
+                                               float temp1[4][4], temp2[4][4], temp3[4][4];
+                                               float (*mat)[4] = NULL, (*imat)[4] = NULL;
+                                               float max_x, max_y, max_z, min_x, min_y, min_z;
+                                               float expx, expy, expz;
+
+                                               max_x = max_y = max_z = -3.4e38;
+                                               min_x = min_y = min_z =  3.4e38;
+
+                                               /* too big stiffness seems only ugly due to linear interpolation
+                                                * no need to have possibility for too big stiffness */
+                                               if (ml->s > 10.0f) ml->s = 10.0f;
+
+                                               /* Rotation of MetaElem is stored in quat */
+                                               quat_to_mat4(temp3, ml->quat);
+
+                                               /* Translation of MetaElem */
+                                               unit_m4(temp2);
+                                               temp2[3][0] = ml->x;
+                                               temp2[3][1] = ml->y;
+                                               temp2[3][2] = ml->z;
+
+                                               mul_m4_m4m4(temp1, temp2, temp3);
+
+                                               /* make a copy because of duplicates */
+                                               process->mainb[a] = new_pgn_element(process, sizeof(MetaElem));
+                                               *(process->mainb[a]) = *ml;
+                                               process->mainb[a]->bb = new_pgn_element(process, sizeof(BoundBox));
+
+                                               mat = new_pgn_element(process, 4 * 4 * sizeof(float));
+                                               imat = new_pgn_element(process, 4 * 4 * sizeof(float));
+
+                                               /* mat is the matrix to transform from mball into the basis-mball */
+                                               invert_m4_m4(obinv, obmat);
+                                               mul_m4_m4m4(temp2, obinv, bob->obmat);
+                                               /* MetaBall transformation */
+                                               mul_m4_m4m4(mat, temp2, temp1);
+
+                                               invert_m4_m4(imat, mat);
+
+                                               process->mainb[a]->rad2 = ml->rad * ml->rad;
+
+                                               process->mainb[a]->mat = (float *) mat;
+                                               process->mainb[a]->imat = (float *) imat;
+
+                                               if (!MB_TYPE_SIZE_SQUARED(ml->type)) {
+                                                       expx = ml->expx;
+                                                       expy = ml->expy;
+                                                       expz = ml->expz;
+                                               }
+                                               else {
+                                                       expx = ml->expx * ml->expx;
+                                                       expy = ml->expy * ml->expy;
+                                                       expz = ml->expz * ml->expz;
+                                               }
+
+                                               /* untransformed Bounding Box of MetaElem */
+                                               /* TODO, its possible the elem type has been changed and the exp* values can use a fallback */
+                                               copy_v3_fl3(process->mainb[a]->bb->vec[0], -expx, -expy, -expz);  /* 0 */
+                                               copy_v3_fl3(process->mainb[a]->bb->vec[1], +expx, -expy, -expz);  /* 1 */
+                                               copy_v3_fl3(process->mainb[a]->bb->vec[2], +expx, +expy, -expz);  /* 2 */
+                                               copy_v3_fl3(process->mainb[a]->bb->vec[3], -expx, +expy, -expz);  /* 3 */
+                                               copy_v3_fl3(process->mainb[a]->bb->vec[4], -expx, -expy, +expz);  /* 4 */
+                                               copy_v3_fl3(process->mainb[a]->bb->vec[5], +expx, -expy, +expz);  /* 5 */
+                                               copy_v3_fl3(process->mainb[a]->bb->vec[6], +expx, +expy, +expz);  /* 6 */
+                                               copy_v3_fl3(process->mainb[a]->bb->vec[7], -expx, +expy, +expz);  /* 7 */
+
+                                               /* transformation of Metalem bb */
+                                               for (i = 0; i < 8; i++)
+                                                       mul_m4_v3((float (*)[4])mat, process->mainb[a]->bb->vec[i]);
+
+                                               /* find max and min of transformed bb */
+                                               for (i = 0; i < 8; i++) {
+                                                       /* find maximums */
+                                                       if (process->mainb[a]->bb->vec[i][0] > max_x) max_x = process->mainb[a]->bb->vec[i][0];
+                                                       if (process->mainb[a]->bb->vec[i][1] > max_y) max_y = process->mainb[a]->bb->vec[i][1];
+                                                       if (process->mainb[a]->bb->vec[i][2] > max_z) max_z = process->mainb[a]->bb->vec[i][2];
+                                                       /* find  minimums */
+                                                       if (process->mainb[a]->bb->vec[i][0] < min_x) min_x = process->mainb[a]->bb->vec[i][0];
+                                                       if (process->mainb[a]->bb->vec[i][1] < min_y) min_y = process->mainb[a]->bb->vec[i][1];
+                                                       if (process->mainb[a]->bb->vec[i][2] < min_z) min_z = process->mainb[a]->bb->vec[i][2];
+                                               }
+                                       
+                                               /* create "new" bb, only point 0 and 6, which are
+                                                * necessary for octal tree filling */
+                                               process->mainb[a]->bb->vec[0][0] = min_x - ml->rad;
+                                               process->mainb[a]->bb->vec[0][1] = min_y - ml->rad;
+                                               process->mainb[a]->bb->vec[0][2] = min_z - ml->rad;
+
+                                               process->mainb[a]->bb->vec[6][0] = max_x + ml->rad;
+                                               process->mainb[a]->bb->vec[6][1] = max_y + ml->rad;
+                                               process->mainb[a]->bb->vec[6][2] = max_z + ml->rad;
+
+                                               a++;
+                                       }
+                                       ml = ml->next;
+                               }
+                       }
+               }
+       }
+
+       
+       /* totsize (= 'manhattan' radius) */
+       totsize = 0.0;
+       for (a = 0; a < process->totelem; a++) {
+               
+               vec[0] = process->mainb[a]->x + process->mainb[a]->rad + process->mainb[a]->expx;
+               vec[1] = process->mainb[a]->y + process->mainb[a]->rad + process->mainb[a]->expy;
+               vec[2] = process->mainb[a]->z + process->mainb[a]->rad + process->mainb[a]->expz;
+
+               calc_mballco(process->mainb[a], vec);
+       
+               size = fabsf(vec[0]);
+               if (size > totsize) totsize = size;
+               size = fabsf(vec[1]);
+               if (size > totsize) totsize = size;
+               size = fabsf(vec[2]);
+               if (size > totsize) totsize = size;
+
+               vec[0] = process->mainb[a]->x - process->mainb[a]->rad;
+               vec[1] = process->mainb[a]->y - process->mainb[a]->rad;
+               vec[2] = process->mainb[a]->z - process->mainb[a]->rad;
+                               
+               calc_mballco(process->mainb[a], vec);
+       
+               size = fabsf(vec[0]);
+               if (size > totsize) totsize = size;
+               size = fabsf(vec[1]);
+               if (size > totsize) totsize = size;
+               size = fabsf(vec[2]);
+               if (size > totsize) totsize = size;
+       }
+
+       for (a = 0; a < process->totelem; a++) {
+               process->thresh += densfunc(process->mainb[a], 2.0f * totsize, 2.0f * totsize, 2.0f * totsize);
+       }
+
+       return totsize;
+}
+
+/* if MetaElem lies in node, then node includes MetaElem pointer (ml_p)
+ * pointing at MetaElem (ml)
+ */
+static void fill_metaball_octal_node(octal_node *node, MetaElem *ml, short i)
+{
+       ml_pointer *ml_p;
+
+       ml_p = MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
+       ml_p->ml = ml;
+       BLI_addtail(&(node->nodes[i]->elems), ml_p);
+       node->count++;
+       
+       if ((ml->flag & MB_NEGATIVE) == 0) {
+               node->nodes[i]->pos++;
+       }
+       else {
+               node->nodes[i]->neg++;
+       }
+}
+
+/* Node is subdivided as is illustrated on the following figure:
+ * 
+ *      +------+------+
+ *     /      /      /|
+ *    +------+------+ |
+ *   /      /      /| +
+ *  +------+------+ |/|
+ *  |      |      | + |
+ *  |      |      |/| +
+ *  +------+------+ |/
+ *  |      |      | +
+ *  |      |      |/
+ *  +------+------+
+ *  
+ */
+static void subdivide_metaball_octal_node(octal_node *node, float size_x, float size_y, float size_z, short depth)
+{
+       MetaElem *ml;
+       ml_pointer *ml_p;
+       float x, y, z;
+       int a, i;
+
+       /* create new nodes */
+       for (a = 0; a < 8; a++) {
+               node->nodes[a] = MEM_mallocN(sizeof(octal_node), "octal_node");
+               for (i = 0; i < 8; i++)
+                       node->nodes[a]->nodes[i] = NULL;
+               node->nodes[a]->parent = node;
+               BLI_listbase_clear(&node->nodes[a]->elems);
+               node->nodes[a]->count = 0;
+               node->nodes[a]->neg = 0;
+               node->nodes[a]->pos = 0;
+       }
+
+       size_x /= 2;
+       size_y /= 2;
+       size_z /= 2;
+       
+       /* center of node */
+       node->x = x = node->x_min + size_x;
+       node->y = y = node->y_min + size_y;
+       node->z = z = node->z_min + size_z;
+
+       /* setting up of border points of new nodes */
+       node->nodes[0]->x_min = node->x_min;
+       node->nodes[0]->y_min = node->y_min;
+       node->nodes[0]->z_min = node->z_min;
+       node->nodes[0]->x = node->nodes[0]->x_min + size_x / 2;
+       node->nodes[0]->y = node->nodes[0]->y_min + size_y / 2;
+       node->nodes[0]->z = node->nodes[0]->z_min + size_z / 2;
+       
+       node->nodes[1]->x_min = x;
+       node->nodes[1]->y_min = node->y_min;
+       node->nodes[1]->z_min = node->z_min;
+       node->nodes[1]->x = node->nodes[1]->x_min + size_x / 2;
+       node->nodes[1]->y = node->nodes[1]->y_min + size_y / 2;
+       node->nodes[1]->z = node->nodes[1]->z_min + size_z / 2;
+
+       node->nodes[2]->x_min = x;
+       node->nodes[2]->y_min = y;
+       node->nodes[2]->z_min = node->z_min;
+       node->nodes[2]->x = node->nodes[2]->x_min + size_x / 2;
+       node->nodes[2]->y = node->nodes[2]->y_min + size_y / 2;
+       node->nodes[2]->z = node->nodes[2]->z_min + size_z / 2;
+
+       node->nodes[3]->x_min = node->x_min;
+       node->nodes[3]->y_min = y;
+       node->nodes[3]->z_min = node->z_min;
+       node->nodes[3]->x = node->nodes[3]->x_min + size_x / 2;
+       node->nodes[3]->y = node->nodes[3]->y_min + size_y / 2;
+       node->nodes[3]->z = node->nodes[3]->z_min + size_z / 2;
+
+       node->nodes[4]->x_min = node->x_min;
+       node->nodes[4]->y_min = node->y_min;
+       node->nodes[4]->z_min = z;
+       node->nodes[4]->x = node->nodes[4]->x_min + size_x / 2;
+       node->nodes[4]->y = node->nodes[4]->y_min + size_y / 2;
+       node->nodes[4]->z = node->nodes[4]->z_min + size_z / 2;
+       
+       node->nodes[5]->x_min = x;
+       node->nodes[5]->y_min = node->y_min;
+       node->nodes[5]->z_min = z;
+       node->nodes[5]->x = node->nodes[5]->x_min + size_x / 2;
+       node->nodes[5]->y = node->nodes[5]->y_min + size_y / 2;
+       node->nodes[5]->z = node->nodes[5]->z_min + size_z / 2;
+
+       node->nodes[6]->x_min = x;
+       node->nodes[6]->y_min = y;
+       node->nodes[6]->z_min = z;
+       node->nodes[6]->x = node->nodes[6]->x_min + size_x / 2;
+       node->nodes[6]->y = node->nodes[6]->y_min + size_y / 2;
+       node->nodes[6]->z = node->nodes[6]->z_min + size_z / 2;
+
+       node->nodes[7]->x_min = node->x_min;
+       node->nodes[7]->y_min = y;
+       node->nodes[7]->z_min = z;
+       node->nodes[7]->x = node->nodes[7]->x_min + size_x / 2;
+       node->nodes[7]->y = node->nodes[7]->y_min + size_y / 2;
+       node->nodes[7]->z = node->nodes[7]->z_min + size_z / 2;
+
+       ml_p = node->elems.first;
+       
+       /* setting up references of MetaElems for new nodes */
+       while (ml_p) {
+               ml = ml_p->ml;
+               if (ml->bb->vec[0][2] < z) {
+                       if (ml->bb->vec[0][1] < y) {
+                               /* vec[0][0] lies in first octant */
+                               if (ml->bb->vec[0][0] < x) {
+                                       /* ml belongs to the (0)1st node */
+                                       fill_metaball_octal_node(node, ml, 0);
+
+                                       /* ml belongs to the (3)4th node */
+                                       if (ml->bb->vec[6][1] >= y) {
+                                               fill_metaball_octal_node(node, ml, 3);
+
+                                               /* ml belongs to the (7)8th node */
+                                               if (ml->bb->vec[6][2] >= z) {
+                                                       fill_metaball_octal_node(node, ml, 7);
+                                               }
+                                       }
+       
+                                       /* ml belongs to the (1)2nd node */
+                                       if (ml->bb->vec[6][0] >= x) {
+                                               fill_metaball_octal_node(node, ml, 1);
+
+                                               /* ml belongs to the (5)6th node */
+                                               if (ml->bb->vec[6][2] >= z) {
+                                                       fill_metaball_octal_node(node, ml, 5);
+                                               }
+                                       }
+
+                                       /* ml belongs to the (2)3th node */
+                                       if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
+                                               fill_metaball_octal_node(node, ml, 2);
+                                               
+                                               /* ml belong to the (6)7th node */
+                                               if (ml->bb->vec[6][2] >= z) {
+                                                       fill_metaball_octal_node(node, ml, 6);
+                                               }
+                                               
+                                       }
+                       
+                                       /* ml belongs to the (4)5th node too */
+                                       if (ml->bb->vec[6][2] >= z) {
+                                               fill_metaball_octal_node(node, ml, 4);
+                                       }
+
+                                       
+                                       
+                               }
+                               /* vec[0][0] is in the (1)second octant */
+                               else {
+                                       /* ml belong to the (1)2nd node */
+                                       fill_metaball_octal_node(node, ml, 1);
+
+                                       /* ml belongs to the (2)3th node */
+                                       if (ml->bb->vec[6][1] >= y) {
+                                               fill_metaball_octal_node(node, ml, 2);
+
+                                               /* ml belongs to the (6)7th node */
+                                               if (ml->bb->vec[6][2] >= z) {
+                                                       fill_metaball_octal_node(node, ml, 6);
+                                               }
+                                               
+                                       }
+                                       
+                                       /* ml belongs to the (5)6th node */
+                                       if (ml->bb->vec[6][2] >= z) {
+                                               fill_metaball_octal_node(node, ml, 5);
+                                       }
+                               }
+                       }
+                       else {
+                               /* vec[0][0] is in the (3)4th octant */
+                               if (ml->bb->vec[0][0] < x) {
+                                       /* ml belongs to the (3)4nd node */
+                                       fill_metaball_octal_node(node, ml, 3);
+                                       
+                                       /* ml belongs to the (7)8th node */
+                                       if (ml->bb->vec[6][2] >= z) {
+                                               fill_metaball_octal_node(node, ml, 7);
+                                       }
+                               
+
+                                       /* ml belongs to the (2)3th node */
+                                       if (ml->bb->vec[6][0] >= x) {
+                                               fill_metaball_octal_node(node, ml, 2);
+                                       
+                                               /* ml belongs to the (6)7th node */
+                                               if (ml->bb->vec[6][2] >= z) {
+                                                       fill_metaball_octal_node(node, ml, 6);
+                                               }
+                                       }
+                               }
+
+                       }
+
+                       /* vec[0][0] is in the (2)3th octant */
+                       if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
+                               /* ml belongs to the (2)3th node */
+                               fill_metaball_octal_node(node, ml, 2);
+                               
+                               /* ml belongs to the (6)7th node */
+                               if (ml->bb->vec[6][2] >= z) {
+                                       fill_metaball_octal_node(node, ml, 6);
+                               }
+                       }
+               }
+               else {
+                       if (ml->bb->vec[0][1] < y) {
+                               /* vec[0][0] lies in (4)5th octant */
+                               if (ml->bb->vec[0][0] < x) {
+                                       /* ml belongs to the (4)5th node */
+                                       fill_metaball_octal_node(node, ml, 4);
+
+                                       if (ml->bb->vec[6][0] >= x) {
+                                               fill_metaball_octal_node(node, ml, 5);
+                                       }
+
+                                       if (ml->bb->vec[6][1] >= y) {
+                                               fill_metaball_octal_node(node, ml, 7);
+                                       }
+                                       
+                                       if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
+                                               fill_metaball_octal_node(node, ml, 6);
+                                       }
+                               }
+                               /* vec[0][0] lies in (5)6th octant */
+                               else {
+                                       fill_metaball_octal_node(node, ml, 5);
+
+                                       if (ml->bb->vec[6][1] >= y) {
+                                               fill_metaball_octal_node(node, ml, 6);
+                                       }
+                               }
+                       }
+                       else {
+                               /* vec[0][0] lies in (7)8th octant */
+                               if (ml->bb->vec[0][0] < x) {
+                                       fill_metaball_octal_node(node, ml, 7);
+
+                                       if (ml->bb->vec[6][0] >= x) {
+                                               fill_metaball_octal_node(node, ml, 6);
+                                       }
+                               }
+
+                       }
+                       
+                       /* vec[0][0] lies in (6)7th octant */
+                       if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
+                               fill_metaball_octal_node(node, ml, 6);
+                       }
+               }
+               ml_p = ml_p->next;
+       }
+
+       /* free references of MetaElems for curent node (it is not needed anymore) */
+       BLI_freelistN(&node->elems);
+
+       depth--;
+       
+       if (depth > 0) {
+               for (a = 0; a < 8; a++) {
+                       if (node->nodes[a]->count > 0) /* if node is not empty, then it is subdivided */
+                               subdivide_metaball_octal_node(node->nodes[a], size_x, size_y, size_z, depth);
+               }
+       }
+}
+
+/* free all octal nodes recursively */
+static void free_metaball_octal_node(octal_node *node)
+{
+       int a;
+       for (a = 0; a < 8; a++) {
+               if (node->nodes[a] != NULL) free_metaball_octal_node(node->nodes[a]);
+       }
+       BLI_freelistN(&node->elems);
+       MEM_freeN(node);
+}
+
+/* If scene include more than one MetaElem, then octree is used */
+static void init_metaball_octal_tree(PROCESS *process, int depth)
+{
+       struct octal_node *node;
+       ml_pointer *ml_p;
+       float size[3];
+       int a;
+       
+       process->metaball_tree = MEM_mallocN(sizeof(octal_tree), "metaball_octal_tree");
+       process->metaball_tree->first = node = MEM_mallocN(sizeof(octal_node), "metaball_octal_node");
+       /* maximal depth of octree */
+       process->metaball_tree->depth = depth;
+
+       process->metaball_tree->neg = node->neg = 0;
+       process->metaball_tree->pos = node->pos = 0;
+       
+       BLI_listbase_clear(&node->elems);
+       node->count = 0;
+
+       for (a = 0; a < 8; a++)
+               node->nodes[a] = NULL;
+
+       node->x_min = node->y_min = node->z_min = FLT_MAX;
+       node->x_max = node->y_max = node->z_max = -FLT_MAX;
+
+       /* size of octal tree scene */
+       for (a = 0; a < process->totelem; a++) {
+               if (process->mainb[a]->bb->vec[0][0] < node->x_min) node->x_min = process->mainb[a]->bb->vec[0][0];
+               if (process->mainb[a]->bb->vec[0][1] < node->y_min) node->y_min = process->mainb[a]->bb->vec[0][1];
+               if (process->mainb[a]->bb->vec[0][2] < node->z_min) node->z_min = process->mainb[a]->bb->vec[0][2];
+
+               if (process->mainb[a]->bb->vec[6][0] > node->x_max) node->x_max = process->mainb[a]->bb->vec[6][0];
+               if (process->mainb[a]->bb->vec[6][1] > node->y_max) node->y_max = process->mainb[a]->bb->vec[6][1];
+               if (process->mainb[a]->bb->vec[6][2] > node->z_max) node->z_max = process->mainb[a]->bb->vec[6][2];
+
+               ml_p = MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
+               ml_p->ml = process->mainb[a];
+               BLI_addtail(&node->elems, ml_p);
+
+               if ((process->mainb[a]->flag & MB_NEGATIVE) == 0) {
+                       /* number of positive MetaElem in scene */
+                       process->metaball_tree->pos++;
+               }
+               else {
+                       /* number of negative MetaElem in scene */
+                       process->metaball_tree->neg++;
+               }
+       }
+
+       /* size of first node */
+       size[0] = node->x_max - node->x_min;
+       size[1] = node->y_max - node->y_min;
+       size[2] = node->z_max - node->z_min;
+
+       /* first node is subdivided recursively */
+       subdivide_metaball_octal_node(node, size[0], size[1], size[2], process->metaball_tree->depth);
+}
+
+static void mball_count(EvaluationContext *eval_ctx, PROCESS *process, Scene *scene, Object *basis)
+{
+       Scene *sce_iter = scene;
+       Base *base;
+       Object *ob, *bob = basis;
+       MetaElem *ml = NULL;
+       int basisnr, obnr;
+       char basisname[MAX_ID_NAME], obname[MAX_ID_NAME];
+       SceneBaseIter iter;
+
+       BLI_split_name_num(basisname, &basisnr, basis->id.name + 2, '.');
+       process->totelem = 0;
+
+       BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 0, NULL, NULL);
+       while (BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 1, &base, &ob)) {
+               if (ob->type == OB_MBALL) {
+                       if (ob == bob) {
+                               MetaBall *mb = ob->data;
+
+                               /* if bob object is in edit mode, then dynamic list of all MetaElems
+                                * is stored in editelems */
+                               if (mb->editelems) ml = mb->editelems->first;
+                               /* if bob object is in object mode */
+                               else ml = mb->elems.first;
+                       }
+                       else {
+                               BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
+
+                               /* object ob has to be in same "group" ... it means, that it has to have
+                                * same base of its name */
+                               if (STREQ(obname, basisname)) {
+                                       MetaBall *mb = ob->data;
+
+                                       /* if object is in edit mode, then dynamic list of all MetaElems
+                                        * is stored in editelems */
+                                       if (mb->editelems) ml = mb->editelems->first;
+                                       /* if bob object is in object mode */
+                                       else ml = mb->elems.first;
+                               }
+                       }
+
+                       for ( ; ml; ml = ml->next) {
+                               if (!(ml->flag & MB_HIDE)) {
+                                       process->totelem++;
+                               }
+                       }
+               }
+       }
+}
+
+void BKE_mball_polygonize(EvaluationContext *eval_ctx, Scene *scene, Object *ob, ListBase *dispbase)
+{
+       MetaBall *mb;
+       DispList *dl;
+       int a, nr_cubes;
+       float *co, *no, totsize, width;
+       PROCESS process = {0};
+
+       mb = ob->data;
+
+       mball_count(eval_ctx, &process, scene, ob);
+
+       if (process.totelem == 0) return;
+       if ((eval_ctx->mode != DAG_EVAL_RENDER) && (mb->flag == MB_UPDATE_NEVER)) return;
+       if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_FAST) return;
+
+       process.thresh = mb->thresh;
+
+       /* total number of MetaElems (totelem) is precomputed in find_basis_mball() function */
+       process.mainb = MEM_mallocN(sizeof(void *) * process.totelem, "mainb");
+       
+       /* initialize all mainb (MetaElems) */
+       totsize = init_meta(eval_ctx, &process, scene, ob);
+
+       /* if scene includes more than one MetaElem, then octal tree optimization is used */
+       if ((process.totelem >    1) && (process.totelem <=   64)) init_metaball_octal_tree(&process, 1);
+       if ((process.totelem >   64) && (process.totelem <=  128)) init_metaball_octal_tree(&process, 2);
+       if ((process.totelem >  128) && (process.totelem <=  512)) init_metaball_octal_tree(&process, 3);
+       if ((process.totelem >  512) && (process.totelem <= 1024)) init_metaball_octal_tree(&process, 4);
+       if (process.totelem  > 1024)                               init_metaball_octal_tree(&process, 5);
+
+       /* don't polygonize metaballs with too high resolution (base mball to small)
+        * note: Eps was 0.0001f but this was giving problems for blood animation for durian, using 0.00001f */
+       if (process.metaball_tree) {
+               if (ob->size[0] <= 0.00001f * (process.metaball_tree->first->x_max - process.metaball_tree->first->x_min) ||
+                   ob->size[1] <= 0.00001f * (process.metaball_tree->first->y_max - process.metaball_tree->first->y_min) ||
+                   ob->size[2] <= 0.00001f * (process.metaball_tree->first->z_max - process.metaball_tree->first->z_min))
+               {
+                       new_pgn_element(&process, -1); /* free values created by init_meta */
+
+                       MEM_freeN(process.mainb);
+
+                       /* free tree */
+                       free_metaball_octal_node(process.metaball_tree->first);
+                       MEM_freeN(process.metaball_tree);
+
+                       return;
+               }
+       }
+
+       /* width is size per polygonize cube */
+       if (eval_ctx->mode == DAG_EVAL_RENDER) {
+               width = mb->rendersize;
+       }
+       else {
+               width = mb->wiresize;
+               if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_HALFRES) {
+                       width *= 2;
+               }
+       }
+       /* nr_cubes is just for safety, minimum is totsize */
+       nr_cubes = (int)(0.5f + totsize / width);
+
+       /* init process */
+       process.function = metaball;
+       process.size = width;
+       process.bounds = nr_cubes;
+       process.cubes = NULL;
+       process.delta = width / (float)(RES * RES);
+
+       polygonize(&process, mb);
+       
+       MEM_freeN(process.mainb);
+
+       /* free octal tree */
+       if (process.totelem > 1) {
+               free_metaball_octal_node(process.metaball_tree->first);
+               MEM_freeN(process.metaball_tree);
+               process.metaball_tree = NULL;
+       }
+
+       if (process.curindex) {
+               VERTEX *ptr = process.vertices.ptr;
+
+               dl = MEM_callocN(sizeof(DispList), "mbaldisp");
+               BLI_addtail(dispbase, dl);
+               dl->type = DL_INDEX4;
+               dl->nr = process.vertices.count;
+               dl->parts = process.curindex;
+
+               dl->index = process.indices;
+               process.indices = NULL;
+
+               a = process.vertices.count;
+               dl->verts = co = MEM_mallocN(sizeof(float[3]) * a, "mballverts");
+               dl->nors = no = MEM_mallocN(sizeof(float[3]) * a, "mballnors");
+
+               for (a = 0; a < process.vertices.count; ptr++, a++, no += 3, co += 3) {
+                       copy_v3_v3(co, ptr->co);
+                       copy_v3_v3(no, ptr->no);
+               }
+       }
+
+       freepolygonize(&process);
+}