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