Metaball tessellation optimization (Octree to BVH)
authorKrzysztof Recko <yetioszek@gmail.com>
Tue, 7 Apr 2015 02:44:39 +0000 (12:44 +1000)
committerCampbell Barton <ideasman42@gmail.com>
Tue, 7 Apr 2015 03:19:50 +0000 (13:19 +1000)
Speedup is non-linear, 2x-10x faster is quite normal.
Patch T43678.

- Switched from an Octree to BVH.
- Finding first points of surface no longer "wastes" density function evaluation: every result is cached.
- Use MemArena instead of using own memory management.
- Correct calculation of metaelem bounding box.
- Remove mball_count(): mballs are now counted "on the go".

source/blender/blenkernel/BKE_mball.h
source/blender/blenkernel/BKE_mball_tessellate.h
source/blender/blenkernel/intern/mball_tessellate.c
source/blender/windowmanager/intern/wm_init_exit.c

index 70f932f3292e4374ca0ebff8d2d4255328e29654..321cbbbd708051199711d39eda5dbe7547b1860c 100644 (file)
@@ -45,8 +45,6 @@ struct MetaBall *BKE_mball_copy(struct MetaBall *mb);
 
 void BKE_mball_make_local(struct MetaBall *mb);
 
-void BKE_mball_cubeTable_free(void);
-
 bool BKE_mball_is_basis_for(struct Object *ob1, struct Object *ob2);
 bool BKE_mball_is_basis(struct Object *ob);
 struct Object *BKE_mball_basis_find(struct Scene *scene, struct Object *ob);
index a3d3e1934b73eae74612916dcd1e6ef4038f691d..361f31b704cdd1fa78339388aa3f5bd88e2ef231 100644 (file)
@@ -31,4 +31,6 @@ void BKE_mball_polygonize(
         struct EvaluationContext *eval_ctx, struct Scene *scene,
         struct Object *ob, struct ListBase *dispbase);
 
+void BKE_mball_cubeTable_free(void);
+
 #endif  /* __BKE_MBALL_TESSELLATE_H__ */
index 04024943310de6648b251c32fce15eb2530f8952..080a8cead7bbaeb5b5e91b23a056d641ffda8794 100644 (file)
 #include "BLI_path_util.h"
 #include "BLI_math.h"
 #include "BLI_utildefines.h"
+#include "BLI_memarena.h"
 
 #include "BKE_global.h"
 
 #include "BKE_depsgraph.h"
 #include "BKE_scene.h"
 #include "BKE_displist.h"
-#include "BKE_mball.h"
 #include "BKE_mball_tessellate.h"  /* own include */
 
-/* Data types */
+#include "BLI_strict_flags.h"
 
-typedef struct vertex {         /* surface vertex */
-       float co[3];  /* position and surface normal */
-       float no[3];
-} VERTEX;
+/* experimental (faster) normal calculation */
+// #define USE_ACCUM_NORMAL
 
-typedef struct vertices {       /* list of vertices in polygonization */
-       int count, max;             /* # vertices, max # allowed */
-       VERTEX *ptr;                /* dynamically allocated */
-} VERTICES;
+/* Data types */
 
 typedef struct corner {         /* corner of a cube */
        int i, j, k;                /* (i, j, k) is index within lattice */
@@ -102,74 +97,161 @@ typedef struct intlists {       /* list of 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 Box {                   /* an AABB with pointer to metalelem */
+       float min[3], max[3];
+       const MetaElem *ml;
+} Box;
+
+typedef struct MetaballBVHNode {       /* BVH node */
+       Box bb[2];                                              /* AABB of children */
+       struct MetaballBVHNode *child[2];
+} MetaballBVHNode;
+
+typedef struct process {        /* parameters, storage */
+       float thresh, size;                     /* mball threshold, single cube size */
+       float delta;                            /* small delta for calculating normals */
+       unsigned int converge_res;      /* converge procedure resolution (more = slower) */
+
+       MetaElem **mainb;                       /* array of all metaelems */
+       unsigned int totelem, mem;      /* number of metaelems */
+
+       MetaballBVHNode metaball_bvh; /* The simplest bvh */
+       Box allbb;                   /* Bounding box of all metaelems */
 
-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 */
+       MetaballBVHNode **bvh_queue; /* Queue used during bvh traversal */
+       unsigned int bvh_queue_size;
+
+       CUBES *cubes;               /* stack of cubes waiting for polygonization */
        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 (*indices)[4];          /* output indices */
+       unsigned int totindex;          /* size of memory allocated for indices */
+       unsigned int curindex;          /* number of currently added indices */
+
+       float (*co)[3], (*no)[3];   /* surface vertices - positions and normals */
+       unsigned int totvertex;         /* memory size */
+       unsigned int curvertex;         /* currently added vertices */
 
-       int pgn_offset;
-       struct pgn_elements *pgn_current;
-       ListBase pgn_list;
+       /* memory allocation from common pool */
+       MemArena *pgn_elements;
 } 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);
+static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2);
+static void add_cube(PROCESS *process, int i, int j, int k);
+static void make_face(PROCESS *process, int i1, int i2, int i3, int i4);
+static void converge(PROCESS *process, const CORNER *c1, const CORNER *c2, float r_p[3]);
+
+/* ******************* SIMPLE BVH ********************* */
+
+static void make_box_union(const BoundBox *a, const Box *b, Box *r_out)
+{
+       r_out->min[0] = min_ff(a->vec[0][0], b->min[0]);
+       r_out->min[1] = min_ff(a->vec[0][1], b->min[1]);
+       r_out->min[2] = min_ff(a->vec[0][2], b->min[2]);
+
+       r_out->max[0] = max_ff(a->vec[6][0], b->max[0]);
+       r_out->max[1] = max_ff(a->vec[6][1], b->max[1]);
+       r_out->max[2] = max_ff(a->vec[6][2], b->max[2]);
+}
+
+static void make_box_from_metaelem(Box *r, const MetaElem *ml)
+{
+       copy_v3_v3(r->max, ml->bb->vec[6]);
+       copy_v3_v3(r->min, ml->bb->vec[0]);
+       r->ml = ml;
+}
+
+/**
+ * Partitions part of mainb array [start, end) along axis s. Returns i,
+ * where centroids of elements in the [start, i) segment lie "on the right side" of div,
+ * and elements in the [i, end) segment lie "on the left"
+ */
+static unsigned int partition_mainb(MetaElem **mainb, unsigned int start, unsigned int end, unsigned int s, float div)
+{
+       unsigned int i = start, j = end - 1;
+       div *= 2.0f;
+
+       while (1) {
+               while (i < j && div > (mainb[i]->bb->vec[6][s] + mainb[i]->bb->vec[0][s])) i++;
+               while (j > i && div < (mainb[j]->bb->vec[6][s] + mainb[j]->bb->vec[0][s])) j--;
+
+               if (i >= j)
+                       break;
+
+               SWAP(MetaElem *, mainb[i], mainb[j]);
+               i++;
+               j--;
+       }
+
+       if (i == start) {
+               i++;
+       }
+
+       return i;
+}
+
+/**
+ * Recursively builds a BVH, dividing elements along the middle of the longest axis of allbox.
+ */
+static void build_bvh_spatial(
+        PROCESS *process, MetaballBVHNode *node,
+        unsigned int start, unsigned int end, const Box *allbox)
+{
+       unsigned int part, j, s;
+       float dim[3], div;
+
+       /* Maximum bvh queue size is number of nodes which are made, equals calls to this function. */
+       process->bvh_queue_size++;
 
+       dim[0] = allbox->max[0] - allbox->min[0];
+       dim[1] = allbox->max[1] - allbox->min[1];
+       dim[2] = allbox->max[2] - allbox->min[2];
+
+       s = 0;
+       if (dim[1] > dim[0] && dim[1] > dim[2]) s = 1;
+       else if (dim[2] > dim[1] && dim[2] > dim[0]) s = 2;
+
+       div = allbox->min[s] + (dim[s] / 2.0f);
+
+       part = partition_mainb(process->mainb, start, end, s, div);
+
+       make_box_from_metaelem(&node->bb[0], process->mainb[start]);
+       node->child[0] = NULL;
+
+       if (part > start + 1) {
+               for (j = start; j < part; j++) {
+                       make_box_union(process->mainb[j]->bb, &node->bb[0], &node->bb[0]);
+               }
+
+               node->child[0] = BLI_memarena_alloc(process->pgn_elements, sizeof(MetaballBVHNode));
+               build_bvh_spatial(process, node->child[0], start, part, &node->bb[0]);
+       }
+
+       node->child[1] = NULL;
+       if (part < end) {
+               make_box_from_metaelem(&node->bb[1], process->mainb[part]);
+
+               if (part < end - 1) {
+                       for (j = part; j < end; j++) {
+                               make_box_union(process->mainb[j]->bb, &node->bb[1], &node->bb[1]);
+                       }
+
+                       node->child[1] = BLI_memarena_alloc(process->pgn_elements, sizeof(MetaballBVHNode));
+                       build_bvh_spatial(process, node->child[1], part, end, &node->bb[1]);
+               }
+       }
+       else {
+               INIT_MINMAX(node->bb[1].min, node->bb[1].max);
+       }
+}
 
 /* ******************** ARITH ************************* */
 
-/* BASED AT CODE (but mostly rewritten) :
+/**
+ * BASED AT CODE (but mostly rewritten) :
  * C code from the article
  * "An Implicit Surface Polygonizer"
  * by Jules Bloomenthal, jbloom@beauty.gmu.edu
@@ -178,9 +260,8 @@ static void converge(PROCESS *process, const float p1[3], const float p2[3], flo
  * 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    */
+ * any and all purposes, provided that this notice appears in all copies.
+ */
 
 #define L   0  /* left direction:      -x, -i */
 #define R   1  /* right direction:     +x, +i */
@@ -197,8 +278,10 @@ static void converge(PROCESS *process, const float p1[3], const float p2[3], flo
 #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) */
+/**
+ * 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) */
@@ -206,19 +289,19 @@ static void converge(PROCESS *process, const float p1[3], const float p2[3], flo
 #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 */
+// #define FLIP(i, bit) ((i) ^ 1 << (bit)) /* flip the given bit of i */
 
+/* ******************** DENSITY COPMPUTATION ********************* */
 
-/* **************** 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)
+/**
+ * Computes density from given metaball at given position.
+ * Metaball equation is: ``(1 - r^2 / R^2)^3 * s``
+ *
+ * r = distance from center
+ * R = metaball radius
+ * s - metaball stiffness
+ */
+static float densfunc(const MetaElem *ball, float x, float y, float z)
 {
        float dist2;
        float dvec[3] = {x, y, z};
@@ -229,37 +312,26 @@ static float densfunc(MetaElem *ball, float x, float y, float z)
                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_CUBE:
+                       if      (dvec[2] > ball->expz)  dvec[2] -= ball->expz;
+                       else if (dvec[2] < -ball->expz) dvec[2] += ball->expz;
+                       else                            dvec[2] = 0.0;
+                       /* fall through */
                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;
+                       /* fall through */
+               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_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:
@@ -277,208 +349,107 @@ static float densfunc(MetaElem *ball, float x, float y, float z)
                        else if (dvec[2] < -ball->len) dvec[2] += ball->len;
                        else                           dvec[2] = 0.0;
                        break;
-               /* *** end deprecated *** */
+                       /* *** end deprecated *** */
        }
 
-       dist2 = 1.0f - (len_squared_v3(dvec) / ball->rad2);
+       /* ball->rad2 is inverse of squared rad */
+       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);
+       /* ball->s is negative if metaball is negative */
+       return (dist2 < 0.0f) ? 0.0f : (ball->s * dist2 * dist2 * dist2);
 }
 
+/**
+ * Computes density at given position form all metaballs which contain this point in their box.
+ * Traverses BVH using a queue.
+ */
 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);
+       int i;
+       float dens = 0.0f;
+       unsigned int front = 0, back = 0;
+       MetaballBVHNode *node;
+
+       process->bvh_queue[front++] = &process->metaball_bvh;
+
+       while (front != back) {
+               node = process->bvh_queue[back++];
+
+               for (i = 0; i < 2; i++) {
+                       if ((node->bb[i].min[0] <= x) && (node->bb[i].max[0] >= x) &&
+                           (node->bb[i].min[1] <= y) && (node->bb[i].max[1] >= y) &&
+                           (node->bb[i].min[2] <= z) && (node->bb[i].max[2] >= z))
+                       {
+                               if (node->child[i])     process->bvh_queue[front++] = node->child[i];
+                               else dens += densfunc(node->bb[i].ml, 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)
+/**
+ * Adds face to indices, expands memory if needed.
+ */
+static void make_face(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;
+       int *cur;
+
+#ifdef USE_ACCUM_NORMAL
+       float n[3];
+#endif
+
+       if (UNLIKELY(process->totindex == process->curindex)) {
+               process->totindex += 4096;
+               process->indices = MEM_reallocN(process->indices, sizeof(int[4]) * process->totindex);
        }
-       
-       cur = process->indices + 4 * process->curindex;
+
+       cur = process->indices[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)
+
+       if (i4 == 0) {
                cur[3] = i3;
-       else 
+       }
+       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;
+
+#ifdef USE_ACCUM_NORMAL
+       if (i4 == 0) {
+               normal_tri_v3(n, process->co[i1], process->co[i2], process->co[i3]);
+               accumulate_vertex_normals(
+                       process->no[i1], process->no[i2], process->no[i3], NULL, n,
+                       process->co[i1], process->co[i2], process->co[i3], 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;
-               }
+       else {
+               normal_quad_v3(n, process->co[i1], process->co[i2], process->co[i3], process->co[i4]);
+               accumulate_vertex_normals(
+                       process->no[i1], process->no[i2], process->no[i3], process->no[i4], n,
+                       process->co[i1], process->co[i2], process->co[i3], process->co[i4]);
        }
-       
-       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;
+#endif
+
 }
 
+/* Frees allocated memory */
 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);
-       }
+       if (process->corners) MEM_freeN(process->corners);
+       if (process->edges) MEM_freeN(process->edges);
+       if (process->centers) MEM_freeN(process->centers);
+       if (process->mainb) MEM_freeN(process->mainb);
+       if (process->bvh_queue) MEM_freeN(process->bvh_queue);
+       if (process->pgn_elements) BLI_memarena_free(process->pgn_elements);
 }
 
+/* **************** POLYGONIZATION ************************ */
+
 /**** Cubical Polygonization (optional) ****/
 
 #define LB  0  /* left bottom edge     */
@@ -495,6 +466,7 @@ static void freepolygonize(PROCESS *process)
 #define TF  11 /* top far edge */
 
 static INTLISTS *cubetable[256];
+static char faces[256];
 
 /* edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
 static int corner1[12] = {
@@ -504,77 +476,89 @@ 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
+       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
+       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)
+/**
+ * triangulate the cube directly, without decomposition
+ */
+static void docube(PROCESS *process, CUBE *cube)
 {
        INTLISTS *polys;
        CORNER *c1, *c2;
        int i, index = 0, count, indexar[8];
-       
+
+       /* Determine which case cube falls into. */
        for (i = 0; i < 8; i++) {
                if (cube->corners[i]->value > 0.0f) {
                        index += (1 << i);
                }
        }
-       
+
+       /* Using faces[] table, adds neighbouring cube if surface intersects face in this direction. */
+       if (MB_BIT(faces[index], 0)) add_cube(process, cube->i - 1, cube->j, cube->k);
+       if (MB_BIT(faces[index], 1)) add_cube(process, cube->i + 1, cube->j, cube->k);
+       if (MB_BIT(faces[index], 2)) add_cube(process, cube->i, cube->j - 1, cube->k);
+       if (MB_BIT(faces[index], 3)) add_cube(process, cube->i, cube->j + 1, cube->k);
+       if (MB_BIT(faces[index], 4)) add_cube(process, cube->i, cube->j, cube->k - 1);
+       if (MB_BIT(faces[index], 5)) add_cube(process, cube->i, cube->j, cube->k + 1);
+
+       /* Using cubetable[], determines polygons for output. */
        for (polys = cubetable[index]; polys; polys = polys->next) {
                INTLIST *edges;
-               
+
                count = 0;
-               
+               /* Sets needed vertex id's lying on the edges. */
                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);
+
+                       indexar[count] = vertid(process, c1, c2);
                        count++;
                }
+
+               /* Adds faces to output. */
                if (count > 2) {
                        switch (count) {
                                case 3:
-                                       accum_mballfaces(process, indexar[2], indexar[1], indexar[0], 0);
+                                       make_face(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]);
+                                       if (indexar[0] == 0) make_face(process, indexar[0], indexar[3], indexar[2], indexar[1]);
+                                       else make_face(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);
+                                       if (indexar[0] == 0) make_face(process, indexar[0], indexar[3], indexar[2], indexar[1]);
+                                       else make_face(process, indexar[3], indexar[2], indexar[1], indexar[0]);
+
+                                       make_face(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]);
+                                               make_face(process, indexar[0], indexar[3], indexar[2], indexar[1]);
+                                               make_face(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]);
+                                               make_face(process, indexar[3], indexar[2], indexar[1], indexar[0]);
+                                               make_face(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]);
+                                               make_face(process, indexar[0], indexar[3], indexar[2], indexar[1]);
+                                               make_face(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]);
+                                               make_face(process, indexar[3], indexar[2], indexar[1], indexar[0]);
+                                               make_face(process, indexar[5], indexar[4], indexar[3], indexar[0]);
                                        }
-                               
-                                       accum_mballfaces(process, indexar[6], indexar[5], indexar[0], 0);
+
+                                       make_face(process, indexar[6], indexar[5], indexar[0], 0);
 
                                        break;
                        }
@@ -582,63 +566,10 @@ static void docube(PROCESS *process, CUBE *cube, MetaBall *mb)
        }
 }
 
-
-/* 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 */
-
+/**
+ * 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 */
@@ -648,32 +579,33 @@ static CORNER *setcorner(PROCESS *process, int i, int j, int k)
        /* 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 = BLI_memarena_alloc(process->pgn_elements, sizeof(CORNER));
 
-       c->i = i; 
+       c->i = i;
        c->co[0] = ((float)i - 0.5f) * process->size;
-       c->j = j; 
+       c->j = j;
        c->co[1] = ((float)j - 0.5f) * process->size;
-       c->k = k; 
+       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->value = metaball(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 */
-
+/**
+ * return next clockwise edge from given edge around given face
+ */
 static int nextcwedge(int edge, int face)
 {
        switch (edge) {
@@ -705,18 +637,18 @@ static int nextcwedge(int edge, int face)
        return 0;
 }
 
-
-/* otherface: return face adjoining edge that is not the given face */
-
+/**
+ * \return the 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 */
-
+/**
+ * create the 256 entry table for cubical polygonization
+ */
 static void makecubetable(void)
 {
        static bool is_done = false;
@@ -731,22 +663,22 @@ static void makecubetable(void)
                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");
+                               INTLISTS *lists = 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 = 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);
                                        }
@@ -756,6 +688,23 @@ static void makecubetable(void)
                                cubetable[i] = lists;
                        }
        }
+
+       for (i = 0; i < 256; i++) {
+               INTLISTS *polys;
+               faces[i] = 0;
+               for (polys = cubetable[i]; polys; polys = polys->next) {
+                       INTLIST *edges;
+
+                       for (edges = polys->list; edges; edges = edges->next) {
+                               if (edges->i == LB || edges->i == LT || edges->i == LN || edges->i == LF) faces[i] |= 1 << L;
+                               if (edges->i == RB || edges->i == RT || edges->i == RN || edges->i == RF) faces[i] |= 1 << R;
+                               if (edges->i == LB || edges->i == RB || edges->i == BN || edges->i == BF) faces[i] |= 1 << B;
+                               if (edges->i == LT || edges->i == RT || edges->i == TN || edges->i == TF) faces[i] |= 1 << T;
+                               if (edges->i == LN || edges->i == RN || edges->i == BN || edges->i == TN) faces[i] |= 1 << N;
+                               if (edges->i == LF || edges->i == RF || edges->i == BF || edges->i == TF) faces[i] |= 1 << F;
+                       }
+               }
+       }
 }
 
 void BKE_mball_cubeTable_free(void)
@@ -768,14 +717,14 @@ void BKE_mball_cubeTable_free(void)
                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;
                }
@@ -785,9 +734,9 @@ void BKE_mball_cubeTable_free(void)
 
 /**** Storage ****/
 
-/* setcenter: set (i, j, k) entry of table[]
- * return 1 if already set; otherwise, set and return 0 */
-
+/**
+ * Inserts cube at lattice i, j, k into hash table, marking it as "done"
+ */
 static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const int j, const int k)
 {
        int index;
@@ -799,30 +748,29 @@ static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const i
        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 = BLI_memarena_alloc(process->pgn_elements, 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)
+/**
+ * Sets vid of vertex lying on given edge.
+ */
+static void setedge(
+        PROCESS *process,
+        int i1, int j1, int k1,
+        int i2, int j2, int k2,
+        int vid)
 {
-       unsigned int index;
+       int index;
        EDGELIST *newe;
-       
+
        if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
                int t = i1;
                i1 = i2;
@@ -835,27 +783,28 @@ static void setedge(PROCESS *process,
                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 = BLI_memarena_alloc(process->pgn_elements, sizeof(EDGELIST));
+
+       newe->i1 = i1;
+       newe->j1 = j1;
        newe->k1 = k1;
-       newe->i2 = i2; 
-       newe->j2 = j2; 
+       newe->i2 = i2;
+       newe->j2 = j2;
        newe->k2 = k2;
        newe->vid = vid;
-       newe->next = table[index];
-       table[index] = newe;
+       newe->next = process->edges[index];
+       process->edges[index] = newe;
 }
 
-
-/* getedge: return vertex id for edge; return -1 if not set */
-
+/**
+ * \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;
@@ -878,62 +827,52 @@ static int getedge(EDGELIST *table[],
        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)
+/**
+ * Adds a vertex, expands memory if needed.
+ */
+static void addtovertices(PROCESS *process, const float v[3], const float no[3])
 {
-       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;
+       if (process->curvertex == process->totvertex) {
+               process->totvertex += 4096;
+               process->co = MEM_reallocN(process->co, process->totvertex * sizeof(float[3]));
+               process->no = MEM_reallocN(process->no, process->totvertex * sizeof(float[3]));
        }
-       vertices->ptr[vertices->count++] = v;
-}
 
-/* vnormal: compute unit length surface normal at point */
+       copy_v3_v3(process->co[process->curvertex], v);
+       copy_v3_v3(process->no[process->curvertex], no);
 
+       process->curvertex++;
+}
+
+#ifndef USE_ACCUM_NORMAL
+/**
+ * Computes normal from density field at given point.
+ *
+ * \note Doesn't do normalization!
+ */
 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]);
+       const float delta = process->delta;
+       const float f = metaball(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;
+       r_no[0] = metaball(process, point[0] + delta, point[1], point[2]) - f;
+       r_no[1] = metaball(process, point[0], point[1] + delta, point[2]) - f;
+       r_no[2] = metaball(process, point[0], point[1], point[2] + delta) - f;
 
-#if 1
-       normalize_v3(r_no);
-#else
+#if 0
        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);
@@ -941,355 +880,242 @@ static void vnormal(PROCESS *process, const float point[3], float r_no[3])
        }
 #endif
 }
+#endif  /* USE_ACCUM_NORMAL */
 
-
-static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2, MetaBall *mb)
+/**
+ * \return the id of vertex between two corners.
+ *
+ * If it wasn't previously computed, does #converge() and adds vertex to process.
+ */
+static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2)
 {
-       VERTEX v;
+       float v[3], no[3];
        int vid = getedge(process->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
 
-       if (vid != -1) {
-               return vid;  /* previously computed */
-       }
+       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);
+       converge(process, c1, c2, v);  /* position */
+
+#ifdef USE_ACCUM_NORMAL
+       zero_v3(no);
+#else
+       vnormal(process, v, no);
+#endif
+
+       addtovertices(process, v, no);            /* save vertex */
+       vid = (int)process->curvertex - 1;
+       setedge(process, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
 
-       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)
+/**
+ * Given two corners, computes approximation of surface intersection point between them.
+ * In case of small threshold, do bisection.
+ */
+static void converge(PROCESS *process, const CORNER *c1, const CORNER *c2, float r_p[3])
 {
-       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;
+       float tmp, dens;
+       unsigned int i;
+       float c1_value, c1_co[3];
+       float c2_value, c2_co[3];
+
+       if (c1->value < c2->value) {
+               c1_value = c2->value;
+               copy_v3_v3(c1_co, c2->co);
+               c2_value = c1->value;
+               copy_v3_v3(c2_co, c1->co);
        }
        else {
-               copy_v3_v3(pos, p1);
-               copy_v3_v3(neg, p2);
-               positive = v1;
-               negative = v2;
+               c1_value = c1->value;
+               copy_v3_v3(c1_co, c1->co);
+               c2_value = c2->value;
+               copy_v3_v3(c2_co, c2->co);
        }
 
-       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;
-               }
-       }
+       for (i = 0; i < process->converge_res; i++) {
+               interp_v3_v3v3(r_p, c1_co, c2_co, 0.5f);
+               dens = metaball(process, r_p[0], r_p[1], r_p[2]);
 
-       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 (dens > 0.0f) {
+                       c1_value = dens;
+                       copy_v3_v3(c1_co, r_p);
                }
-       }
-
-       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];
+               else {
+                       c2_value = dens;
+                       copy_v3_v3(c2_co, r_p);
                }
        }
 
-       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];
-               }
-       }
+       tmp = -c1_value / (c2_value - c1_value);
+       interp_v3_v3v3(r_p, c1_co, c2_co, tmp);
+}
 
-       /* This is necessary to find start point */
-       while (1) {
-               mid_v3_v3v3(&p[0], pos, neg);
+/**
+ * Adds cube at given lattice position to cube stack of process.
+ */
+static void add_cube(PROCESS *process, int i, int j, int k)
+{
+       CUBES *ncube;
+       int n;
 
-               if (i++ == RES) {
-                       return;
-               }
+       /* test if cube has been found before */
+       if (setcenter(process, process->centers, i, j, k) == 0) {
+               /* push cube on stack: */
+               ncube = BLI_memarena_alloc(process->pgn_elements, sizeof(CUBES));
+               ncube->next = process->cubes;
+               process->cubes = ncube;
 
-               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]);
-               }
+               ncube->cube.i = i;
+               ncube->cube.j = j;
+               ncube->cube.k = k;
+
+               /* set corners of initial cube: */
+               for (n = 0; n < 8; n++)
+                       ncube->cube.corners[n] = setcorner(process, i + MB_BIT(n, 2), j + MB_BIT(n, 1), k + MB_BIT(n, 0));
        }
 }
 
-/* ************************************** */
-static void add_cube(PROCESS *process, int i, int j, int k, int count)
+static void next_lattice(int r[3], const float pos[3], const float size)
 {
-       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));
-                               }
-                       }
+       r[0] = (int)ceil((pos[0] / size) + 0.5f);
+       r[1] = (int)ceil((pos[1] / size) + 0.5f);
+       r[2] = (int)ceil((pos[2] / size) + 0.5f);
 }
-
-
-static void find_first_points(PROCESS *process, MetaBall *mb, int a)
+static void prev_lattice(int r[3], const float pos[3], const float size)
 {
-       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;
-                       }
+       next_lattice(r, pos, size);
+       r[0]--; r[1]--; r[2]--;
+}
+static void closest_latice(int r[3], const float pos[3], const float size)
+{
+       r[0] = (int)floorf(pos[0] / size + 1.0f);
+       r[1] = (int)floorf(pos[1] / size + 1.0f);
+       r[2] = (int)floorf(pos[2] / size + 1.0f);
+}
 
-                       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;
+/**
+ * Find at most 26 cubes to start polygonization from.
+ */
+static void find_first_points(PROCESS *process, const unsigned int em)
+{
+       const MetaElem *ml;
+       int center[3], lbn[3], rtf[3], it[3], dir[3], add[3];
+       float tmp[3], a, b;
+
+       ml = process->mainb[em];
+
+       mid_v3_v3v3(tmp, ml->bb->vec[0], ml->bb->vec[6]);
+       closest_latice(center, tmp, process->size);
+       prev_lattice(lbn, ml->bb->vec[0], process->size);
+       next_lattice(rtf, ml->bb->vec[6], process->size);
+
+       for (dir[0] = -1; dir[0] <= 1; dir[0]++) {
+               for (dir[1] = -1; dir[1] <= 1; dir[1]++) {
+                       for (dir[2] = -1; dir[2] <= 1; dir[2]++) {
+                               if (dir[0] == 0 && dir[1] == 0 && dir[2] == 0) {
+                                       continue;
                                }
-                       
-                               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;
 
-                                               }
+                               copy_v3_v3_int(it, center);
+
+                               b = setcorner(process, it[0], it[1], it[2])->value;
+                               do {
+                                       it[0] += dir[0];
+                                       it[1] += dir[1];
+                                       it[2] += dir[2];
+                                       a = b;
+                                       b = setcorner(process, it[0], it[1], it[2])->value;
+
+                                       if (a * b < 0.0f) {
+                                               add[0] = it[0] - dir[0];
+                                               add[1] = it[1] - dir[1];
+                                               add[2] = it[2] - dir[2];
+                                               DO_MIN(it, add);
+                                               add_cube(process, add[0], add[1], add[2]);
+                                               break;
                                        }
-                               }
+                               } while ((it[0] > lbn[0]) && (it[1] > lbn[1]) && (it[2] > lbn[2]) &&
+                                        (it[0] < rtf[0]) && (it[1] < rtf[1]) && (it[2] < rtf[2]));
                        }
                }
        }
 }
 
-static void polygonize(PROCESS *process, MetaBall *mb)
+/**
+ * The main polygonization proc.
+ * Allocates memory, makes cubetable,
+ * finds starting surface points
+ * and processes cubes on the stack until none left.
+ */
+static void polygonize(PROCESS *process)
 {
        CUBE c;
-       int a;
-
-       process->vertices.count = process->vertices.max = 0;
-       process->vertices.ptr = NULL;
+       unsigned int i;
 
-       /* 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();
+       process->bvh_queue = MEM_callocN(sizeof(MetaballBVHNode *) * process->bvh_queue_size, "Metaball BVH Queue");
 
-       for (a = 0; a < process->totelem; a++) {
+       makecubetable();
 
-               /* try to find 8 points on the surface for each MetaElem */
-               find_first_points(process, mb, a);
+       for (i = 0; i < process->totelem; i++) {
+               find_first_points(process, i);
        }
 
-       /* polygonize all MetaElems of current MetaBall */
-       while (process->cubes != NULL) { /* process active cubes till none left */
+       while (process->cubes != NULL) {
                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);
+
+               docube(process, &c);
        }
 }
 
-static float init_meta(EvaluationContext *eval_ctx, PROCESS *process, Scene *scene, Object *ob)    /* return totsize */
+/**
+ * Iterates over ALL objects in the scene and all of its sets, including
+ * making all duplis(not only metas). Copies metas to mainb array.
+ * Computes bounding boxes for building BVH. */
+static void init_meta(EvaluationContext *eval_ctx, PROCESS *process, Scene *scene, Object *ob)
 {
        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;
+       const MetaElem *ml;
+       float obinv[4][4], obmat[4][4];
+       unsigned int i;
+       int 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;
                                }
@@ -1312,562 +1138,124 @@ static float init_meta(EvaluationContext *eval_ctx, PROCESS *process, Scene *sce
                        }
 
                        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 pos[4][4], rot[4][4];
                                                float expx, expy, expz;
+                                               float tempmin[3], tempmax[3];
 
-                                               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);
+                                               MetaElem *new_ml;
 
                                                /* 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));
+                                               new_ml = BLI_memarena_alloc(process->pgn_elements, sizeof(MetaElem));
+                                               *(new_ml) = *ml;
+                                               new_ml->bb = BLI_memarena_alloc(process->pgn_elements, sizeof(BoundBox));
+                                               new_ml->mat = BLI_memarena_alloc(process->pgn_elements, 4 * 4 * sizeof(float));
+                                               new_ml->imat = BLI_memarena_alloc(process->pgn_elements, 4 * 4 * sizeof(float));
 
-                                               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);
+                                               /* too big stiffness seems only ugly due to linear interpolation
+                                               * no need to have possibility for too big stiffness */
+                                               if (ml->s > 10.0f) new_ml->s = 10.0f;
+                                               else new_ml->s = ml->s;
 
-                                               process->mainb[a]->rad2 = ml->rad * ml->rad;
+                                               /* if metaball is negative, set stiffness negative */
+                                               if (new_ml->flag & MB_NEGATIVE) new_ml->s = -new_ml->s;
 
-                                               process->mainb[a]->mat = (float *) mat;
-                                               process->mainb[a]->imat = (float *) imat;
+                                               /* Translation of MetaElem */
+                                               unit_m4(pos);
+                                               pos[3][0] = ml->x;
+                                               pos[3][1] = ml->y;
+                                               pos[3][2] = ml->z;
 
-                                               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;
+                                               /* Rotation of MetaElem is stored in quat */
+                                               quat_to_mat4(rot, ml->quat);
+
+                                               /* basis object space -> world -> ml object space -> position -> rotation -> ml local space */
+                                               mul_m4_series((float(*)[4])new_ml->mat, obinv, bob->obmat, pos, rot);
+                                               /* ml local space -> basis object space */
+                                               invert_m4_m4((float(*)[4])new_ml->imat, (float(*)[4])new_ml->mat);
+
+                                               /* rad2 is inverse of squared radius */
+                                               new_ml->rad2 = 1 / (ml->rad * ml->rad);
+
+                                               /* initial dimensions = radius */
+                                               expx = ml->rad;
+                                               expy = ml->rad;
+                                               expz = ml->rad;
+
+                                               switch (ml->type) {
+                                                       case MB_BALL:
+                                                               break;
+                                                       case MB_CUBE: /* cube is "expanded" by expz, expy and expx */
+                                                               expz += ml->expz;
+                                                               /* fall through */
+                                                       case MB_PLANE: /* plane is "expanded" by expy and expx */
+                                                               expy += ml->expy;
+                                                               /* fall through */
+                                                       case MB_TUBE: /* tube is "expanded" by expx */
+                                                               expx += ml->expx;
+                                                               break;
+                                                       case MB_ELIPSOID: /* ellipsoid is "stretched" by exp* */
+                                                               expx *= ml->expx;
+                                                               expy *= ml->expy;
+                                                               expz *= ml->expz;
+                                                               break;
                                                }
 
                                                /* 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 */
+                                               copy_v3_fl3(new_ml->bb->vec[0], -expx, -expy, -expz);  /* 0 */
+                                               copy_v3_fl3(new_ml->bb->vec[1], +expx, -expy, -expz);  /* 1 */
+                                               copy_v3_fl3(new_ml->bb->vec[2], +expx, +expy, -expz);  /* 2 */
+                                               copy_v3_fl3(new_ml->bb->vec[3], -expx, +expy, -expz);  /* 3 */
+                                               copy_v3_fl3(new_ml->bb->vec[4], -expx, -expy, +expz);  /* 4 */
+                                               copy_v3_fl3(new_ml->bb->vec[5], +expx, -expy, +expz);  /* 5 */
+                                               copy_v3_fl3(new_ml->bb->vec[6], +expx, +expy, +expz);  /* 6 */
+                                               copy_v3_fl3(new_ml->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]);
+                                                       mul_m4_v3((float(*)[4])new_ml->mat, new_ml->bb->vec[i]);
 
                                                /* find max and min of transformed bb */
+                                               INIT_MINMAX(tempmin, tempmax);
                                                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];
+                                                       DO_MINMAX(new_ml->bb->vec[i], tempmin, tempmax);
                                                }
-                                       
-                                               /* 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;
-       }
+                                               /* set only point 0 and 6 - AABB of Metaelem */
+                                               copy_v3_v3(new_ml->bb->vec[0], tempmin);
+                                               copy_v3_v3(new_ml->bb->vec[6], tempmax);
 
-       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);
-                                               }
-                                       }
+                                               /* add new_ml to mainb[] */
+                                               if (process->totelem == process->mem) {
+                                                       MetaElem **newelem;
+                                                       process->mem = process->mem * 2 + 10;
+                                                       newelem = MEM_mallocN(sizeof(MetaElem *) * process->mem, "metaballs");
 
-                                       /* 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);
+                                                       memcpy(newelem, process->mainb, sizeof(MetaElem *) * process->totelem);
+                                                       if (process->mainb) MEM_freeN(process->mainb);
+                                                       process->mainb = newelem;
                                                }
-                                               
-                                       }
-                       
-                                       /* 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);
+                                               process->mainb[process->totelem++] = new_ml;
                                        }
+                                       ml = ml->next;
                                }
-
                        }
-                       
-                       /* 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++;
-                               }
-                       }
-               }
+       /* compute AABB of all Metaelems */
+       if (process->totelem > 0) {
+               copy_v3_v3(process->allbb.min, process->mainb[0]->bb->vec[0]);
+               copy_v3_v3(process->allbb.max, process->mainb[0]->bb->vec[6]);
+               for (i = 1; i < process->totelem; i++)
+                       make_box_union(process->mainb[i]->bb, &process->allbb, &process->allbb);
        }
 }
 
@@ -1875,102 +1263,66 @@ void BKE_mball_polygonize(EvaluationContext *eval_ctx, Scene *scene, Object *ob,
 {
        MetaBall *mb;
        DispList *dl;
-       int a, nr_cubes;
-       float *co, *no, totsize, width;
+       unsigned int a;
        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);
+       if      (process.thresh < 0.001f) process.converge_res = 16;
+       else if (process.thresh < 0.01f)  process.converge_res = 8;
+       else if (process.thresh < 0.1f)   process.converge_res = 4;
+       else                              process.converge_res = 2;
 
-                       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;
 
-       /* width is size per polygonize cube */
        if (eval_ctx->mode == DAG_EVAL_RENDER) {
-               width = mb->rendersize;
+               process.size = mb->rendersize;
        }
        else {
-               width = mb->wiresize;
+               process.size = mb->wiresize;
                if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_HALFRES) {
-                       width *= 2;
+                       process.size *= 2.0f;
                }
        }
-       /* 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;
+       process.delta = process.size * 0.001f;
 
-               dl = MEM_callocN(sizeof(DispList), "mbaldisp");
-               BLI_addtail(dispbase, dl);
-               dl->type = DL_INDEX4;
-               dl->nr = process.vertices.count;
-               dl->parts = process.curindex;
+       process.pgn_elements = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "Metaball memarena");
 
-               dl->index = process.indices;
-               process.indices = NULL;
+       /* initialize all mainb (MetaElems) */
+       init_meta(eval_ctx, &process, scene, ob);
 
-               a = process.vertices.count;
-               dl->verts = co = MEM_mallocN(sizeof(float[3]) * a, "mballverts");
-               dl->nors = no = MEM_mallocN(sizeof(float[3]) * a, "mballnors");
+       if (process.totelem > 0) {
+               build_bvh_spatial(&process, &process.metaball_bvh, 0, process.totelem, &process.allbb);
 
-               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);
+               /* 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 (ob->size[0] > 0.00001f * (process.allbb.max[0] - process.allbb.min[0]) ||
+                   ob->size[1] > 0.00001f * (process.allbb.max[1] - process.allbb.min[1]) ||
+                   ob->size[2] > 0.00001f * (process.allbb.max[2] - process.allbb.min[2]))
+               {
+                       polygonize(&process);
+
+                       /* add resulting surface to displist */
+                       if (process.curindex) {
+                               dl = MEM_callocN(sizeof(DispList), "mballdisp");
+                               BLI_addtail(dispbase, dl);
+                               dl->type = DL_INDEX4;
+                               dl->nr = (int)process.curvertex;
+                               dl->parts = (int)process.curindex;
+
+                               dl->index = (int *)process.indices;
+
+                               for (a = 0; a < process.curvertex; a++) {
+                                       normalize_v3(process.no[a]);
+                               }
+
+                               dl->verts = (float *)process.co;
+                               dl->nors = (float *)process.no;
+                       }
                }
        }
 
index f1f73563decb87071e7a654d6250443d4a24b8b5..4f019c0a913a96957e8d7bebc504edb8320d4062 100644 (file)
@@ -60,7 +60,7 @@
 #include "BKE_global.h"
 #include "BKE_library.h"
 #include "BKE_main.h"
-#include "BKE_mball.h"
+#include "BKE_mball_tessellate.h"
 #include "BKE_node.h"
 #include "BKE_report.h"