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