Threaded object update and EvaluationContext
[blender.git] / source / blender / blenkernel / intern / mball.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19  * All rights reserved.
20  *
21  * Contributor(s): Jiri Hnidek <jiri.hnidek@vslib.cz>.
22  *
23  * ***** END GPL LICENSE BLOCK *****
24  *
25  * MetaBalls are created from a single Object (with a name without number in it),
26  * here the DispList and BoundBox also is located.
27  * All objects with the same name (but with a number in it) are added to this.
28  *
29  * texture coordinates are patched within the displist
30  */
31
32 /** \file blender/blenkernel/intern/mball.c
33  *  \ingroup bke
34  */
35
36 #include <stdio.h>
37 #include <string.h>
38 #include <math.h>
39 #include <stdlib.h>
40 #include <ctype.h>
41 #include <float.h>
42 #include "MEM_guardedalloc.h"
43
44 #include "DNA_material_types.h"
45 #include "DNA_object_types.h"
46 #include "DNA_meta_types.h"
47 #include "DNA_scene_types.h"
48
49
50 #include "BLI_blenlib.h"
51 #include "BLI_math.h"
52 #include "BLI_utildefines.h"
53
54 #include "BKE_global.h"
55 #include "BKE_main.h"
56
57 /*  #include "BKE_object.h" */
58 #include "BKE_animsys.h"
59 #include "BKE_curve.h"
60 #include "BKE_depsgraph.h"
61 #include "BKE_scene.h"
62 #include "BKE_library.h"
63 #include "BKE_displist.h"
64 #include "BKE_mball.h"
65 #include "BKE_object.h"
66 #include "BKE_material.h"
67
68 /* Data types */
69
70 typedef struct vertex {         /* surface vertex */
71         float co[3];  /* position and surface normal */
72         float no[3];
73 } VERTEX;
74
75 typedef struct vertices {       /* list of vertices in polygonization */
76         int count, max;             /* # vertices, max # allowed */
77         VERTEX *ptr;                /* dynamically allocated */
78 } VERTICES;
79
80 typedef struct corner {         /* corner of a cube */
81         int i, j, k;                /* (i, j, k) is index within lattice */
82         float co[3], value;       /* location and function value */
83         struct corner *next;
84 } CORNER;
85
86 typedef struct cube {           /* partitioning cell (cube) */
87         int i, j, k;                /* lattice location of cube */
88         CORNER *corners[8];         /* eight corners */
89 } CUBE;
90
91 typedef struct cubes {          /* linked list of cubes acting as stack */
92         CUBE cube;                  /* a single cube */
93         struct cubes *next;         /* remaining elements */
94 } CUBES;
95
96 typedef struct centerlist {     /* list of cube locations */
97         int i, j, k;                /* cube location */
98         struct centerlist *next;    /* remaining elements */
99 } CENTERLIST;
100
101 typedef struct edgelist {       /* list of edges */
102         int i1, j1, k1, i2, j2, k2; /* edge corner ids */
103         int vid;                    /* vertex id */
104         struct edgelist *next;      /* remaining elements */
105 } EDGELIST;
106
107 typedef struct intlist {        /* list of integers */
108         int i;                      /* an integer */
109         struct intlist *next;       /* remaining elements */
110 } INTLIST;
111
112 typedef struct intlists {       /* list of list of integers */
113         INTLIST *list;              /* a list of integers */
114         struct intlists *next;      /* remaining elements */
115 } INTLISTS;
116
117 /* dividing scene using octal tree makes polygonisation faster */
118 typedef struct ml_pointer {
119         struct ml_pointer *next, *prev;
120         struct MetaElem *ml;
121 } ml_pointer;
122
123 typedef struct octal_node {
124         struct octal_node *nodes[8];/* children of current node */
125         struct octal_node *parent;  /* parent of current node */
126         struct ListBase elems;      /* ListBase of MetaElem pointers (ml_pointer) */
127         float x_min, y_min, z_min;  /* 1st border point */
128         float x_max, y_max, z_max;  /* 7th border point */
129         float x, y, z;              /* center of node */
130         int pos, neg;               /* number of positive and negative MetaElements in the node */
131         int count;                  /* number of MetaElems, which belongs to the node */
132 } octal_node;
133
134 typedef struct octal_tree {
135         struct octal_node *first;   /* first node */
136         int pos, neg;               /* number of positive and negative MetaElements in the scene */
137         short depth;                /* number of scene subdivision */
138 } octal_tree;
139
140 struct pgn_elements {
141         struct pgn_elements *next, *prev;
142         char *data;
143 };
144
145 typedef struct process {        /* parameters, function, storage */
146         /* ** old G_mb contents ** */
147         float thresh;
148         int totelem;
149         MetaElem **mainb;
150         octal_tree *metaball_tree;
151
152         /* ** old process contents ** */
153
154         /* what happens here? floats, I think. */
155         /*  float (*function)(void);     */     /* implicit surface function */
156         float (*function)(struct process *, float, float, float);
157         float size, delta;          /* cube size, normal delta */
158         int bounds;                 /* cube range within lattice */
159         CUBES *cubes;               /* active cubes */
160         VERTICES vertices;          /* surface vertices */
161         CENTERLIST **centers;       /* cube center hash table */
162         CORNER **corners;           /* corner value hash table */
163         EDGELIST **edges;           /* edge and vertex id hash table */
164
165         /* Runtime things */
166         int *indices;
167         int totindex, curindex;
168
169         int pgn_offset;
170         struct pgn_elements *pgn_current;
171         ListBase pgn_list;
172 } PROCESS;
173
174 /* Forward declarations */
175 static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2, MetaBall *mb);
176 static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const int j, const int k);
177 static CORNER *setcorner(PROCESS *process, int i, int j, int k);
178 static void converge(PROCESS *process, const float p1[3], const float p2[3], float v1, float v2,
179                      float p[3], MetaBall *mb, int f);
180
181 /* Functions */
182
183 void BKE_mball_unlink(MetaBall *mb)
184 {
185         int a;
186         
187         for (a = 0; a < mb->totcol; a++) {
188                 if (mb->mat[a]) mb->mat[a]->id.us--;
189                 mb->mat[a] = NULL;
190         }
191 }
192
193
194 /* do not free mball itself */
195 void BKE_mball_free(MetaBall *mb)
196 {
197         BKE_mball_unlink(mb);
198         
199         if (mb->adt) {
200                 BKE_free_animdata((ID *)mb);
201                 mb->adt = NULL;
202         }
203         if (mb->mat) MEM_freeN(mb->mat);
204         BLI_freelistN(&mb->elems);
205         if (mb->disp.first) BKE_displist_free(&mb->disp);
206 }
207
208 MetaBall *BKE_mball_add(Main *bmain, const char *name)
209 {
210         MetaBall *mb;
211         
212         mb = BKE_libblock_alloc(&bmain->mball, ID_MB, name);
213         
214         mb->size[0] = mb->size[1] = mb->size[2] = 1.0;
215         mb->texflag = MB_AUTOSPACE;
216         
217         mb->wiresize = 0.4f;
218         mb->rendersize = 0.2f;
219         mb->thresh = 0.6f;
220         
221         return mb;
222 }
223
224 MetaBall *BKE_mball_copy(MetaBall *mb)
225 {
226         MetaBall *mbn;
227         int a;
228         
229         mbn = BKE_libblock_copy(&mb->id);
230
231         BLI_duplicatelist(&mbn->elems, &mb->elems);
232         
233         mbn->mat = MEM_dupallocN(mb->mat);
234         for (a = 0; a < mbn->totcol; a++) {
235                 id_us_plus((ID *)mbn->mat[a]);
236         }
237
238         mbn->editelems = NULL;
239         mbn->lastelem = NULL;
240         
241         return mbn;
242 }
243
244 static void extern_local_mball(MetaBall *mb)
245 {
246         if (mb->mat) {
247                 extern_local_matarar(mb->mat, mb->totcol);
248         }
249 }
250
251 void BKE_mball_make_local(MetaBall *mb)
252 {
253         Main *bmain = G.main;
254         Object *ob;
255         int is_local = FALSE, is_lib = FALSE;
256
257         /* - only lib users: do nothing
258          * - only local users: set flag
259          * - mixed: make copy
260          */
261         
262         if (mb->id.lib == NULL) return;
263         if (mb->id.us == 1) {
264                 id_clear_lib_data(bmain, &mb->id);
265                 extern_local_mball(mb);
266                 
267                 return;
268         }
269
270         for (ob = G.main->object.first; ob && ELEM(0, is_lib, is_local); ob = ob->id.next) {
271                 if (ob->data == mb) {
272                         if (ob->id.lib) is_lib = TRUE;
273                         else is_local = TRUE;
274                 }
275         }
276         
277         if (is_local && is_lib == FALSE) {
278                 id_clear_lib_data(bmain, &mb->id);
279                 extern_local_mball(mb);
280         }
281         else if (is_local && is_lib) {
282                 MetaBall *mb_new = BKE_mball_copy(mb);
283                 mb_new->id.us = 0;
284
285                 /* Remap paths of new ID using old library as base. */
286                 BKE_id_lib_local_paths(bmain, mb->id.lib, &mb_new->id);
287
288                 for (ob = G.main->object.first; ob; ob = ob->id.next) {
289                         if (ob->data == mb) {
290                                 if (ob->id.lib == NULL) {
291                                         ob->data = mb_new;
292                                         mb_new->id.us++;
293                                         mb->id.us--;
294                                 }
295                         }
296                 }
297         }
298 }
299
300 /* most simple meta-element adding function
301  * don't do context manipulation here (rna uses) */
302 MetaElem *BKE_mball_element_add(MetaBall *mb, const int type)
303 {
304         MetaElem *ml = MEM_callocN(sizeof(MetaElem), "metaelem");
305
306         unit_qt(ml->quat);
307
308         ml->rad = 2.0;
309         ml->s = 2.0;
310         ml->flag = MB_SCALE_RAD;
311
312         switch (type) {
313                 case MB_BALL:
314                         ml->type = MB_BALL;
315                         ml->expx = ml->expy = ml->expz = 1.0;
316
317                         break;
318                 case MB_TUBE:
319                         ml->type = MB_TUBE;
320                         ml->expx = ml->expy = ml->expz = 1.0;
321
322                         break;
323                 case MB_PLANE:
324                         ml->type = MB_PLANE;
325                         ml->expx = ml->expy = ml->expz = 1.0;
326
327                         break;
328                 case MB_ELIPSOID:
329                         ml->type = MB_ELIPSOID;
330                         ml->expx = 1.2f;
331                         ml->expy = 0.8f;
332                         ml->expz = 1.0;
333
334                         break;
335                 case MB_CUBE:
336                         ml->type = MB_CUBE;
337                         ml->expx = ml->expy = ml->expz = 1.0;
338
339                         break;
340                 default:
341                         break;
342         }
343
344         BLI_addtail(&mb->elems, ml);
345
346         return ml;
347 }
348 /** Compute bounding box of all MetaElems/MetaBalls.
349  *
350  * Bounding box is computed from polygonized surface. Object *ob is
351  * basic MetaBall (usually with name Meta). All other MetaBalls (with
352  * names Meta.001, Meta.002, etc) are included in this Bounding Box.
353  */
354 void BKE_mball_texspace_calc(Object *ob)
355 {
356         DispList *dl;
357         BoundBox *bb;
358         float *data, min[3], max[3] /*, loc[3], size[3] */;
359         int tot, do_it = FALSE;
360
361         if (ob->bb == NULL) ob->bb = MEM_callocN(sizeof(BoundBox), "mb boundbox");
362         bb = ob->bb;
363         
364         /* Weird one, this. */
365 /*      INIT_MINMAX(min, max); */
366         (min)[0] = (min)[1] = (min)[2] = 1.0e30f;
367         (max)[0] = (max)[1] = (max)[2] = -1.0e30f;
368
369         dl = ob->curve_cache->disp.first;
370         while (dl) {
371                 tot = dl->nr;
372                 if (tot) do_it = TRUE;
373                 data = dl->verts;
374                 while (tot--) {
375                         /* Also weird... but longer. From utildefines. */
376                         minmax_v3v3_v3(min, max, data);
377                         data += 3;
378                 }
379                 dl = dl->next;
380         }
381
382         if (!do_it) {
383                 min[0] = min[1] = min[2] = -1.0f;
384                 max[0] = max[1] = max[2] = 1.0f;
385         }
386 #if 0
387         loc[0] = (min[0] + max[0]) / 2.0f;
388         loc[1] = (min[1] + max[1]) / 2.0f;
389         loc[2] = (min[2] + max[2]) / 2.0f;
390
391         size[0] = (max[0] - min[0]) / 2.0f;
392         size[1] = (max[1] - min[1]) / 2.0f;
393         size[2] = (max[2] - min[2]) / 2.0f;
394 #endif
395         BKE_boundbox_init_from_minmax(bb, min, max);
396 }
397
398 float *BKE_mball_make_orco(Object *ob, ListBase *dispbase)
399 {
400         BoundBox *bb;
401         DispList *dl;
402         float *data, *orco, *orcodata;
403         float loc[3], size[3];
404         int a;
405
406         /* restore size and loc */
407         bb = ob->bb;
408         loc[0] = (bb->vec[0][0] + bb->vec[4][0]) / 2.0f;
409         size[0] = bb->vec[4][0] - loc[0];
410         loc[1] = (bb->vec[0][1] + bb->vec[2][1]) / 2.0f;
411         size[1] = bb->vec[2][1] - loc[1];
412         loc[2] = (bb->vec[0][2] + bb->vec[1][2]) / 2.0f;
413         size[2] = bb->vec[1][2] - loc[2];
414
415         dl = dispbase->first;
416         orcodata = MEM_mallocN(sizeof(float) * 3 * dl->nr, "MballOrco");
417
418         data = dl->verts;
419         orco = orcodata;
420         a = dl->nr;
421         while (a--) {
422                 orco[0] = (data[0] - loc[0]) / size[0];
423                 orco[1] = (data[1] - loc[1]) / size[1];
424                 orco[2] = (data[2] - loc[2]) / size[2];
425
426                 data += 3;
427                 orco += 3;
428         }
429
430         return orcodata;
431 }
432
433 /* Note on mball basis stuff 2.5x (this is a can of worms)
434  * This really needs a rewrite/refactor its totally broken in anything other then basic cases
435  * Multiple Scenes + Set Scenes & mixing mball basis SHOULD work but fails to update the depsgraph on rename
436  * and linking into scenes or removal of basis mball. so take care when changing this code.
437  * 
438  * Main idiot thing here is that the system returns find_basis_mball() objects which fail a is_basis_mball() test.
439  *
440  * Not only that but the depsgraph and their areas depend on this behavior!, so making small fixes here isn't worth it.
441  * - Campbell
442  */
443
444
445 /** \brief Test, if Object *ob is basic MetaBall.
446  *
447  * It test last character of Object ID name. If last character
448  * is digit it return 0, else it return 1.
449  */
450 bool BKE_mball_is_basis(Object *ob)
451 {
452         /* just a quick test */
453         const int len = strlen(ob->id.name);
454         return (!isdigit(ob->id.name[len - 1]));
455 }
456
457 /* return nonzero if ob1 is a basis mball for ob */
458 bool BKE_mball_is_basis_for(Object *ob1, Object *ob2)
459 {
460         int basis1nr, basis2nr;
461         char basis1name[MAX_ID_NAME], basis2name[MAX_ID_NAME];
462
463         BLI_split_name_num(basis1name, &basis1nr, ob1->id.name + 2, '.');
464         BLI_split_name_num(basis2name, &basis2nr, ob2->id.name + 2, '.');
465
466         if (!strcmp(basis1name, basis2name)) {
467                 return BKE_mball_is_basis(ob1);
468         }
469         else {
470                 return false;
471         }
472 }
473
474 /* \brief copy some properties from object to other metaball object with same base name
475  *
476  * When some properties (wiresize, threshold, update flags) of metaball are changed, then this properties
477  * are copied to all metaballs in same "group" (metaballs with same base name: MBall,
478  * MBall.001, MBall.002, etc). The most important is to copy properties to the base metaball,
479  * because this metaball influence polygonisation of metaballs. */
480 void BKE_mball_properties_copy(Scene *scene, Object *active_object)
481 {
482         Scene *sce_iter = scene;
483         Base *base;
484         Object *ob;
485         MetaBall *active_mball = (MetaBall *)active_object->data;
486         int basisnr, obnr;
487         char basisname[MAX_ID_NAME], obname[MAX_ID_NAME];
488         SceneBaseIter iter;
489         EvaluationContext *eval_ctx = G.main->eval_ctx;
490
491         BLI_split_name_num(basisname, &basisnr, active_object->id.name + 2, '.');
492
493         /* XXX recursion check, see scene.c, just too simple code this BKE_scene_base_iter_next() */
494         if (F_ERROR == BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 0, NULL, NULL))
495                 return;
496         
497         while (BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 1, &base, &ob)) {
498                 if (ob->type == OB_MBALL) {
499                         if (ob != active_object) {
500                                 BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
501
502                                 /* Object ob has to be in same "group" ... it means, that it has to have
503                                  * same base of its name */
504                                 if (strcmp(obname, basisname) == 0) {
505                                         MetaBall *mb = ob->data;
506
507                                         /* Copy properties from selected/edited metaball */
508                                         mb->wiresize = active_mball->wiresize;
509                                         mb->rendersize = active_mball->rendersize;
510                                         mb->thresh = active_mball->thresh;
511                                         mb->flag = active_mball->flag;
512                                 }
513                         }
514                 }
515         }
516 }
517
518 /** \brief This function finds basic MetaBall.
519  *
520  * Basic MetaBall doesn't include any number at the end of
521  * its name. All MetaBalls with same base of name can be
522  * blended. MetaBalls with different basic name can't be
523  * blended.
524  *
525  * warning!, is_basis_mball() can fail on returned object, see long note above.
526  */
527 Object *BKE_mball_basis_find(Scene *scene, Object *basis)
528 {
529         Scene *sce_iter = scene;
530         Base *base;
531         Object *ob, *bob = basis;
532         int basisnr, obnr;
533         char basisname[MAX_ID_NAME], obname[MAX_ID_NAME];
534         SceneBaseIter iter;
535         EvaluationContext *eval_ctx = G.main->eval_ctx;
536
537         BLI_split_name_num(basisname, &basisnr, basis->id.name + 2, '.');
538
539         /* XXX recursion check, see scene.c, just too simple code this BKE_scene_base_iter_next() */
540         if (F_ERROR == BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 0, NULL, NULL))
541                 return NULL;
542
543         while (BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 1, &base, &ob)) {
544                 if (ob->type == OB_MBALL) {
545                         if (ob != bob) {
546                                 BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
547
548                                 /* object ob has to be in same "group" ... it means, that it has to have
549                                  * same base of its name */
550                                 if (strcmp(obname, basisname) == 0) {
551                                         if (obnr < basisnr) {
552                                                 if (!(ob->flag & OB_FROMDUPLI)) {
553                                                         basis = ob;
554                                                         basisnr = obnr;
555                                                 }
556                                         }
557                                 }
558                         }
559                 }
560         }
561
562         return basis;
563 }
564
565
566 /* ******************** ARITH ************************* */
567
568 /* BASED AT CODE (but mostly rewritten) :
569  * C code from the article
570  * "An Implicit Surface Polygonizer"
571  * by Jules Bloomenthal, jbloom@beauty.gmu.edu
572  * in "Graphics Gems IV", Academic Press, 1994
573  *
574  * Authored by Jules Bloomenthal, Xerox PARC.
575  * Copyright (c) Xerox Corporation, 1991.  All rights reserved.
576  * Permission is granted to reproduce, use and distribute this code for
577  * any and all purposes, provided that this notice appears in all copies. */
578
579 #define RES 12 /* # converge iterations    */
580
581 #define L   0  /* left direction:       -x, -i */
582 #define R   1  /* right direction:      +x, +i */
583 #define B   2  /* bottom direction: -y, -j */
584 #define T   3  /* top direction:        +y, +j */
585 #define N   4  /* near direction:       -z, -k */
586 #define F   5  /* far direction:        +z, +k */
587 #define LBN 0  /* left bottom near corner  */
588 #define LBF 1  /* left bottom far corner   */
589 #define LTN 2  /* left top near corner     */
590 #define LTF 3  /* left top far corner      */
591 #define RBN 4  /* right bottom near corner */
592 #define RBF 5  /* right bottom far corner  */
593 #define RTN 6  /* right top near corner    */
594 #define RTF 7  /* right top far corner     */
595
596 /* the LBN corner of cube (i, j, k), corresponds with location
597  * (i-0.5)*size, (j-0.5)*size, (k-0.5)*size) */
598
599 #define HASHBIT     (5)
600 #define HASHSIZE    (size_t)(1 << (3 * HASHBIT))   /*! < hash table size (32768) */
601
602 #define HASH(i, j, k) ((((( (i) & 31) << 5) | ( (j) & 31)) << 5) | ( (k) & 31) )
603
604 #define MB_BIT(i, bit) (((i) >> (bit)) & 1)
605 #define FLIP(i, bit) ((i) ^ 1 << (bit)) /* flip the given bit of i */
606
607
608 /* **************** POLYGONIZATION ************************ */
609
610 static void calc_mballco(MetaElem *ml, float vec[3])
611 {
612         if (ml->mat) {
613                 mul_m4_v3((float (*)[4])ml->mat, vec);
614         }
615 }
616
617 static float densfunc(MetaElem *ball, float x, float y, float z)
618 {
619         float dist2;
620         float dvec[3] = {x, y, z};
621
622         mul_m4_v3((float (*)[4])ball->imat, dvec);
623
624         switch (ball->type) {
625                 case MB_BALL:
626                         /* do nothing */
627                         break;
628                 case MB_TUBE:
629                         if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
630                         else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
631                         else                            dvec[0] = 0.0;
632                         break;
633                 case MB_PLANE:
634                         if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
635                         else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
636                         else                            dvec[0] = 0.0;
637                         if      (dvec[1] >  ball->expy) dvec[1] -= ball->expy;
638                         else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
639                         else                            dvec[1] = 0.0;
640                         break;
641                 case MB_ELIPSOID:
642                         dvec[0] /= ball->expx;
643                         dvec[1] /= ball->expy;
644                         dvec[2] /= ball->expz;
645                         break;
646                 case MB_CUBE:
647                         if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
648                         else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
649                         else                            dvec[0] = 0.0;
650
651                         if      (dvec[1] >  ball->expy) dvec[1] -= ball->expy;
652                         else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
653                         else                            dvec[1] = 0.0;
654
655                         if      (dvec[2] >  ball->expz) dvec[2] -= ball->expz;
656                         else if (dvec[2] < -ball->expz) dvec[2] += ball->expz;
657                         else                            dvec[2] = 0.0;
658                         break;
659
660                 /* *** deprecated, could be removed?, do-versioned at least *** */
661                 case MB_TUBEX:
662                         if      (dvec[0] >  ball->len) dvec[0] -= ball->len;
663                         else if (dvec[0] < -ball->len) dvec[0] += ball->len;
664                         else                           dvec[0] = 0.0;
665                         break;
666                 case MB_TUBEY:
667                         if      (dvec[1] >  ball->len) dvec[1] -= ball->len;
668                         else if (dvec[1] < -ball->len) dvec[1] += ball->len;
669                         else                           dvec[1] = 0.0;
670                         break;
671                 case MB_TUBEZ:
672                         if      (dvec[2] >  ball->len) dvec[2] -= ball->len;
673                         else if (dvec[2] < -ball->len) dvec[2] += ball->len;
674                         else                           dvec[2] = 0.0;
675                         break;
676                 /* *** end deprecated *** */
677         }
678
679         dist2 = 1.0f - (len_squared_v3(dvec) / ball->rad2);
680
681         if ((ball->flag & MB_NEGATIVE) == 0) {
682                 return (dist2 < 0.0f) ? -0.5f : (ball->s * dist2 * dist2 * dist2) - 0.5f;
683         }
684         else {
685                 return (dist2 < 0.0f) ? 0.5f : 0.5f - (ball->s * dist2 * dist2 * dist2);
686         }
687 }
688
689 static octal_node *find_metaball_octal_node(octal_node *node, float x, float y, float z, short depth)
690 {
691         if (!depth) return node;
692         
693         if (z < node->z) {
694                 if (y < node->y) {
695                         if (x < node->x) {
696                                 if (node->nodes[0])
697                                         return find_metaball_octal_node(node->nodes[0], x, y, z, depth--);
698                                 else
699                                         return node;
700                         }
701                         else {
702                                 if (node->nodes[1])
703                                         return find_metaball_octal_node(node->nodes[1], x, y, z, depth--);
704                                 else
705                                         return node;
706                         }
707                 }
708                 else {
709                         if (x < node->x) {
710                                 if (node->nodes[3])
711                                         return find_metaball_octal_node(node->nodes[3], x, y, z, depth--);
712                                 else
713                                         return node;
714                         }
715                         else {
716                                 if (node->nodes[2])
717                                         return find_metaball_octal_node(node->nodes[2], x, y, z, depth--);
718                                 else
719                                         return node;
720                         }
721                 }
722         }
723         else {
724                 if (y < node->y) {
725                         if (x < node->x) {
726                                 if (node->nodes[4])
727                                         return find_metaball_octal_node(node->nodes[4], x, y, z, depth--);
728                                 else
729                                         return node;
730                         }
731                         else {
732                                 if (node->nodes[5])
733                                         return find_metaball_octal_node(node->nodes[5], x, y, z, depth--);
734                                 else
735                                         return node;
736                         }
737                 }
738                 else {
739                         if (x < node->x) {
740                                 if (node->nodes[7])
741                                         return find_metaball_octal_node(node->nodes[7], x, y, z, depth--);
742                                 else
743                                         return node;
744                         }
745                         else {
746                                 if (node->nodes[6])
747                                         return find_metaball_octal_node(node->nodes[6], x, y, z, depth--);
748                                 else
749                                         return node;
750                         }
751                 }
752         }
753         
754         return node;
755 }
756
757 static float metaball(PROCESS *process, float x, float y, float z)
758 /*  float x, y, z; */
759 {
760         octal_tree *metaball_tree = process->metaball_tree;
761         struct octal_node *node;
762         struct ml_pointer *ml_p;
763         float dens = 0;
764         int a;
765         
766         if (process->totelem > 1) {
767                 node = find_metaball_octal_node(metaball_tree->first, x, y, z, metaball_tree->depth);
768                 if (node) {
769                         for (ml_p = node->elems.first; ml_p; ml_p = ml_p->next) {
770                                 dens += densfunc(ml_p->ml, x, y, z);
771                         }
772
773                         dens += -0.5f * (metaball_tree->pos - node->pos);
774                         dens +=  0.5f * (metaball_tree->neg - node->neg);
775                 }
776                 else {
777                         for (a = 0; a < process->totelem; a++) {
778                                 dens += densfunc(process->mainb[a], x, y, z);
779                         }
780                 }
781         }
782         else {
783                 dens += densfunc(process->mainb[0], x, y, z);
784         }
785
786         return process->thresh - dens;
787 }
788
789 /* ******************************************** */
790
791 static void accum_mballfaces(PROCESS *process, int i1, int i2, int i3, int i4)
792 {
793         int *newi, *cur;
794         /* static int i = 0; I would like to delete altogether, but I don't dare to, yet */
795
796         if (process->totindex == process->curindex) {
797                 process->totindex += 256;
798                 newi = MEM_mallocN(4 * sizeof(int) * process->totindex, "vertindex");
799                 
800                 if (process->indices) {
801                         memcpy(newi, process->indices, 4 * sizeof(int) * (process->totindex - 256));
802                         MEM_freeN(process->indices);
803                 }
804                 process->indices = newi;
805         }
806         
807         cur = process->indices + 4 * process->curindex;
808
809         /* displists now support array drawing, we treat tri's as fake quad */
810         
811         cur[0] = i1;
812         cur[1] = i2;
813         cur[2] = i3;
814         if (i4 == 0)
815                 cur[3] = i3;
816         else 
817                 cur[3] = i4;
818         
819         process->curindex++;
820
821 }
822
823 /* ******************* MEMORY MANAGEMENT *********************** */
824 static void *new_pgn_element(PROCESS *process, int size)
825 {
826         /* during polygonize 1000s of elements are allocated
827          * and never freed in between. Freeing only done at the end.
828          */
829         int blocksize = 16384;
830         void *adr;
831         
832         if (size > 10000 || size == 0) {
833                 printf("incorrect use of new_pgn_element\n");
834         }
835         else if (size == -1) {
836                 struct pgn_elements *cur = process->pgn_list.first;
837                 while (cur) {
838                         MEM_freeN(cur->data);
839                         cur = cur->next;
840                 }
841                 BLI_freelistN(&process->pgn_list);
842                 
843                 return NULL;
844         }
845         
846         size = 4 * ( (size + 3) / 4);
847         
848         if (process->pgn_current) {
849                 if (size + process->pgn_offset < blocksize) {
850                         adr = (void *) (process->pgn_current->data + process->pgn_offset);
851                         process->pgn_offset += size;
852                         return adr;
853                 }
854         }
855         
856         process->pgn_current = MEM_callocN(sizeof(struct pgn_elements), "newpgn");
857         process->pgn_current->data = MEM_callocN(blocksize, "newpgn");
858         BLI_addtail(&process->pgn_list, process->pgn_current);
859         
860         process->pgn_offset = size;
861         return process->pgn_current->data;
862 }
863
864 static void freepolygonize(PROCESS *process)
865 {
866         MEM_freeN(process->corners);
867         MEM_freeN(process->edges);
868         MEM_freeN(process->centers);
869
870         new_pgn_element(process, -1);
871
872         if (process->vertices.ptr) {
873                 MEM_freeN(process->vertices.ptr);
874         }
875 }
876
877 /**** Cubical Polygonization (optional) ****/
878
879 #define LB  0  /* left bottom edge      */
880 #define LT  1  /* left top edge */
881 #define LN  2  /* left near edge        */
882 #define LF  3  /* left far edge */
883 #define RB  4  /* right bottom edge */
884 #define RT  5  /* right top edge        */
885 #define RN  6  /* right near edge       */
886 #define RF  7  /* right far edge        */
887 #define BN  8  /* bottom near edge      */
888 #define BF  9  /* bottom far edge       */
889 #define TN  10 /* top near edge */
890 #define TF  11 /* top far edge  */
891
892 static INTLISTS *cubetable[256];
893
894 /* edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
895 static int corner1[12] = {
896         LBN, LTN, LBN, LBF, RBN, RTN, RBN, RBF, LBN, LBF, LTN, LTF
897 };
898 static int corner2[12] = {
899         LBF, LTF, LTN, LTF, RBF, RTF, RTN, RTF, RBN, RBF, RTN, RTF
900 };
901 static int leftface[12] = {
902         B,  L,  L,  F,  R,  T,  N,  R,  N,  B,  T,  F
903 };
904 /* face on left when going corner1 to corner2 */
905 static int rightface[12] = {
906         L,  T,  N,  L,  B,  R,  R,  F,  B,  F,  N,  T
907 };
908 /* face on right when going corner1 to corner2 */
909
910
911 /* docube: triangulate the cube directly, without decomposition */
912
913 static void docube(PROCESS *process, CUBE *cube, MetaBall *mb)
914 {
915         INTLISTS *polys;
916         CORNER *c1, *c2;
917         int i, index = 0, count, indexar[8];
918         
919         for (i = 0; i < 8; i++) if (cube->corners[i]->value > 0.0f) index += (1 << i);
920         
921         for (polys = cubetable[index]; polys; polys = polys->next) {
922                 INTLIST *edges;
923                 
924                 count = 0;
925                 
926                 for (edges = polys->list; edges; edges = edges->next) {
927                         c1 = cube->corners[corner1[edges->i]];
928                         c2 = cube->corners[corner2[edges->i]];
929                         
930                         indexar[count] = vertid(process, c1, c2, mb);
931                         count++;
932                 }
933                 if (count > 2) {
934                         switch (count) {
935                                 case 3:
936                                         accum_mballfaces(process, indexar[2], indexar[1], indexar[0], 0);
937                                         break;
938                                 case 4:
939                                         if (indexar[0] == 0) accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
940                                         else accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
941                                         break;
942                                 case 5:
943                                         if (indexar[0] == 0) accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
944                                         else accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
945                                 
946                                         accum_mballfaces(process, indexar[4], indexar[3], indexar[0], 0);
947                                         break;
948                                 case 6:
949                                         if (indexar[0] == 0) {
950                                                 accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
951                                                 accum_mballfaces(process, indexar[0], indexar[5], indexar[4], indexar[3]);
952                                         }
953                                         else {
954                                                 accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
955                                                 accum_mballfaces(process, indexar[5], indexar[4], indexar[3], indexar[0]);
956                                         }
957                                         break;
958                                 case 7:
959                                         if (indexar[0] == 0) {
960                                                 accum_mballfaces(process, indexar[0], indexar[3], indexar[2], indexar[1]);
961                                                 accum_mballfaces(process, indexar[0], indexar[5], indexar[4], indexar[3]);
962                                         }
963                                         else {
964                                                 accum_mballfaces(process, indexar[3], indexar[2], indexar[1], indexar[0]);
965                                                 accum_mballfaces(process, indexar[5], indexar[4], indexar[3], indexar[0]);
966                                         }
967                                 
968                                         accum_mballfaces(process, indexar[6], indexar[5], indexar[0], 0);
969
970                                         break;
971                         }
972                 }
973         }
974 }
975
976
977 /* testface: given cube at lattice (i, j, k), and four corners of face,
978  * if surface crosses face, compute other four corners of adjacent cube
979  * and add new cube to cube stack */
980
981 static void testface(PROCESS *process, int i, int j, int k, CUBE *old, int bit, int c1, int c2, int c3, int c4)
982 {
983         CUBE newc;
984         CUBES *oldcubes = process->cubes;
985         CORNER *corn1, *corn2, *corn3, *corn4;
986         int n, pos;
987
988         corn1 = old->corners[c1];
989         corn2 = old->corners[c2];
990         corn3 = old->corners[c3];
991         corn4 = old->corners[c4];
992         
993         pos = corn1->value > 0.0f ? 1 : 0;
994
995         /* test if no surface crossing */
996         if ( (corn2->value > 0) == pos && (corn3->value > 0) == pos && (corn4->value > 0) == pos) return;
997         /* test if cube out of bounds */
998         /*if ( abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) return;*/
999         /* test if already visited (always as last) */
1000         if (setcenter(process, process->centers, i, j, k)) {
1001                 return;
1002         }
1003
1004         /* create new cube and add cube to top of stack: */
1005         process->cubes = (CUBES *) new_pgn_element(process, sizeof(CUBES));
1006         process->cubes->next = oldcubes;
1007         
1008         newc.i = i;
1009         newc.j = j;
1010         newc.k = k;
1011         for (n = 0; n < 8; n++) newc.corners[n] = NULL;
1012         
1013         newc.corners[FLIP(c1, bit)] = corn1;
1014         newc.corners[FLIP(c2, bit)] = corn2;
1015         newc.corners[FLIP(c3, bit)] = corn3;
1016         newc.corners[FLIP(c4, bit)] = corn4;
1017
1018         if (newc.corners[0] == NULL) newc.corners[0] = setcorner(process, i, j, k);
1019         if (newc.corners[1] == NULL) newc.corners[1] = setcorner(process, i, j, k + 1);
1020         if (newc.corners[2] == NULL) newc.corners[2] = setcorner(process, i, j + 1, k);
1021         if (newc.corners[3] == NULL) newc.corners[3] = setcorner(process, i, j + 1, k + 1);
1022         if (newc.corners[4] == NULL) newc.corners[4] = setcorner(process, i + 1, j, k);
1023         if (newc.corners[5] == NULL) newc.corners[5] = setcorner(process, i + 1, j, k + 1);
1024         if (newc.corners[6] == NULL) newc.corners[6] = setcorner(process, i + 1, j + 1, k);
1025         if (newc.corners[7] == NULL) newc.corners[7] = setcorner(process, i + 1, j + 1, k + 1);
1026
1027         process->cubes->cube = newc;
1028 }
1029
1030 /* setcorner: return corner with the given lattice location
1031  * set (and cache) its function value */
1032
1033 static CORNER *setcorner(PROCESS *process, int i, int j, int k)
1034 {
1035         /* for speed, do corner value caching here */
1036         CORNER *c;
1037         int index;
1038
1039         /* does corner exist? */
1040         index = HASH(i, j, k);
1041         c = process->corners[index];
1042         
1043         for (; c != NULL; c = c->next) {
1044                 if (c->i == i && c->j == j && c->k == k) {
1045                         return c;
1046                 }
1047         }
1048
1049         c = (CORNER *) new_pgn_element(process, sizeof(CORNER));
1050
1051         c->i = i; 
1052         c->co[0] = ((float)i - 0.5f) * process->size;
1053         c->j = j; 
1054         c->co[1] = ((float)j - 0.5f) * process->size;
1055         c->k = k; 
1056         c->co[2] = ((float)k - 0.5f) * process->size;
1057         c->value = process->function(process, c->co[0], c->co[1], c->co[2]);
1058         
1059         c->next = process->corners[index];
1060         process->corners[index] = c;
1061         
1062         return c;
1063 }
1064
1065
1066 /* nextcwedge: return next clockwise edge from given edge around given face */
1067
1068 static int nextcwedge(int edge, int face)
1069 {
1070         switch (edge) {
1071                 case LB:
1072                         return (face == L) ? LF : BN;
1073                 case LT:
1074                         return (face == L) ? LN : TF;
1075                 case LN:
1076                         return (face == L) ? LB : TN;
1077                 case LF:
1078                         return (face == L) ? LT : BF;
1079                 case RB:
1080                         return (face == R) ? RN : BF;
1081                 case RT:
1082                         return (face == R) ? RF : TN;
1083                 case RN:
1084                         return (face == R) ? RT : BN;
1085                 case RF:
1086                         return (face == R) ? RB : TF;
1087                 case BN:
1088                         return (face == B) ? RB : LN;
1089                 case BF:
1090                         return (face == B) ? LB : RF;
1091                 case TN:
1092                         return (face == T) ? LT : RN;
1093                 case TF:
1094                         return (face == T) ? RT : LF;
1095         }
1096         return 0;
1097 }
1098
1099
1100 /* otherface: return face adjoining edge that is not the given face */
1101
1102 static int otherface(int edge, int face)
1103 {
1104         int other = leftface[edge];
1105         return face == other ? rightface[edge] : other;
1106 }
1107
1108
1109 /* makecubetable: create the 256 entry table for cubical polygonization */
1110
1111 static void makecubetable(void)
1112 {
1113         static int is_done = FALSE;
1114         int i, e, c, done[12], pos[8];
1115
1116         if (is_done) return;
1117         is_done = TRUE;
1118
1119         for (i = 0; i < 256; i++) {
1120                 for (e = 0; e < 12; e++) done[e] = 0;
1121                 for (c = 0; c < 8; c++) pos[c] = MB_BIT(i, c);
1122                 for (e = 0; e < 12; e++)
1123                         if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
1124                                 INTLIST *ints = NULL;
1125                                 INTLISTS *lists = (INTLISTS *) MEM_callocN(sizeof(INTLISTS), "mball_intlist");
1126                                 int start = e, edge = e;
1127                                 
1128                                 /* get face that is to right of edge from pos to neg corner: */
1129                                 int face = pos[corner1[e]] ? rightface[e] : leftface[e];
1130                                 
1131                                 while (1) {
1132                                         edge = nextcwedge(edge, face);
1133                                         done[edge] = 1;
1134                                         if (pos[corner1[edge]] != pos[corner2[edge]]) {
1135                                                 INTLIST *tmp = ints;
1136                                                 
1137                                                 ints = (INTLIST *) MEM_callocN(sizeof(INTLIST), "mball_intlist");
1138                                                 ints->i = edge;
1139                                                 ints->next = tmp; /* add edge to head of list */
1140                                                 
1141                                                 if (edge == start) break;
1142                                                 face = otherface(edge, face);
1143                                         }
1144                                 }
1145                                 lists->list = ints; /* add ints to head of table entry */
1146                                 lists->next = cubetable[i];
1147                                 cubetable[i] = lists;
1148                         }
1149         }
1150 }
1151
1152 void BKE_mball_cubeTable_free(void)
1153 {
1154         int i;
1155         INTLISTS *lists, *nlists;
1156         INTLIST *ints, *nints;
1157
1158         for (i = 0; i < 256; i++) {
1159                 lists = cubetable[i];
1160                 while (lists) {
1161                         nlists = lists->next;
1162                         
1163                         ints = lists->list;
1164                         while (ints) {
1165                                 nints = ints->next;
1166                                 MEM_freeN(ints);
1167                                 ints = nints;
1168                         }
1169                         
1170                         MEM_freeN(lists);
1171                         lists = nlists;
1172                 }
1173                 cubetable[i] = NULL;
1174         }
1175 }
1176
1177 /**** Storage ****/
1178
1179 /* setcenter: set (i, j, k) entry of table[]
1180  * return 1 if already set; otherwise, set and return 0 */
1181
1182 static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const int j, const int k)
1183 {
1184         int index;
1185         CENTERLIST *newc, *l, *q;
1186
1187         index = HASH(i, j, k);
1188         q = table[index];
1189
1190         for (l = q; l != NULL; l = l->next) {
1191                 if (l->i == i && l->j == j && l->k == k) return 1;
1192         }
1193         
1194         newc = (CENTERLIST *) new_pgn_element(process, sizeof(CENTERLIST));
1195         newc->i = i; 
1196         newc->j = j; 
1197         newc->k = k; 
1198         newc->next = q;
1199         table[index] = newc;
1200         
1201         return 0;
1202 }
1203
1204
1205 /* setedge: set vertex id for edge */
1206
1207 static void setedge(PROCESS *process,
1208                     EDGELIST *table[],
1209                     int i1, int j1,
1210                     int k1, int i2,
1211                     int j2, int k2,
1212                     int vid)
1213 {
1214         unsigned int index;
1215         EDGELIST *newe;
1216         
1217         if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
1218                 int t = i1;
1219                 i1 = i2;
1220                 i2 = t;
1221                 t = j1;
1222                 j1 = j2;
1223                 j2 = t;
1224                 t = k1;
1225                 k1 = k2;
1226                 k2 = t;
1227         }
1228         index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
1229         newe = (EDGELIST *) new_pgn_element(process, sizeof(EDGELIST));
1230         newe->i1 = i1; 
1231         newe->j1 = j1; 
1232         newe->k1 = k1;
1233         newe->i2 = i2; 
1234         newe->j2 = j2; 
1235         newe->k2 = k2;
1236         newe->vid = vid;
1237         newe->next = table[index];
1238         table[index] = newe;
1239 }
1240
1241
1242 /* getedge: return vertex id for edge; return -1 if not set */
1243
1244 static int getedge(EDGELIST *table[],
1245                    int i1, int j1, int k1,
1246                    int i2, int j2, int k2)
1247 {
1248         EDGELIST *q;
1249         
1250         if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
1251                 int t = i1;
1252                 i1 = i2;
1253                 i2 = t;
1254                 t = j1;
1255                 j1 = j2;
1256                 j2 = t;
1257                 t = k1;
1258                 k1 = k2;
1259                 k2 = t;
1260         }
1261         q = table[HASH(i1, j1, k1) + HASH(i2, j2, k2)];
1262         for (; q != NULL; q = q->next) {
1263                 if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 &&
1264                     q->i2 == i2 && q->j2 == j2 && q->k2 == k2)
1265                 {
1266                         return q->vid;
1267                 }
1268         }
1269         return -1;
1270 }
1271
1272
1273 /**** Vertices ****/
1274
1275 #undef R
1276
1277
1278
1279 /* vertid: return index for vertex on edge:
1280  * c1->value and c2->value are presumed of different sign
1281  * return saved index if any; else compute vertex and save */
1282
1283 /* addtovertices: add v to sequence of vertices */
1284
1285 static void addtovertices(VERTICES *vertices, VERTEX v)
1286 {
1287         if (vertices->count == vertices->max) {
1288                 int i;
1289                 VERTEX *newv;
1290                 vertices->max = vertices->count == 0 ? 10 : 2 * vertices->count;
1291                 newv = (VERTEX *) MEM_callocN(vertices->max * sizeof(VERTEX), "addtovertices");
1292                 
1293                 for (i = 0; i < vertices->count; i++) newv[i] = vertices->ptr[i];
1294                 
1295                 if (vertices->ptr != NULL) MEM_freeN(vertices->ptr);
1296                 vertices->ptr = newv;
1297         }
1298         vertices->ptr[vertices->count++] = v;
1299 }
1300
1301 /* vnormal: compute unit length surface normal at point */
1302
1303 static void vnormal(PROCESS *process, const float point[3], float r_no[3])
1304 {
1305         const float delta = 0.2f * process->delta;
1306         const float f = process->function(process, point[0], point[1], point[2]);
1307
1308         r_no[0] = process->function(process, point[0] + delta, point[1], point[2]) - f;
1309         r_no[1] = process->function(process, point[0], point[1] + delta, point[2]) - f;
1310         r_no[2] = process->function(process, point[0], point[1], point[2] + delta) - f;
1311
1312 #if 1
1313         normalize_v3(r_no);
1314 #else
1315         f = normalize_v3(r_no);
1316         
1317         if (0) {
1318                 float tvec[3];
1319                 
1320                 delta *= 2.0f;
1321                 
1322                 f = process->function(process, point[0], point[1], point[2]);
1323         
1324                 tvec[0] = process->function(process, point[0] + delta, point[1], point[2]) - f;
1325                 tvec[1] = process->function(process, point[0], point[1] + delta, point[2]) - f;
1326                 tvec[2] = process->function(process, point[0], point[1], point[2] + delta) - f;
1327         
1328                 if (normalize_v3(tvec) != 0.0f) {
1329                         add_v3_v3(r_no, tvec);
1330                         normalize_v3(r_no);
1331                 }
1332         }
1333 #endif
1334 }
1335
1336
1337 static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2, MetaBall *mb)
1338 {
1339         VERTEX v;
1340         int vid = getedge(process->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
1341
1342         if (vid != -1) {
1343                 return vid;  /* previously computed */
1344         }
1345
1346         converge(process, c1->co, c2->co, c1->value, c2->value, v.co, mb, 1); /* position */
1347         vnormal(process, v.co, v.no);
1348
1349         addtovertices(&process->vertices, v);            /* save vertex */
1350         vid = process->vertices.count - 1;
1351         setedge(process, process->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
1352         
1353         return vid;
1354 }
1355
1356
1357 /* converge: from two points of differing sign, converge to zero crossing */
1358 /* watch it: p1 and p2 are used to calculate */
1359 static void converge(PROCESS *process, const float p1[3], const float p2[3], float v1, float v2,
1360                      float p[3], MetaBall *mb, int f)
1361 {
1362         int i = 0;
1363         float pos[3], neg[3];
1364         float positive = 0.0f, negative = 0.0f;
1365         float dvec[3];
1366         
1367         if (v1 < 0) {
1368                 copy_v3_v3(pos, p2);
1369                 copy_v3_v3(neg, p1);
1370                 positive = v2;
1371                 negative = v1;
1372         }
1373         else {
1374                 copy_v3_v3(pos, p1);
1375                 copy_v3_v3(neg, p2);
1376                 positive = v1;
1377                 negative = v2;
1378         }
1379
1380         sub_v3_v3v3(dvec, pos, neg);
1381
1382 /* Approximation by linear interpolation is faster then binary subdivision,
1383  * but it results sometimes (mb->thresh < 0.2) into the strange results */
1384         if ((mb->thresh > 0.2f) && (f == 1)) {
1385                 if ((dvec[1] == 0.0f) && (dvec[2] == 0.0f)) {
1386                         p[0] = neg[0] - negative * dvec[0] / (positive - negative);
1387                         p[1] = neg[1];
1388                         p[2] = neg[2];
1389                         return;
1390                 }
1391                 if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
1392                         p[0] = neg[0];
1393                         p[1] = neg[1] - negative * dvec[1] / (positive - negative);
1394                         p[2] = neg[2];
1395                         return;
1396                 }
1397                 if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
1398                         p[0] = neg[0];
1399                         p[1] = neg[1];
1400                         p[2] = neg[2] - negative * dvec[2] / (positive - negative);
1401                         return;
1402                 }
1403         }
1404
1405         if ((dvec[1] == 0.0f) && (dvec[2] == 0.0f)) {
1406                 p[1] = neg[1];
1407                 p[2] = neg[2];
1408                 while (1) {
1409                         if (i++ == RES) return;
1410                         p[0] = 0.5f * (pos[0] + neg[0]);
1411                         if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[0] = p[0];
1412                         else                                                       neg[0] = p[0];
1413                 }
1414         }
1415
1416         if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
1417                 p[0] = neg[0];
1418                 p[2] = neg[2];
1419                 while (1) {
1420                         if (i++ == RES) return;
1421                         p[1] = 0.5f * (pos[1] + neg[1]);
1422                         if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[1] = p[1];
1423                         else                                                       neg[1] = p[1];
1424                 }
1425         }
1426
1427         if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
1428                 p[0] = neg[0];
1429                 p[1] = neg[1];
1430                 while (1) {
1431                         if (i++ == RES) return;
1432                         p[2] = 0.5f * (pos[2] + neg[2]);
1433                         if ((process->function(process, p[0], p[1], p[2])) > 0.0f) pos[2] = p[2];
1434                         else                                                       neg[2] = p[2];
1435                 }
1436         }
1437
1438         /* This is necessary to find start point */
1439         while (1) {
1440                 mid_v3_v3v3(&p[0], pos, neg);
1441
1442                 if (i++ == RES) {
1443                         return;
1444                 }
1445
1446                 if ((process->function(process, p[0], p[1], p[2])) > 0.0f) {
1447                         copy_v3_v3(pos, &p[0]);
1448                 }
1449                 else {
1450                         copy_v3_v3(neg, &p[0]);
1451                 }
1452         }
1453 }
1454
1455 /* ************************************** */
1456 static void add_cube(PROCESS *process, int i, int j, int k, int count)
1457 {
1458         CUBES *ncube;
1459         int n;
1460         int a, b, c;
1461
1462         /* hmmm, not only one, but eight cube will be added on the stack 
1463          * ... */
1464         for (a = i - 1; a < i + count; a++)
1465                 for (b = j - 1; b < j + count; b++)
1466                         for (c = k - 1; c < k + count; c++) {
1467                                 /* test if cube has been found before */
1468                                 if (setcenter(process, process->centers, a, b, c) == 0) {
1469                                         /* push cube on stack: */
1470                                         ncube = (CUBES *) new_pgn_element(process, sizeof(CUBES));
1471                                         ncube->next = process->cubes;
1472                                         process->cubes = ncube;
1473
1474                                         ncube->cube.i = a;
1475                                         ncube->cube.j = b;
1476                                         ncube->cube.k = c;
1477
1478                                         /* set corners of initial cube: */
1479                                         for (n = 0; n < 8; n++)
1480                                                 ncube->cube.corners[n] = setcorner(process, a + MB_BIT(n, 2), b + MB_BIT(n, 1), c + MB_BIT(n, 0));
1481                                 }
1482                         }
1483 }
1484
1485
1486 static void find_first_points(PROCESS *process, MetaBall *mb, int a)
1487 {
1488         MetaElem *ml;
1489         float f;
1490
1491         ml = process->mainb[a];
1492         f = 1.0f - (mb->thresh / ml->s);
1493
1494         /* Skip, when Stiffness of MetaElement is too small ... MetaElement can't be
1495          * visible alone ... but still can influence others MetaElements :-) */
1496         if (f > 0.0f) {
1497                 float IN[3] = {0.0f}, OUT[3] = {0.0f}, in[3] = {0.0f}, out[3];
1498                 int i, j, k, c_i, c_j, c_k;
1499                 int index[3] = {1, 0, -1};
1500                 float in_v /*, out_v*/;
1501                 float workp[3];
1502                 float dvec[3];
1503                 float tmp_v, workp_v, max_len, nx, ny, nz, max_dim;
1504
1505                 calc_mballco(ml, in);
1506                 in_v = process->function(process, in[0], in[1], in[2]);
1507
1508                 for (i = 0; i < 3; i++) {
1509                         switch (ml->type) {
1510                                 case MB_BALL:
1511                                         OUT[0] = out[0] = IN[0] + index[i] * ml->rad;
1512                                         break;
1513                                 case MB_TUBE:
1514                                 case MB_PLANE:
1515                                 case MB_ELIPSOID:
1516                                 case MB_CUBE:
1517                                         OUT[0] = out[0] = IN[0] + index[i] * (ml->expx + ml->rad);
1518                                         break;
1519                         }
1520
1521                         for (j = 0; j < 3; j++) {
1522                                 switch (ml->type) {
1523                                         case MB_BALL:
1524                                                 OUT[1] = out[1] = IN[1] + index[j] * ml->rad;
1525                                                 break;
1526                                         case MB_TUBE:
1527                                         case MB_PLANE:
1528                                         case MB_ELIPSOID:
1529                                         case MB_CUBE:
1530                                                 OUT[1] = out[1] = IN[1] + index[j] * (ml->expy + ml->rad);
1531                                                 break;
1532                                 }
1533                         
1534                                 for (k = 0; k < 3; k++) {
1535                                         out[0] = OUT[0];
1536                                         out[1] = OUT[1];
1537                                         switch (ml->type) {
1538                                                 case MB_BALL:
1539                                                 case MB_TUBE:
1540                                                 case MB_PLANE:
1541                                                         out[2] = IN[2] + index[k] * ml->rad;
1542                                                         break;
1543                                                 case MB_ELIPSOID:
1544                                                 case MB_CUBE:
1545                                                         out[2] = IN[2] + index[k] * (ml->expz + ml->rad);
1546                                                         break;
1547                                         }
1548
1549                                         calc_mballco(ml, out);
1550
1551                                         /*out_v = process->function(out[0], out[1], out[2]);*/ /*UNUSED*/
1552
1553                                         /* find "first points" on Implicit Surface of MetaElemnt ml */
1554                                         copy_v3_v3(workp, in);
1555                                         workp_v = in_v;
1556                                         max_len = len_v3v3(out, in);
1557
1558                                         nx = abs((out[0] - in[0]) / process->size);
1559                                         ny = abs((out[1] - in[1]) / process->size);
1560                                         nz = abs((out[2] - in[2]) / process->size);
1561                                         
1562                                         max_dim = max_fff(nx, ny, nz);
1563                                         if (max_dim != 0.0f) {
1564                                                 float len = 0.0f;
1565
1566                                                 dvec[0] = (out[0] - in[0]) / max_dim;
1567                                                 dvec[1] = (out[1] - in[1]) / max_dim;
1568                                                 dvec[2] = (out[2] - in[2]) / max_dim;
1569
1570                                                 while (len <= max_len) {
1571                                                         add_v3_v3(workp, dvec);
1572
1573                                                         /* compute value of implicite function */
1574                                                         tmp_v = process->function(process, workp[0], workp[1], workp[2]);
1575                                                         /* add cube to the stack, when value of implicite function crosses zero value */
1576                                                         if ((tmp_v < 0.0f && workp_v >= 0.0f) || (tmp_v > 0.0f && workp_v <= 0.0f)) {
1577
1578                                                                 /* indexes of CUBE, which includes "first point" */
1579                                                                 c_i = (int)floor(workp[0] / process->size);
1580                                                                 c_j = (int)floor(workp[1] / process->size);
1581                                                                 c_k = (int)floor(workp[2] / process->size);
1582                                                                 
1583                                                                 /* add CUBE (with indexes c_i, c_j, c_k) to the stack,
1584                                                                  * this cube includes found point of Implicit Surface */
1585                                                                 if ((ml->flag & MB_NEGATIVE) == 0) {
1586                                                                         add_cube(process, c_i, c_j, c_k, 1);
1587                                                                 }
1588                                                                 else {
1589                                                                         add_cube(process, c_i, c_j, c_k, 2);
1590                                                                 }
1591                                                         }
1592                                                         len = len_v3v3(workp, in);
1593                                                         workp_v = tmp_v;
1594
1595                                                 }
1596                                         }
1597                                 }
1598                         }
1599                 }
1600         }
1601 }
1602
1603 static void polygonize(PROCESS *process, MetaBall *mb)
1604 {
1605         CUBE c;
1606         int a;
1607
1608         process->vertices.count = process->vertices.max = 0;
1609         process->vertices.ptr = NULL;
1610
1611         /* allocate hash tables and build cube polygon table: */
1612         process->centers = MEM_callocN(HASHSIZE * sizeof(CENTERLIST *), "mbproc->centers");
1613         process->corners = MEM_callocN(HASHSIZE * sizeof(CORNER *), "mbproc->corners");
1614         process->edges = MEM_callocN(2 * HASHSIZE * sizeof(EDGELIST *), "mbproc->edges");
1615         makecubetable();
1616
1617         for (a = 0; a < process->totelem; a++) {
1618
1619                 /* try to find 8 points on the surface for each MetaElem */
1620                 find_first_points(process, mb, a);
1621         }
1622
1623         /* polygonize all MetaElems of current MetaBall */
1624         while (process->cubes != NULL) { /* process active cubes till none left */
1625                 c = process->cubes->cube;
1626
1627                 /* polygonize the cube directly: */
1628                 docube(process, &c, mb);
1629                 
1630                 /* pop current cube from stack */
1631                 process->cubes = process->cubes->next;
1632                 
1633                 /* test six face directions, maybe add to stack: */
1634                 testface(process, c.i - 1, c.j, c.k, &c, 2, LBN, LBF, LTN, LTF);
1635                 testface(process, c.i + 1, c.j, c.k, &c, 2, RBN, RBF, RTN, RTF);
1636                 testface(process, c.i, c.j - 1, c.k, &c, 1, LBN, LBF, RBN, RBF);
1637                 testface(process, c.i, c.j + 1, c.k, &c, 1, LTN, LTF, RTN, RTF);
1638                 testface(process, c.i, c.j, c.k - 1, &c, 0, LBN, LTN, RBN, RTN);
1639                 testface(process, c.i, c.j, c.k + 1, &c, 0, LBF, LTF, RBF, RTF);
1640         }
1641 }
1642
1643 static float init_meta(EvaluationContext *eval_ctx, PROCESS *process, Scene *scene, Object *ob)    /* return totsize */
1644 {
1645         Scene *sce_iter = scene;
1646         Base *base;
1647         Object *bob;
1648         MetaBall *mb;
1649         MetaElem *ml;
1650         float size, totsize, obinv[4][4], obmat[4][4], vec[3];
1651         //float max = 0.0f;
1652         int a, obnr, zero_size = 0;
1653         char obname[MAX_ID_NAME];
1654         SceneBaseIter iter;
1655
1656         copy_m4_m4(obmat, ob->obmat);   /* to cope with duplicators from BKE_scene_base_iter_next */
1657         invert_m4_m4(obinv, ob->obmat);
1658         a = 0;
1659         
1660         BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
1661         
1662         /* make main array */
1663         BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 0, NULL, NULL);
1664         while (BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 1, &base, &bob)) {
1665
1666                 if (bob->type == OB_MBALL) {
1667                         zero_size = 0;
1668                         ml = NULL;
1669
1670                         if (bob == ob && (base->flag & OB_FROMDUPLI) == 0) {
1671                                 mb = ob->data;
1672         
1673                                 if (mb->editelems) ml = mb->editelems->first;
1674                                 else ml = mb->elems.first;
1675                         }
1676                         else {
1677                                 char name[MAX_ID_NAME];
1678                                 int nr;
1679                                 
1680                                 BLI_split_name_num(name, &nr, bob->id.name + 2, '.');
1681                                 if (strcmp(obname, name) == 0) {
1682                                         mb = bob->data;
1683                                         
1684                                         if (mb->editelems) ml = mb->editelems->first;
1685                                         else ml = mb->elems.first;
1686                                 }
1687                         }
1688
1689                         /* when metaball object has zero scale, then MetaElem to this MetaBall
1690                          * will not be put to mainb array */
1691                         if (has_zero_axis_m4(bob->obmat)) {
1692                                 zero_size = 1;
1693                         }
1694                         else if (bob->parent) {
1695                                 struct Object *pob = bob->parent;
1696                                 while (pob) {
1697                                         if (has_zero_axis_m4(pob->obmat)) {
1698                                                 zero_size = 1;
1699                                                 break;
1700                                         }
1701                                         pob = pob->parent;
1702                                 }
1703                         }
1704
1705                         if (zero_size) {
1706                                 unsigned int ml_count = 0;
1707                                 while (ml) {
1708                                         ml_count++;
1709                                         ml = ml->next;
1710                                 }
1711                                 process->totelem -= ml_count;
1712                         }
1713                         else {
1714                                 while (ml) {
1715                                         if (!(ml->flag & MB_HIDE)) {
1716                                                 int i;
1717                                                 float temp1[4][4], temp2[4][4], temp3[4][4];
1718                                                 float (*mat)[4] = NULL, (*imat)[4] = NULL;
1719                                                 float max_x, max_y, max_z, min_x, min_y, min_z;
1720                                                 float expx, expy, expz;
1721
1722                                                 max_x = max_y = max_z = -3.4e38;
1723                                                 min_x = min_y = min_z =  3.4e38;
1724
1725                                                 /* too big stiffness seems only ugly due to linear interpolation
1726                                                  * no need to have possibility for too big stiffness */
1727                                                 if (ml->s > 10.0f) ml->s = 10.0f;
1728
1729                                                 /* Rotation of MetaElem is stored in quat */
1730                                                 quat_to_mat4(temp3, ml->quat);
1731
1732                                                 /* Translation of MetaElem */
1733                                                 unit_m4(temp2);
1734                                                 temp2[3][0] = ml->x;
1735                                                 temp2[3][1] = ml->y;
1736                                                 temp2[3][2] = ml->z;
1737
1738                                                 mul_m4_m4m4(temp1, temp2, temp3);
1739
1740                                                 /* make a copy because of duplicates */
1741                                                 process->mainb[a] = new_pgn_element(process, sizeof(MetaElem));
1742                                                 *(process->mainb[a]) = *ml;
1743                                                 process->mainb[a]->bb = new_pgn_element(process, sizeof(BoundBox));
1744
1745                                                 mat = new_pgn_element(process, 4 * 4 * sizeof(float));
1746                                                 imat = new_pgn_element(process, 4 * 4 * sizeof(float));
1747
1748                                                 /* mat is the matrix to transform from mball into the basis-mball */
1749                                                 invert_m4_m4(obinv, obmat);
1750                                                 mul_m4_m4m4(temp2, obinv, bob->obmat);
1751                                                 /* MetaBall transformation */
1752                                                 mul_m4_m4m4(mat, temp2, temp1);
1753
1754                                                 invert_m4_m4(imat, mat);
1755
1756                                                 process->mainb[a]->rad2 = ml->rad * ml->rad;
1757
1758                                                 process->mainb[a]->mat = (float *) mat;
1759                                                 process->mainb[a]->imat = (float *) imat;
1760
1761                                                 if (!MB_TYPE_SIZE_SQUARED(ml->type)) {
1762                                                         expx = ml->expx;
1763                                                         expy = ml->expy;
1764                                                         expz = ml->expz;
1765                                                 }
1766                                                 else {
1767                                                         expx = ml->expx * ml->expx;
1768                                                         expy = ml->expy * ml->expy;
1769                                                         expz = ml->expz * ml->expz;
1770                                                 }
1771
1772                                                 /* untransformed Bounding Box of MetaElem */
1773                                                 /* TODO, its possible the elem type has been changed and the exp* values can use a fallback */
1774                                                 copy_v3_fl3(process->mainb[a]->bb->vec[0], -expx, -expy, -expz);  /* 0 */
1775                                                 copy_v3_fl3(process->mainb[a]->bb->vec[1], +expx, -expy, -expz);  /* 1 */
1776                                                 copy_v3_fl3(process->mainb[a]->bb->vec[2], +expx, +expy, -expz);  /* 2 */
1777                                                 copy_v3_fl3(process->mainb[a]->bb->vec[3], -expx, +expy, -expz);  /* 3 */
1778                                                 copy_v3_fl3(process->mainb[a]->bb->vec[4], -expx, -expy, +expz);  /* 4 */
1779                                                 copy_v3_fl3(process->mainb[a]->bb->vec[5], +expx, -expy, +expz);  /* 5 */
1780                                                 copy_v3_fl3(process->mainb[a]->bb->vec[6], +expx, +expy, +expz);  /* 6 */
1781                                                 copy_v3_fl3(process->mainb[a]->bb->vec[7], -expx, +expy, +expz);  /* 7 */
1782
1783                                                 /* transformation of Metalem bb */
1784                                                 for (i = 0; i < 8; i++)
1785                                                         mul_m4_v3((float (*)[4])mat, process->mainb[a]->bb->vec[i]);
1786
1787                                                 /* find max and min of transformed bb */
1788                                                 for (i = 0; i < 8; i++) {
1789                                                         /* find maximums */
1790                                                         if (process->mainb[a]->bb->vec[i][0] > max_x) max_x = process->mainb[a]->bb->vec[i][0];
1791                                                         if (process->mainb[a]->bb->vec[i][1] > max_y) max_y = process->mainb[a]->bb->vec[i][1];
1792                                                         if (process->mainb[a]->bb->vec[i][2] > max_z) max_z = process->mainb[a]->bb->vec[i][2];
1793                                                         /* find  minimums */
1794                                                         if (process->mainb[a]->bb->vec[i][0] < min_x) min_x = process->mainb[a]->bb->vec[i][0];
1795                                                         if (process->mainb[a]->bb->vec[i][1] < min_y) min_y = process->mainb[a]->bb->vec[i][1];
1796                                                         if (process->mainb[a]->bb->vec[i][2] < min_z) min_z = process->mainb[a]->bb->vec[i][2];
1797                                                 }
1798                                         
1799                                                 /* create "new" bb, only point 0 and 6, which are
1800                                                  * necessary for octal tree filling */
1801                                                 process->mainb[a]->bb->vec[0][0] = min_x - ml->rad;
1802                                                 process->mainb[a]->bb->vec[0][1] = min_y - ml->rad;
1803                                                 process->mainb[a]->bb->vec[0][2] = min_z - ml->rad;
1804
1805                                                 process->mainb[a]->bb->vec[6][0] = max_x + ml->rad;
1806                                                 process->mainb[a]->bb->vec[6][1] = max_y + ml->rad;
1807                                                 process->mainb[a]->bb->vec[6][2] = max_z + ml->rad;
1808
1809                                                 a++;
1810                                         }
1811                                         ml = ml->next;
1812                                 }
1813                         }
1814                 }
1815         }
1816
1817         
1818         /* totsize (= 'manhattan' radius) */
1819         totsize = 0.0;
1820         for (a = 0; a < process->totelem; a++) {
1821                 
1822                 vec[0] = process->mainb[a]->x + process->mainb[a]->rad + process->mainb[a]->expx;
1823                 vec[1] = process->mainb[a]->y + process->mainb[a]->rad + process->mainb[a]->expy;
1824                 vec[2] = process->mainb[a]->z + process->mainb[a]->rad + process->mainb[a]->expz;
1825
1826                 calc_mballco(process->mainb[a], vec);
1827         
1828                 size = fabsf(vec[0]);
1829                 if (size > totsize) totsize = size;
1830                 size = fabsf(vec[1]);
1831                 if (size > totsize) totsize = size;
1832                 size = fabsf(vec[2]);
1833                 if (size > totsize) totsize = size;
1834
1835                 vec[0] = process->mainb[a]->x - process->mainb[a]->rad;
1836                 vec[1] = process->mainb[a]->y - process->mainb[a]->rad;
1837                 vec[2] = process->mainb[a]->z - process->mainb[a]->rad;
1838                                 
1839                 calc_mballco(process->mainb[a], vec);
1840         
1841                 size = fabsf(vec[0]);
1842                 if (size > totsize) totsize = size;
1843                 size = fabsf(vec[1]);
1844                 if (size > totsize) totsize = size;
1845                 size = fabsf(vec[2]);
1846                 if (size > totsize) totsize = size;
1847         }
1848
1849         for (a = 0; a < process->totelem; a++) {
1850                 process->thresh += densfunc(process->mainb[a], 2.0f * totsize, 2.0f * totsize, 2.0f * totsize);
1851         }
1852
1853         return totsize;
1854 }
1855
1856 /* if MetaElem lies in node, then node includes MetaElem pointer (ml_p)
1857  * pointing at MetaElem (ml)
1858  */
1859 static void fill_metaball_octal_node(octal_node *node, MetaElem *ml, short i)
1860 {
1861         ml_pointer *ml_p;
1862
1863         ml_p = MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
1864         ml_p->ml = ml;
1865         BLI_addtail(&(node->nodes[i]->elems), ml_p);
1866         node->count++;
1867         
1868         if ((ml->flag & MB_NEGATIVE) == 0) {
1869                 node->nodes[i]->pos++;
1870         }
1871         else {
1872                 node->nodes[i]->neg++;
1873         }
1874 }
1875
1876 /* Node is subdivided as is illustrated on the following figure:
1877  * 
1878  *      +------+------+
1879  *     /      /      /|
1880  *    +------+------+ |
1881  *   /      /      /| +
1882  *  +------+------+ |/|
1883  *  |      |      | + |
1884  *  |      |      |/| +
1885  *  +------+------+ |/
1886  *  |      |      | +
1887  *  |      |      |/
1888  *  +------+------+
1889  *  
1890  */
1891 static void subdivide_metaball_octal_node(octal_node *node, float size_x, float size_y, float size_z, short depth)
1892 {
1893         MetaElem *ml;
1894         ml_pointer *ml_p;
1895         float x, y, z;
1896         int a, i;
1897
1898         /* create new nodes */
1899         for (a = 0; a < 8; a++) {
1900                 node->nodes[a] = MEM_mallocN(sizeof(octal_node), "octal_node");
1901                 for (i = 0; i < 8; i++)
1902                         node->nodes[a]->nodes[i] = NULL;
1903                 node->nodes[a]->parent = node;
1904                 node->nodes[a]->elems.first = NULL;
1905                 node->nodes[a]->elems.last = NULL;
1906                 node->nodes[a]->count = 0;
1907                 node->nodes[a]->neg = 0;
1908                 node->nodes[a]->pos = 0;
1909         }
1910
1911         size_x /= 2;
1912         size_y /= 2;
1913         size_z /= 2;
1914         
1915         /* center of node */
1916         node->x = x = node->x_min + size_x;
1917         node->y = y = node->y_min + size_y;
1918         node->z = z = node->z_min + size_z;
1919
1920         /* setting up of border points of new nodes */
1921         node->nodes[0]->x_min = node->x_min;
1922         node->nodes[0]->y_min = node->y_min;
1923         node->nodes[0]->z_min = node->z_min;
1924         node->nodes[0]->x = node->nodes[0]->x_min + size_x / 2;
1925         node->nodes[0]->y = node->nodes[0]->y_min + size_y / 2;
1926         node->nodes[0]->z = node->nodes[0]->z_min + size_z / 2;
1927         
1928         node->nodes[1]->x_min = x;
1929         node->nodes[1]->y_min = node->y_min;
1930         node->nodes[1]->z_min = node->z_min;
1931         node->nodes[1]->x = node->nodes[1]->x_min + size_x / 2;
1932         node->nodes[1]->y = node->nodes[1]->y_min + size_y / 2;
1933         node->nodes[1]->z = node->nodes[1]->z_min + size_z / 2;
1934
1935         node->nodes[2]->x_min = x;
1936         node->nodes[2]->y_min = y;
1937         node->nodes[2]->z_min = node->z_min;
1938         node->nodes[2]->x = node->nodes[2]->x_min + size_x / 2;
1939         node->nodes[2]->y = node->nodes[2]->y_min + size_y / 2;
1940         node->nodes[2]->z = node->nodes[2]->z_min + size_z / 2;
1941
1942         node->nodes[3]->x_min = node->x_min;
1943         node->nodes[3]->y_min = y;
1944         node->nodes[3]->z_min = node->z_min;
1945         node->nodes[3]->x = node->nodes[3]->x_min + size_x / 2;
1946         node->nodes[3]->y = node->nodes[3]->y_min + size_y / 2;
1947         node->nodes[3]->z = node->nodes[3]->z_min + size_z / 2;
1948
1949         node->nodes[4]->x_min = node->x_min;
1950         node->nodes[4]->y_min = node->y_min;
1951         node->nodes[4]->z_min = z;
1952         node->nodes[4]->x = node->nodes[4]->x_min + size_x / 2;
1953         node->nodes[4]->y = node->nodes[4]->y_min + size_y / 2;
1954         node->nodes[4]->z = node->nodes[4]->z_min + size_z / 2;
1955         
1956         node->nodes[5]->x_min = x;
1957         node->nodes[5]->y_min = node->y_min;
1958         node->nodes[5]->z_min = z;
1959         node->nodes[5]->x = node->nodes[5]->x_min + size_x / 2;
1960         node->nodes[5]->y = node->nodes[5]->y_min + size_y / 2;
1961         node->nodes[5]->z = node->nodes[5]->z_min + size_z / 2;
1962
1963         node->nodes[6]->x_min = x;
1964         node->nodes[6]->y_min = y;
1965         node->nodes[6]->z_min = z;
1966         node->nodes[6]->x = node->nodes[6]->x_min + size_x / 2;
1967         node->nodes[6]->y = node->nodes[6]->y_min + size_y / 2;
1968         node->nodes[6]->z = node->nodes[6]->z_min + size_z / 2;
1969
1970         node->nodes[7]->x_min = node->x_min;
1971         node->nodes[7]->y_min = y;
1972         node->nodes[7]->z_min = z;
1973         node->nodes[7]->x = node->nodes[7]->x_min + size_x / 2;
1974         node->nodes[7]->y = node->nodes[7]->y_min + size_y / 2;
1975         node->nodes[7]->z = node->nodes[7]->z_min + size_z / 2;
1976
1977         ml_p = node->elems.first;
1978         
1979         /* setting up references of MetaElems for new nodes */
1980         while (ml_p) {
1981                 ml = ml_p->ml;
1982                 if (ml->bb->vec[0][2] < z) {
1983                         if (ml->bb->vec[0][1] < y) {
1984                                 /* vec[0][0] lies in first octant */
1985                                 if (ml->bb->vec[0][0] < x) {
1986                                         /* ml belongs to the (0)1st node */
1987                                         fill_metaball_octal_node(node, ml, 0);
1988
1989                                         /* ml belongs to the (3)4th node */
1990                                         if (ml->bb->vec[6][1] >= y) {
1991                                                 fill_metaball_octal_node(node, ml, 3);
1992
1993                                                 /* ml belongs to the (7)8th node */
1994                                                 if (ml->bb->vec[6][2] >= z) {
1995                                                         fill_metaball_octal_node(node, ml, 7);
1996                                                 }
1997                                         }
1998         
1999                                         /* ml belongs to the (1)2nd node */
2000                                         if (ml->bb->vec[6][0] >= x) {
2001                                                 fill_metaball_octal_node(node, ml, 1);
2002
2003                                                 /* ml belongs to the (5)6th node */
2004                                                 if (ml->bb->vec[6][2] >= z) {
2005                                                         fill_metaball_octal_node(node, ml, 5);
2006                                                 }
2007                                         }
2008
2009                                         /* ml belongs to the (2)3th node */
2010                                         if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
2011                                                 fill_metaball_octal_node(node, ml, 2);
2012                                                 
2013                                                 /* ml belong to the (6)7th node */
2014                                                 if (ml->bb->vec[6][2] >= z) {
2015                                                         fill_metaball_octal_node(node, ml, 6);
2016                                                 }
2017                                                 
2018                                         }
2019                         
2020                                         /* ml belongs to the (4)5th node too */
2021                                         if (ml->bb->vec[6][2] >= z) {
2022                                                 fill_metaball_octal_node(node, ml, 4);
2023                                         }
2024
2025                                         
2026                                         
2027                                 }
2028                                 /* vec[0][0] is in the (1)second octant */
2029                                 else {
2030                                         /* ml belong to the (1)2nd node */
2031                                         fill_metaball_octal_node(node, ml, 1);
2032
2033                                         /* ml belongs to the (2)3th node */
2034                                         if (ml->bb->vec[6][1] >= y) {
2035                                                 fill_metaball_octal_node(node, ml, 2);
2036
2037                                                 /* ml belongs to the (6)7th node */
2038                                                 if (ml->bb->vec[6][2] >= z) {
2039                                                         fill_metaball_octal_node(node, ml, 6);
2040                                                 }
2041                                                 
2042                                         }
2043                                         
2044                                         /* ml belongs to the (5)6th node */
2045                                         if (ml->bb->vec[6][2] >= z) {
2046                                                 fill_metaball_octal_node(node, ml, 5);
2047                                         }
2048                                 }
2049                         }
2050                         else {
2051                                 /* vec[0][0] is in the (3)4th octant */
2052                                 if (ml->bb->vec[0][0] < x) {
2053                                         /* ml belongs to the (3)4nd node */
2054                                         fill_metaball_octal_node(node, ml, 3);
2055                                         
2056                                         /* ml belongs to the (7)8th node */
2057                                         if (ml->bb->vec[6][2] >= z) {
2058                                                 fill_metaball_octal_node(node, ml, 7);
2059                                         }
2060                                 
2061
2062                                         /* ml belongs to the (2)3th node */
2063                                         if (ml->bb->vec[6][0] >= x) {
2064                                                 fill_metaball_octal_node(node, ml, 2);
2065                                         
2066                                                 /* ml belongs to the (6)7th node */
2067                                                 if (ml->bb->vec[6][2] >= z) {
2068                                                         fill_metaball_octal_node(node, ml, 6);
2069                                                 }
2070                                         }
2071                                 }
2072
2073                         }
2074
2075                         /* vec[0][0] is in the (2)3th octant */
2076                         if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
2077                                 /* ml belongs to the (2)3th node */
2078                                 fill_metaball_octal_node(node, ml, 2);
2079                                 
2080                                 /* ml belongs to the (6)7th node */
2081                                 if (ml->bb->vec[6][2] >= z) {
2082                                         fill_metaball_octal_node(node, ml, 6);
2083                                 }
2084                         }
2085                 }
2086                 else {
2087                         if (ml->bb->vec[0][1] < y) {
2088                                 /* vec[0][0] lies in (4)5th octant */
2089                                 if (ml->bb->vec[0][0] < x) {
2090                                         /* ml belongs to the (4)5th node */
2091                                         fill_metaball_octal_node(node, ml, 4);
2092
2093                                         if (ml->bb->vec[6][0] >= x) {
2094                                                 fill_metaball_octal_node(node, ml, 5);
2095                                         }
2096
2097                                         if (ml->bb->vec[6][1] >= y) {
2098                                                 fill_metaball_octal_node(node, ml, 7);
2099                                         }
2100                                         
2101                                         if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
2102                                                 fill_metaball_octal_node(node, ml, 6);
2103                                         }
2104                                 }
2105                                 /* vec[0][0] lies in (5)6th octant */
2106                                 else {
2107                                         fill_metaball_octal_node(node, ml, 5);
2108
2109                                         if (ml->bb->vec[6][1] >= y) {
2110                                                 fill_metaball_octal_node(node, ml, 6);
2111                                         }
2112                                 }
2113                         }
2114                         else {
2115                                 /* vec[0][0] lies in (7)8th octant */
2116                                 if (ml->bb->vec[0][0] < x) {
2117                                         fill_metaball_octal_node(node, ml, 7);
2118
2119                                         if (ml->bb->vec[6][0] >= x) {
2120                                                 fill_metaball_octal_node(node, ml, 6);
2121                                         }
2122                                 }
2123
2124                         }
2125                         
2126                         /* vec[0][0] lies in (6)7th octant */
2127                         if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
2128                                 fill_metaball_octal_node(node, ml, 6);
2129                         }
2130                 }
2131                 ml_p = ml_p->next;
2132         }
2133
2134         /* free references of MetaElems for curent node (it is not needed anymore) */
2135         BLI_freelistN(&node->elems);
2136
2137         depth--;
2138         
2139         if (depth > 0) {
2140                 for (a = 0; a < 8; a++) {
2141                         if (node->nodes[a]->count > 0) /* if node is not empty, then it is subdivided */
2142                                 subdivide_metaball_octal_node(node->nodes[a], size_x, size_y, size_z, depth);
2143                 }
2144         }
2145 }
2146
2147 /* free all octal nodes recursively */
2148 static void free_metaball_octal_node(octal_node *node)
2149 {
2150         int a;
2151         for (a = 0; a < 8; a++) {
2152                 if (node->nodes[a] != NULL) free_metaball_octal_node(node->nodes[a]);
2153         }
2154         BLI_freelistN(&node->elems);
2155         MEM_freeN(node);
2156 }
2157
2158 /* If scene include more than one MetaElem, then octree is used */
2159 static void init_metaball_octal_tree(PROCESS *process, int depth)
2160 {
2161         struct octal_node *node;
2162         ml_pointer *ml_p;
2163         float size[3];
2164         int a;
2165         
2166         process->metaball_tree = MEM_mallocN(sizeof(octal_tree), "metaball_octal_tree");
2167         process->metaball_tree->first = node = MEM_mallocN(sizeof(octal_node), "metaball_octal_node");
2168         /* maximal depth of octree */
2169         process->metaball_tree->depth = depth;
2170
2171         process->metaball_tree->neg = node->neg = 0;
2172         process->metaball_tree->pos = node->pos = 0;
2173         
2174         node->elems.first = NULL;
2175         node->elems.last = NULL;
2176         node->count = 0;
2177
2178         for (a = 0; a < 8; a++)
2179                 node->nodes[a] = NULL;
2180
2181         node->x_min = node->y_min = node->z_min = FLT_MAX;
2182         node->x_max = node->y_max = node->z_max = -FLT_MAX;
2183
2184         /* size of octal tree scene */
2185         for (a = 0; a < process->totelem; a++) {
2186                 if (process->mainb[a]->bb->vec[0][0] < node->x_min) node->x_min = process->mainb[a]->bb->vec[0][0];
2187                 if (process->mainb[a]->bb->vec[0][1] < node->y_min) node->y_min = process->mainb[a]->bb->vec[0][1];
2188                 if (process->mainb[a]->bb->vec[0][2] < node->z_min) node->z_min = process->mainb[a]->bb->vec[0][2];
2189
2190                 if (process->mainb[a]->bb->vec[6][0] > node->x_max) node->x_max = process->mainb[a]->bb->vec[6][0];
2191                 if (process->mainb[a]->bb->vec[6][1] > node->y_max) node->y_max = process->mainb[a]->bb->vec[6][1];
2192                 if (process->mainb[a]->bb->vec[6][2] > node->z_max) node->z_max = process->mainb[a]->bb->vec[6][2];
2193
2194                 ml_p = MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
2195                 ml_p->ml = process->mainb[a];
2196                 BLI_addtail(&node->elems, ml_p);
2197
2198                 if ((process->mainb[a]->flag & MB_NEGATIVE) == 0) {
2199                         /* number of positive MetaElem in scene */
2200                         process->metaball_tree->pos++;
2201                 }
2202                 else {
2203                         /* number of negative MetaElem in scene */
2204                         process->metaball_tree->neg++;
2205                 }
2206         }
2207
2208         /* size of first node */
2209         size[0] = node->x_max - node->x_min;
2210         size[1] = node->y_max - node->y_min;
2211         size[2] = node->z_max - node->z_min;
2212
2213         /* first node is subdivided recursively */
2214         subdivide_metaball_octal_node(node, size[0], size[1], size[2], process->metaball_tree->depth);
2215 }
2216
2217 static void mball_count(EvaluationContext *eval_ctx, PROCESS *process, Scene *scene, Object *basis)
2218 {
2219         Scene *sce_iter = scene;
2220         Base *base;
2221         Object *ob, *bob = basis;
2222         MetaElem *ml = NULL;
2223         int basisnr, obnr;
2224         char basisname[MAX_ID_NAME], obname[MAX_ID_NAME];
2225         SceneBaseIter iter;
2226
2227         BLI_split_name_num(basisname, &basisnr, basis->id.name + 2, '.');
2228         process->totelem = 0;
2229
2230         /* XXX recursion check, see scene.c, just too simple code this BKE_scene_base_iter_next() */
2231         if (F_ERROR == BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 0, NULL, NULL))
2232                 return;
2233
2234         while (BKE_scene_base_iter_next(eval_ctx, &iter, &sce_iter, 1, &base, &ob)) {
2235                 if (ob->type == OB_MBALL) {
2236                         if (ob == bob) {
2237                                 MetaBall *mb = ob->data;
2238
2239                                 /* if bob object is in edit mode, then dynamic list of all MetaElems
2240                                  * is stored in editelems */
2241                                 if (mb->editelems) ml = mb->editelems->first;
2242                                 /* if bob object is in object mode */
2243                                 else ml = mb->elems.first;
2244                         }
2245                         else {
2246                                 BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
2247
2248                                 /* object ob has to be in same "group" ... it means, that it has to have
2249                                  * same base of its name */
2250                                 if (strcmp(obname, basisname) == 0) {
2251                                         MetaBall *mb = ob->data;
2252
2253                                         /* if object is in edit mode, then dynamic list of all MetaElems
2254                                          * is stored in editelems */
2255                                         if (mb->editelems) ml = mb->editelems->first;
2256                                         /* if bob object is in object mode */
2257                                         else ml = mb->elems.first;
2258                                 }
2259                         }
2260
2261                         for ( ; ml; ml = ml->next) {
2262                                 if (!(ml->flag & MB_HIDE)) {
2263                                         process->totelem++;
2264                                 }
2265                         }
2266                 }
2267         }
2268 }
2269
2270 void BKE_mball_polygonize(EvaluationContext *eval_ctx, Scene *scene, Object *ob, ListBase *dispbase)
2271 {
2272         MetaBall *mb;
2273         DispList *dl;
2274         int a, nr_cubes;
2275         float *co, *no, totsize, width;
2276         PROCESS process = {0};
2277
2278         mb = ob->data;
2279
2280         mball_count(eval_ctx, &process, scene, ob);
2281
2282         if (process.totelem == 0) return;
2283         if ((eval_ctx->for_render == false) && (mb->flag == MB_UPDATE_NEVER)) return;
2284         if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_FAST) return;
2285
2286         process.thresh = mb->thresh;
2287
2288         /* total number of MetaElems (totelem) is precomputed in find_basis_mball() function */
2289         process.mainb = MEM_mallocN(sizeof(void *) * process.totelem, "mainb");
2290         
2291         /* initialize all mainb (MetaElems) */
2292         totsize = init_meta(eval_ctx, &process, scene, ob);
2293
2294         /* if scene includes more than one MetaElem, then octal tree optimization is used */
2295         if ((process.totelem >    1) && (process.totelem <=   64)) init_metaball_octal_tree(&process, 1);
2296         if ((process.totelem >   64) && (process.totelem <=  128)) init_metaball_octal_tree(&process, 2);
2297         if ((process.totelem >  128) && (process.totelem <=  512)) init_metaball_octal_tree(&process, 3);
2298         if ((process.totelem >  512) && (process.totelem <= 1024)) init_metaball_octal_tree(&process, 4);
2299         if (process.totelem  > 1024)                               init_metaball_octal_tree(&process, 5);
2300
2301         /* don't polygonize metaballs with too high resolution (base mball to small)
2302          * note: Eps was 0.0001f but this was giving problems for blood animation for durian, using 0.00001f */
2303         if (process.metaball_tree) {
2304                 if (ob->size[0] <= 0.00001f * (process.metaball_tree->first->x_max - process.metaball_tree->first->x_min) ||
2305                     ob->size[1] <= 0.00001f * (process.metaball_tree->first->y_max - process.metaball_tree->first->y_min) ||
2306                     ob->size[2] <= 0.00001f * (process.metaball_tree->first->z_max - process.metaball_tree->first->z_min))
2307                 {
2308                         new_pgn_element(&process, -1); /* free values created by init_meta */
2309
2310                         MEM_freeN(process.mainb);
2311
2312                         /* free tree */
2313                         free_metaball_octal_node(process.metaball_tree->first);
2314                         MEM_freeN(process.metaball_tree);
2315
2316                         return;
2317                 }
2318         }
2319
2320         /* width is size per polygonize cube */
2321         if (eval_ctx->for_render) {
2322                 width = mb->rendersize;
2323         }
2324         else {
2325                 width = mb->wiresize;
2326                 if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_HALFRES) {
2327                         width *= 2;
2328                 }
2329         }
2330         /* nr_cubes is just for safety, minimum is totsize */
2331         nr_cubes = (int)(0.5f + totsize / width);
2332
2333         /* init process */
2334         process.function = metaball;
2335         process.size = width;
2336         process.bounds = nr_cubes;
2337         process.cubes = NULL;
2338         process.delta = width / (float)(RES * RES);
2339
2340         polygonize(&process, mb);
2341         
2342         MEM_freeN(process.mainb);
2343
2344         /* free octal tree */
2345         if (process.totelem > 1) {
2346                 free_metaball_octal_node(process.metaball_tree->first);
2347                 MEM_freeN(process.metaball_tree);
2348                 process.metaball_tree = NULL;
2349         }
2350
2351         if (process.curindex) {
2352                 VERTEX *ptr = process.vertices.ptr;
2353
2354                 dl = MEM_callocN(sizeof(DispList), "mbaldisp");
2355                 BLI_addtail(dispbase, dl);
2356                 dl->type = DL_INDEX4;
2357                 dl->nr = process.vertices.count;
2358                 dl->parts = process.curindex;
2359
2360                 dl->index = process.indices;
2361                 process.indices = NULL;
2362
2363                 a = process.vertices.count;
2364                 dl->verts = co = MEM_mallocN(sizeof(float) * 3 * a, "mballverts");
2365                 dl->nors = no = MEM_mallocN(sizeof(float) * 3 * a, "mballnors");
2366
2367                 for (a = 0; a < process.vertices.count; ptr++, a++, no += 3, co += 3) {
2368                         copy_v3_v3(co, ptr->co);
2369                         copy_v3_v3(no, ptr->no);
2370                 }
2371         }
2372
2373         freepolygonize(&process);
2374 }
2375
2376 bool BKE_mball_minmax_ex(MetaBall *mb, float min[3], float max[3],
2377                          float obmat[4][4], const short flag)
2378 {
2379         const float scale = obmat ? mat4_to_scale(obmat) : 1.0f;
2380         MetaElem *ml;
2381         bool changed = false;
2382         float centroid[3], vec[3];
2383
2384         INIT_MINMAX(min, max);
2385
2386         for (ml = mb->elems.first; ml; ml = ml->next) {
2387                 if ((ml->flag & flag) == flag) {
2388                         const float scale_mb = (ml->rad * 0.5f) * scale;
2389                         int i;
2390
2391                         if (obmat) {
2392                                 mul_v3_m4v3(centroid, obmat, &ml->x);
2393                         }
2394                         else {
2395                                 copy_v3_v3(centroid, &ml->x);
2396                         }
2397
2398                         /* TODO, non circle shapes cubes etc, probably nobody notices - campbell */
2399                         for (i = -1; i != 3; i += 2) {
2400                                 copy_v3_v3(vec, centroid);
2401                                 add_v3_fl(vec, scale_mb * i);
2402                                 minmax_v3v3_v3(min, max, vec);
2403                         }
2404                         changed = true;
2405                 }
2406         }
2407
2408         return changed;
2409 }
2410
2411
2412 /* basic vertex data functions */
2413 bool BKE_mball_minmax(MetaBall *mb, float min[3], float max[3])
2414 {
2415         MetaElem *ml;
2416
2417         INIT_MINMAX(min, max);
2418
2419         for (ml = mb->elems.first; ml; ml = ml->next) {
2420                 minmax_v3v3_v3(min, max, &ml->x);
2421         }
2422
2423         return (mb->elems.first != NULL);
2424 }
2425
2426 bool BKE_mball_center_median(MetaBall *mb, float r_cent[3])
2427 {
2428         MetaElem *ml;
2429         int total = 0;
2430
2431         zero_v3(r_cent);
2432
2433         for (ml = mb->elems.first; ml; ml = ml->next) {
2434                 add_v3_v3(r_cent, &ml->x);
2435                 total++;
2436         }
2437
2438         if (total) {
2439                 mul_v3_fl(r_cent, 1.0f / (float)total);
2440         }
2441
2442         return (total != 0);
2443 }
2444
2445 bool BKE_mball_center_bounds(MetaBall *mb, float r_cent[3])
2446 {
2447         float min[3], max[3];
2448
2449         if (BKE_mball_minmax(mb, min, max)) {
2450                 mid_v3_v3v3(r_cent, min, max);
2451                 return 1;
2452         }
2453
2454         return 0;
2455 }
2456
2457 void BKE_mball_translate(MetaBall *mb, const float offset[3])
2458 {
2459         MetaElem *ml;
2460
2461         for (ml = mb->elems.first; ml; ml = ml->next) {
2462                 add_v3_v3(&ml->x, offset);
2463         }
2464 }
2465
2466 /* *** select funcs *** */
2467 void BKE_mball_select_all(struct MetaBall *mb)
2468 {
2469         MetaElem *ml;
2470
2471         for (ml = mb->editelems->first; ml; ml = ml->next) {
2472                 ml->flag |= SELECT;
2473         }
2474 }
2475
2476 void BKE_mball_deselect_all(MetaBall *mb)
2477 {
2478         MetaElem *ml;
2479
2480         for (ml = mb->editelems->first; ml; ml = ml->next) {
2481                 ml->flag &= ~SELECT;
2482         }
2483 }
2484
2485 void BKE_mball_select_swap(struct MetaBall *mb)
2486 {
2487         MetaElem *ml;
2488
2489         for (ml = mb->editelems->first; ml; ml = ml->next) {
2490                 ml->flag ^= SELECT;
2491         }
2492 }
2493