8b2e66812c7073714ba911fa6d6a353f75725ba6
[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_metaball_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_metaball_free(MetaBall *mb)
189 {
190         BKE_metaball_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) freedisplist(&mb->disp);
200 }
201
202 MetaBall *BKE_metaball_add(const char *name)
203 {
204         MetaBall *mb;
205         
206         mb= alloc_libblock(&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_metaball_copy(MetaBall *mb)
219 {
220         MetaBall *mbn;
221         int a;
222         
223         mbn= copy_libblock(&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_metaball_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_metaball_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_metaball_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_metaball_tex_space_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         boundbox_set_from_min_max(bb, min, max);
391 }
392
393 float *BKE_metaball_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_metaball_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_metaball_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_metaball_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_metaball_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 next_object() */
486         if (F_ERROR==next_object(&sce_iter, 0, NULL, NULL))
487                 return;
488         
489         while (next_object(&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_metaball_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 next_object() */
532         if (F_ERROR==next_object(&sce_iter, 0, NULL, NULL))
533                 return NULL;
534         
535         while (next_object(&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 static int corner2[12]     = {
924         LBF, LTF, LTN, LTF, RBF, RTF, RTN, RTF, RBN, RBF, RTN, RTF};
925 static int leftface[12]    = {
926         B,  L,  L,  F,  R,  T,  N,  R,  N,  B,  T,  F};
927 /* face on left when going corner1 to corner2 */
928 static int rightface[12]   = {
929         L,  T,  N,  L,  B,  R,  R,  F,  B,  F,  N,  T};
930 /* face on right when going corner1 to corner2 */
931
932
933 /* docube: triangulate the cube directly, without decomposition */
934
935 static void docube(CUBE *cube, PROCESS *p, MetaBall *mb)
936 {
937         INTLISTS *polys;
938         CORNER *c1, *c2;
939         int i, index = 0, count, indexar[8];
940         
941         for (i = 0; i < 8; i++) if (cube->corners[i]->value > 0.0f) index += (1<<i);
942         
943         for (polys = cubetable[index]; polys; polys = polys->next) {
944                 INTLIST *edges;
945                 
946                 count = 0;
947                 
948                 for (edges = polys->list; edges; edges = edges->next) {
949                         c1 = cube->corners[corner1[edges->i]];
950                         c2 = cube->corners[corner2[edges->i]];
951                         
952                         indexar[count] = vertid(c1, c2, p, mb);
953                         count++;
954                 }
955                 if (count>2) {
956                         switch (count) {
957                         case 3:
958                                 accum_mballfaces(indexar[2], indexar[1], indexar[0], 0);
959                                 break;
960                         case 4:
961                                 if (indexar[0]==0) accum_mballfaces(indexar[0], indexar[3], indexar[2], indexar[1]);
962                                 else accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]);
963                                 break;
964                         case 5:
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
968                                 accum_mballfaces(indexar[4], indexar[3], indexar[0], 0);
969                                 break;
970                         case 6:
971                                 if (indexar[0]==0) {
972                                         accum_mballfaces(indexar[0], indexar[3], indexar[2], indexar[1]);
973                                         accum_mballfaces(indexar[0], indexar[5], indexar[4], indexar[3]);
974                                 }
975                                 else {
976                                         accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]);
977                                         accum_mballfaces(indexar[5], indexar[4], indexar[3], indexar[0]);
978                                 }
979                                 break;
980                         case 7:
981                                 if (indexar[0]==0) {
982                                         accum_mballfaces(indexar[0], indexar[3], indexar[2], indexar[1]);
983                                         accum_mballfaces(indexar[0], indexar[5], indexar[4], indexar[3]);
984                                 }
985                                 else {
986                                         accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]);
987                                         accum_mballfaces(indexar[5], indexar[4], indexar[3], indexar[0]);
988                                 }
989                                 
990                                 accum_mballfaces(indexar[6], indexar[5], indexar[0], 0);
991                                 
992                                 break;
993                         }
994                 }
995         }
996 }
997
998
999 /* testface: given cube at lattice (i, j, k), and four corners of face,
1000  * if surface crosses face, compute other four corners of adjacent cube
1001  * and add new cube to cube stack */
1002
1003 static void testface(int i, int j, int k, CUBE* old, int bit, int c1, int c2, int c3, int c4, PROCESS *p)
1004 {
1005         CUBE newc;
1006         CUBES *oldcubes = p->cubes;
1007         CORNER *corn1, *corn2, *corn3, *corn4;
1008         int n, pos;
1009
1010         corn1= old->corners[c1];
1011         corn2= old->corners[c2];
1012         corn3= old->corners[c3];
1013         corn4= old->corners[c4];
1014         
1015         pos = corn1->value > 0.0f ? 1 : 0;
1016
1017         /* test if no surface crossing */
1018         if ( (corn2->value > 0) == pos && (corn3->value > 0) == pos && (corn4->value > 0) == pos) return;
1019         /* test if cube out of bounds */
1020         /*if ( abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) return;*/
1021         /* test if already visited (always as last) */
1022         if (setcenter(p->centers, i, j, k)) return;
1023
1024
1025         /* create new cube and add cube to top of stack: */
1026         p->cubes = (CUBES *) new_pgn_element(sizeof(CUBES));
1027         p->cubes->next = oldcubes;
1028         
1029         newc.i = i;
1030         newc.j = j;
1031         newc.k = k;
1032         for (n = 0; n < 8; n++) newc.corners[n] = NULL;
1033         
1034         newc.corners[FLIP(c1, bit)] = corn1;
1035         newc.corners[FLIP(c2, bit)] = corn2;
1036         newc.corners[FLIP(c3, bit)] = corn3;
1037         newc.corners[FLIP(c4, bit)] = corn4;
1038
1039         if (newc.corners[0]==NULL) newc.corners[0] = setcorner(p, i, j, k);
1040         if (newc.corners[1]==NULL) newc.corners[1] = setcorner(p, i, j, k+1);
1041         if (newc.corners[2]==NULL) newc.corners[2] = setcorner(p, i, j+1, k);
1042         if (newc.corners[3]==NULL) newc.corners[3] = setcorner(p, i, j+1, k+1);
1043         if (newc.corners[4]==NULL) newc.corners[4] = setcorner(p, i+1, j, k);
1044         if (newc.corners[5]==NULL) newc.corners[5] = setcorner(p, i+1, j, k+1);
1045         if (newc.corners[6]==NULL) newc.corners[6] = setcorner(p, i+1, j+1, k);
1046         if (newc.corners[7]==NULL) newc.corners[7] = setcorner(p, i+1, j+1, k+1);
1047
1048         p->cubes->cube= newc;   
1049 }
1050
1051 /* setcorner: return corner with the given lattice location
1052  * set (and cache) its function value */
1053
1054 static CORNER *setcorner (PROCESS* p, int i, int j, int k)
1055 {
1056         /* for speed, do corner value caching here */
1057         CORNER *c;
1058         int index;
1059
1060         /* does corner exist? */
1061         index = HASH(i, j, k);
1062         c = p->corners[index];
1063         
1064         for (; c != NULL; c = c->next) {
1065                 if (c->i == i && c->j == j && c->k == k) {
1066                         return c;
1067                 }
1068         }
1069
1070         c = (CORNER *) new_pgn_element(sizeof(CORNER));
1071
1072         c->i = i; 
1073         c->x = ((float)i-0.5f)*p->size;
1074         c->j = j; 
1075         c->y = ((float)j-0.5f)*p->size;
1076         c->k = k; 
1077         c->z = ((float)k-0.5f)*p->size;
1078         c->value = p->function(c->x, c->y, c->z);
1079         
1080         c->next = p->corners[index];
1081         p->corners[index] = c;
1082         
1083         return c;
1084 }
1085
1086
1087 /* nextcwedge: return next clockwise edge from given edge around given face */
1088
1089 static int nextcwedge (int edge, int face)
1090 {
1091         switch (edge) {
1092         case LB: 
1093                 return (face == L)? LF : BN;
1094         case LT: 
1095                 return (face == L)? LN : TF;
1096         case LN: 
1097                 return (face == L)? LB : TN;
1098         case LF: 
1099                 return (face == L)? LT : BF;
1100         case RB: 
1101                 return (face == R)? RN : BF;
1102         case RT: 
1103                 return (face == R)? RF : TN;
1104         case RN: 
1105                 return (face == R)? RT : BN;
1106         case RF: 
1107                 return (face == R)? RB : TF;
1108         case BN: 
1109                 return (face == B)? RB : LN;
1110         case BF: 
1111                 return (face == B)? LB : RF;
1112         case TN: 
1113                 return (face == T)? LT : RN;
1114         case TF: 
1115                 return (face == T)? RT : LF;
1116         }
1117         return 0;
1118 }
1119
1120
1121 /* otherface: return face adjoining edge that is not the given face */
1122
1123 static int otherface (int edge, int face)
1124 {
1125         int other = leftface[edge];
1126         return face == other? rightface[edge] : other;
1127 }
1128
1129
1130 /* makecubetable: create the 256 entry table for cubical polygonization */
1131
1132 static void makecubetable (void)
1133 {
1134         static int isdone= 0;
1135         int i, e, c, done[12], pos[8];
1136
1137         if (isdone) return;
1138         isdone= 1;
1139
1140         for (i = 0; i < 256; i++) {
1141                 for (e = 0; e < 12; e++) done[e] = 0;
1142                 for (c = 0; c < 8; c++) pos[c] = MB_BIT(i, c);
1143                 for (e = 0; e < 12; e++)
1144                         if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
1145                                 INTLIST *ints = NULL;
1146                                 INTLISTS *lists = (INTLISTS *) MEM_callocN(sizeof(INTLISTS), "mball_intlist");
1147                                 int start = e, edge = e;
1148                                 
1149                                 /* get face that is to right of edge from pos to neg corner: */
1150                                 int face = pos[corner1[e]]? rightface[e] : leftface[e];
1151                                 
1152                                 while (1) {
1153                                         edge = nextcwedge(edge, face);
1154                                         done[edge] = 1;
1155                                         if (pos[corner1[edge]] != pos[corner2[edge]]) {
1156                                                 INTLIST *tmp = ints;
1157                                                 
1158                                                 ints = (INTLIST *) MEM_callocN(sizeof(INTLIST), "mball_intlist");
1159                                                 ints->i = edge;
1160                                                 ints->next = tmp; /* add edge to head of list */
1161                                                 
1162                                                 if (edge == start) break;
1163                                                 face = otherface(edge, face);
1164                                         }
1165                                 }
1166                                 lists->list = ints; /* add ints to head of table entry */
1167                                 lists->next = cubetable[i];
1168                                 cubetable[i] = lists;
1169                         }
1170         }
1171 }
1172
1173 void BKE_metaball_cubeTable_free(void)
1174 {
1175         int i;
1176         INTLISTS *lists, *nlists;
1177         INTLIST *ints, *nints;
1178
1179         for (i = 0; i < 256; i++) {
1180                 lists= cubetable[i];
1181                 while (lists) {
1182                         nlists= lists->next;
1183                         
1184                         ints= lists->list;
1185                         while (ints) {
1186                                 nints= ints->next;
1187                                 MEM_freeN(ints);
1188                                 ints= nints;
1189                         }
1190                         
1191                         MEM_freeN(lists);
1192                         lists= nlists;
1193                 }
1194                 cubetable[i]= NULL;
1195         }
1196 }
1197
1198 /**** Storage ****/
1199
1200 /* setcenter: set (i, j, k) entry of table[]
1201  * return 1 if already set; otherwise, set and return 0 */
1202
1203 static int setcenter(CENTERLIST *table[], int i, int j, int k)
1204 {
1205         int index;
1206         CENTERLIST *newc, *l, *q;
1207
1208         index= HASH(i, j, k);
1209         q= table[index];
1210
1211         for (l = q; l != NULL; l = l->next) {
1212                 if (l->i == i && l->j == j && l->k == k) return 1;
1213         }
1214         
1215         newc = (CENTERLIST *) new_pgn_element(sizeof(CENTERLIST));
1216         newc->i = i; 
1217         newc->j = j; 
1218         newc->k = k; 
1219         newc->next = q;
1220         table[index] = newc;
1221         
1222         return 0;
1223 }
1224
1225
1226 /* setedge: set vertex id for edge */
1227
1228 static void setedge (EDGELIST *table[],
1229                           int i1, int j1,
1230                           int k1, int i2,
1231                           int j2, int k2,
1232                           int vid)
1233 {
1234         unsigned int index;
1235         EDGELIST *newe;
1236         
1237         if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) {
1238                 int t=i1; 
1239                 i1=i2; 
1240                 i2=t; 
1241                 t=j1; 
1242                 j1=j2; 
1243                 j2=t; 
1244                 t=k1; 
1245                 k1=k2; 
1246                 k2=t;
1247         }
1248         index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
1249         newe = (EDGELIST *) new_pgn_element(sizeof(EDGELIST));
1250         newe->i1 = i1; 
1251         newe->j1 = j1; 
1252         newe->k1 = k1;
1253         newe->i2 = i2; 
1254         newe->j2 = j2; 
1255         newe->k2 = k2;
1256         newe->vid = vid;
1257         newe->next = table[index];
1258         table[index] = newe;
1259 }
1260
1261
1262 /* getedge: return vertex id for edge; return -1 if not set */
1263
1264 static int getedge (EDGELIST *table[],
1265                          int i1, int j1, int k1,
1266                          int i2, int j2, int k2)
1267 {
1268         EDGELIST *q;
1269         
1270         if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) {
1271                 int t=i1; 
1272                 i1=i2; 
1273                 i2=t; 
1274                 t=j1; 
1275                 j1=j2; 
1276                 j2=t; 
1277                 t=k1; 
1278                 k1=k2; 
1279                 k2=t;
1280         }
1281         q = table[HASH(i1, j1, k1)+HASH(i2, j2, k2)];
1282         for (; q != NULL; q = q->next) {
1283                 if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 &&
1284                         q->i2 == i2 && q->j2 == j2 && q->k2 == k2)
1285                 {
1286                         return q->vid;
1287                 }
1288         }
1289         return -1;
1290 }
1291
1292
1293 /**** Vertices ****/
1294
1295 #undef R
1296
1297
1298
1299 /* vertid: return index for vertex on edge:
1300  * c1->value and c2->value are presumed of different sign
1301  * return saved index if any; else compute vertex and save */
1302
1303 /* addtovertices: add v to sequence of vertices */
1304
1305 static void addtovertices (VERTICES *vertices, VERTEX v)
1306 {
1307         if (vertices->count == vertices->max) {
1308                 int i;
1309                 VERTEX *newv;
1310                 vertices->max = vertices->count == 0 ? 10 : 2*vertices->count;
1311                 newv = (VERTEX *) MEM_callocN(vertices->max * sizeof(VERTEX), "addtovertices");
1312                 
1313                 for (i = 0; i < vertices->count; i++) newv[i] = vertices->ptr[i];
1314                 
1315                 if (vertices->ptr != NULL) MEM_freeN(vertices->ptr);
1316                 vertices->ptr = newv;
1317         }
1318         vertices->ptr[vertices->count++] = v;
1319 }
1320
1321 /* vnormal: compute unit length surface normal at point */
1322
1323 static void vnormal (MB_POINT *point, PROCESS *p, MB_POINT *v)
1324 {
1325         float delta= 0.2f*p->delta;
1326         float f = p->function(point->x, point->y, point->z);
1327
1328         v->x = p->function(point->x+delta, point->y, point->z)-f;
1329         v->y = p->function(point->x, point->y+delta, point->z)-f;
1330         v->z = p->function(point->x, point->y, point->z+delta)-f;
1331         f = sqrtf(v->x*v->x + v->y*v->y + v->z*v->z);
1332
1333         if (f != 0.0f) {
1334                 v->x /= f; 
1335                 v->y /= f; 
1336                 v->z /= f;
1337         }
1338         
1339         if (FALSE) {
1340                 MB_POINT temp;
1341                 
1342                 delta *= 2.0f;
1343                 
1344                 f = p->function(point->x, point->y, point->z);
1345         
1346                 temp.x = p->function(point->x+delta, point->y, point->z)-f;
1347                 temp.y = p->function(point->x, point->y+delta, point->z)-f;
1348                 temp.z = p->function(point->x, point->y, point->z+delta)-f;
1349                 f = sqrtf(temp.x*temp.x + temp.y*temp.y + temp.z*temp.z);
1350         
1351                 if (f != 0.0f) {
1352                         temp.x /= f; 
1353                         temp.y /= f; 
1354                         temp.z /= f;
1355                         
1356                         v->x+= temp.x;
1357                         v->y+= temp.y;
1358                         v->z+= temp.z;
1359                         
1360                         f = sqrtf(v->x*v->x + v->y*v->y + v->z*v->z);
1361                 
1362                         if (f != 0.0f) {
1363                                 v->x /= f; 
1364                                 v->y /= f; 
1365                                 v->z /= f;
1366                         }
1367                 }
1368         }
1369         
1370 }
1371
1372
1373 static int vertid (CORNER *c1, CORNER *c2, PROCESS *p, MetaBall *mb)
1374 {
1375         VERTEX v;
1376         MB_POINT a, b;
1377         int vid = getedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
1378
1379         if (vid != -1) return vid;                           /* previously computed */
1380         a.x = c1->x; 
1381         a.y = c1->y; 
1382         a.z = c1->z;
1383         b.x = c2->x; 
1384         b.y = c2->y; 
1385         b.z = c2->z;
1386
1387         converge(&a, &b, c1->value, c2->value, p->function, &v.position, mb, 1); /* position */
1388         vnormal(&v.position, p, &v.normal);
1389
1390         addtovertices(&p->vertices, v);                    /* save vertex */
1391         vid = p->vertices.count-1;
1392         setedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
1393         
1394         return vid;
1395 }
1396
1397
1398
1399
1400 /* converge: from two points of differing sign, converge to zero crossing */
1401 /* watch it: p1 and p2 are used to calculate */
1402 static void converge (MB_POINT *p1, MB_POINT *p2, float v1, float v2,
1403                            float (*function)(float, float, float), MB_POINT *p, MetaBall *mb, int f)
1404 {
1405         int i = 0;
1406         MB_POINT pos, neg;
1407         float positive = 0.0f, negative = 0.0f;
1408         float dx = 0.0f, dy = 0.0f, dz = 0.0f;
1409         
1410         if (v1 < 0) {
1411                 pos= *p2;
1412                 neg= *p1;
1413                 positive = v2;
1414                 negative = v1;
1415         }
1416         else {
1417                 pos= *p1;
1418                 neg= *p2;
1419                 positive = v1;
1420                 negative = v2;
1421         }
1422
1423         dx = pos.x - neg.x;
1424         dy = pos.y - neg.y;
1425         dz = pos.z - neg.z;
1426
1427 /* Approximation by linear interpolation is faster then binary subdivision,
1428  * but it results sometimes (mb->thresh < 0.2) into the strange results */
1429         if ((mb->thresh > 0.2f) && (f==1)) {
1430         if ((dy == 0.0f) && (dz == 0.0f)) {
1431                 p->x = neg.x - negative*dx/(positive-negative);
1432                 p->y = neg.y;
1433                 p->z = neg.z;
1434                 return;
1435         }
1436           if ((dx == 0.0f) && (dz == 0.0f)) {
1437                 p->x = neg.x;
1438                 p->y = neg.y - negative*dy/(positive-negative);
1439                 p->z = neg.z;
1440                 return;
1441         }
1442         if ((dx == 0.0f) && (dy == 0.0f)) {
1443                 p->x = neg.x;
1444                 p->y = neg.y;
1445                 p->z = neg.z - negative*dz/(positive-negative);
1446                 return;
1447         }
1448         }
1449
1450         if ((dy == 0.0f) && (dz == 0.0f)) {
1451                 p->y = neg.y;
1452                 p->z = neg.z;
1453                 while (1) {
1454                         if (i++ == RES) return;
1455                         p->x = 0.5f*(pos.x + neg.x);
1456                         if ((function(p->x, p->y, p->z)) > 0.0f)        pos.x = p->x; else neg.x = p->x;
1457                 }
1458         }
1459
1460         if ((dx == 0.0f) && (dz == 0.0f)) {
1461                 p->x = neg.x;
1462                 p->z = neg.z;
1463                 while (1) {
1464                         if (i++ == RES) return;
1465                         p->y = 0.5f*(pos.y + neg.y);
1466                         if ((function(p->x, p->y, p->z)) > 0.0f)        pos.y = p->y; else neg.y = p->y;
1467                 }
1468           }
1469    
1470         if ((dx == 0.0f) && (dy == 0.0f)) {
1471                 p->x = neg.x;
1472                 p->y = neg.y;
1473                 while (1) {
1474                         if (i++ == RES) return;
1475                         p->z = 0.5f*(pos.z + neg.z);
1476                         if ((function(p->x, p->y, p->z)) > 0.0f)        pos.z = p->z; else neg.z = p->z;
1477                 }
1478         }
1479
1480         /* This is necessary to find start point */
1481         while (1) {
1482                 p->x = 0.5f*(pos.x + neg.x);
1483                 p->y = 0.5f*(pos.y + neg.y);
1484                 p->z = 0.5f*(pos.z + neg.z);
1485
1486                 if (i++ == RES) return;
1487    
1488                 if ((function(p->x, p->y, p->z)) > 0.0f) {
1489                         pos.x = p->x;
1490                         pos.y = p->y;
1491                         pos.z = p->z;
1492                 }
1493                 else {
1494                         neg.x = p->x;
1495                         neg.y = p->y;
1496                         neg.z = p->z;
1497                 }
1498         }
1499 }
1500
1501 /* ************************************** */
1502 static void add_cube(PROCESS *mbproc, int i, int j, int k, int count)
1503 {
1504         CUBES *ncube;
1505         int n;
1506         int a, b, c;
1507
1508         /* hmmm, not only one, but eight cube will be added on the stack 
1509          * ... */
1510         for (a=i-1; a<i+count; a++)
1511                 for (b=j-1; b<j+count; b++)
1512                         for (c=k-1; c<k+count; c++) {
1513                                 /* test if cube has been found before */
1514                                 if ( setcenter(mbproc->centers, a, b, c)==0 ) {
1515                                         /* push cube on stack: */
1516                                         ncube= (CUBES *) new_pgn_element(sizeof(CUBES));
1517                                         ncube->next= mbproc->cubes;
1518                                         mbproc->cubes= ncube;
1519
1520                                         ncube->cube.i= a;
1521                                         ncube->cube.j= b;
1522                                         ncube->cube.k= c;
1523
1524                                         /* set corners of initial cube: */
1525                                         for (n = 0; n < 8; n++)
1526                                         ncube->cube.corners[n] = setcorner(mbproc, a+MB_BIT(n, 2), b+MB_BIT(n, 1), c+MB_BIT(n, 0));
1527                                 }
1528                         }
1529 }
1530
1531
1532 static void find_first_points(PROCESS *mbproc, MetaBall *mb, int a)
1533 {
1534         MB_POINT IN, in, OUT, out; /*point;*/
1535         MetaElem *ml;
1536         int i, j, k, c_i, c_j, c_k;
1537         int index[3]={1, 0, -1};
1538         float f =0.0f;
1539         float in_v /*, out_v*/;
1540         MB_POINT workp;
1541         float tmp_v, workp_v, max_len, len, dx, dy, dz, nx, ny, nz, MAXN;
1542
1543         ml = mainb[a];
1544
1545         f = 1-(mb->thresh/ml->s);
1546
1547         /* Skip, when Stiffness of MetaElement is too small ... MetaElement can't be
1548          * visible alone ... but still can influence others MetaElements :-) */
1549         if (f > 0.0f) {
1550                 OUT.x = IN.x = in.x= 0.0;
1551                 OUT.y = IN.y = in.y= 0.0;
1552                 OUT.z = IN.z = in.z= 0.0;
1553
1554                 calc_mballco(ml, (float *)&in);
1555                 in_v = mbproc->function(in.x, in.y, in.z);
1556
1557                 for (i=0;i<3;i++) {
1558                         switch (ml->type) {
1559                                 case MB_BALL:
1560                                         OUT.x = out.x= IN.x + index[i]*ml->rad;
1561                                         break;
1562                                 case MB_TUBE:
1563                                 case MB_PLANE:
1564                                 case MB_ELIPSOID:
1565                                 case MB_CUBE:
1566                                         OUT.x = out.x= IN.x + index[i]*(ml->expx + ml->rad);
1567                                         break;
1568                         }
1569
1570                         for (j=0;j<3;j++) {
1571                                 switch (ml->type) {
1572                                         case MB_BALL:
1573                                                 OUT.y = out.y= IN.y + index[j]*ml->rad;
1574                                                 break;
1575                                         case MB_TUBE:
1576                                         case MB_PLANE:
1577                                         case MB_ELIPSOID:
1578                                         case MB_CUBE:
1579                                                 OUT.y = out.y= IN.y + index[j]*(ml->expy + ml->rad);
1580                                                 break;
1581                                 }
1582                         
1583                                 for (k=0;k<3;k++) {
1584                                         out.x = OUT.x;
1585                                         out.y = OUT.y;
1586                                         switch (ml->type) {
1587                                                 case MB_BALL:
1588                                                 case MB_TUBE:
1589                                                 case MB_PLANE:
1590                                                         out.z= IN.z + index[k]*ml->rad;
1591                                                         break;
1592                                                 case MB_ELIPSOID:
1593                                                 case MB_CUBE:
1594                                                         out.z= IN.z + index[k]*(ml->expz + ml->rad);
1595                                                         break;
1596                                         }
1597
1598                                         calc_mballco(ml, (float *)&out);
1599
1600                                         /*out_v = mbproc->function(out.x, out.y, out.z);*/ /*UNUSED*/
1601
1602                                         /* find "first points" on Implicit Surface of MetaElemnt ml */
1603                                         workp.x = in.x;
1604                                         workp.y = in.y;
1605                                         workp.z = in.z;
1606                                         workp_v = in_v;
1607                                         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));
1608
1609                                         nx = abs((out.x - in.x)/mbproc->size);
1610                                         ny = abs((out.y - in.y)/mbproc->size);
1611                                         nz = abs((out.z - in.z)/mbproc->size);
1612                                         
1613                                         MAXN = MAX3(nx, ny, nz);
1614                                         if (MAXN!=0.0f) {
1615                                                 dx = (out.x - in.x)/MAXN;
1616                                                 dy = (out.y - in.y)/MAXN;
1617                                                 dz = (out.z - in.z)/MAXN;
1618
1619                                                 len = 0.0;
1620                                                 while (len<=max_len) {
1621                                                         workp.x += dx;
1622                                                         workp.y += dy;
1623                                                         workp.z += dz;
1624                                                         /* compute value of implicite function */
1625                                                         tmp_v = mbproc->function(workp.x, workp.y, workp.z);
1626                                                         /* add cube to the stack, when value of implicite function crosses zero value */
1627                                                         if ((tmp_v<0.0f && workp_v>=0.0f)||(tmp_v>0.0f && workp_v<=0.0f)) {
1628
1629                                                                 /* indexes of CUBE, which includes "first point" */
1630                                                                 c_i= (int)floor(workp.x/mbproc->size);
1631                                                                 c_j= (int)floor(workp.y/mbproc->size);
1632                                                                 c_k= (int)floor(workp.z/mbproc->size);
1633                                                                 
1634                                                                 /* add CUBE (with indexes c_i, c_j, c_k) to the stack,
1635                                                                  * this cube includes found point of Implicit Surface */
1636                                                                 if (ml->flag & MB_NEGATIVE)
1637                                                                         add_cube(mbproc, c_i, c_j, c_k, 2);
1638                                                                 else
1639                                                                         add_cube(mbproc, c_i, c_j, c_k, 1);
1640                                                         }
1641                                                         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));
1642                                                         workp_v = tmp_v;
1643
1644                                                 }
1645                                         }
1646                                 }
1647                         }
1648                 }
1649         }
1650 }
1651
1652 static void polygonize(PROCESS *mbproc, MetaBall *mb)
1653 {
1654         CUBE c;
1655         int a;
1656
1657         mbproc->vertices.count = mbproc->vertices.max = 0;
1658         mbproc->vertices.ptr = NULL;
1659
1660         /* allocate hash tables and build cube polygon table: */
1661         mbproc->centers = MEM_callocN(HASHSIZE * sizeof(CENTERLIST *), "mbproc->centers");
1662         mbproc->corners = MEM_callocN(HASHSIZE * sizeof(CORNER *), "mbproc->corners");
1663         mbproc->edges = MEM_callocN(2*HASHSIZE * sizeof(EDGELIST *), "mbproc->edges");
1664         makecubetable();
1665
1666         for (a=0; a<totelem; a++) {
1667
1668                 /* try to find 8 points on the surface for each MetaElem */
1669                 find_first_points(mbproc, mb, a);       
1670         }
1671
1672         /* polygonize all MetaElems of current MetaBall */
1673         while (mbproc->cubes != NULL) { /* process active cubes till none left */
1674                 c = mbproc->cubes->cube;
1675
1676                 /* polygonize the cube directly: */
1677                 docube(&c, mbproc, mb);
1678                 
1679                 /* pop current cube from stack */
1680                 mbproc->cubes = mbproc->cubes->next;
1681                 
1682                 /* test six face directions, maybe add to stack: */
1683                 testface(c.i-1, c.j, c.k, &c, 2, LBN, LBF, LTN, LTF, mbproc);
1684                 testface(c.i+1, c.j, c.k, &c, 2, RBN, RBF, RTN, RTF, mbproc);
1685                 testface(c.i, c.j-1, c.k, &c, 1, LBN, LBF, RBN, RBF, mbproc);
1686                 testface(c.i, c.j+1, c.k, &c, 1, LTN, LTF, RTN, RTF, mbproc);
1687                 testface(c.i, c.j, c.k-1, &c, 0, LBN, LTN, RBN, RTN, mbproc);
1688                 testface(c.i, c.j, c.k+1, &c, 0, LBF, LTF, RBF, RTF, mbproc);
1689         }
1690 }
1691
1692 static float init_meta(Scene *scene, Object *ob)        /* return totsize */
1693 {
1694         Scene *sce_iter= scene;
1695         Base *base;
1696         Object *bob;
1697         MetaBall *mb;
1698         MetaElem *ml;
1699         float size, totsize, obinv[4][4], obmat[4][4], vec[3];
1700         //float max=0.0;
1701         int a, obnr, zero_size=0;
1702         char obname[MAX_ID_NAME];
1703         
1704         copy_m4_m4(obmat, ob->obmat);   /* to cope with duplicators from next_object */
1705         invert_m4_m4(obinv, ob->obmat);
1706         a= 0;
1707         
1708         BLI_split_name_num(obname, &obnr, ob->id.name+2, '.');
1709         
1710         /* make main array */
1711         next_object(&sce_iter, 0, NULL, NULL);
1712         while (next_object(&sce_iter, 1, &base, &bob)) {
1713
1714                 if (bob->type==OB_MBALL) {
1715                         zero_size= 0;
1716                         ml= NULL;
1717
1718                         if (bob==ob && (base->flag & OB_FROMDUPLI)==0) {
1719                                 mb= ob->data;
1720         
1721                                 if (mb->editelems) ml= mb->editelems->first;
1722                                 else ml= mb->elems.first;
1723                         }
1724                         else {
1725                                 char name[MAX_ID_NAME];
1726                                 int nr;
1727                                 
1728                                 BLI_split_name_num(name, &nr, bob->id.name+2, '.');
1729                                 if ( strcmp(obname, name)==0 ) {
1730                                         mb= bob->data;
1731                                         
1732                                         if (mb->editelems) ml= mb->editelems->first;
1733                                         else ml= mb->elems.first;
1734                                 }
1735                         }
1736
1737                         /* when metaball object has zero scale, then MetaElem to this MetaBall
1738                          * will not be put to mainb array */
1739                         if (bob->size[0]==0.0f || bob->size[1]==0.0f || bob->size[2]==0.0f) {
1740                                 zero_size= 1;
1741                         }
1742                         else if (bob->parent) {
1743                                 struct Object *pob=bob->parent;
1744                                 while (pob) {
1745                                         if (pob->size[0]==0.0f || pob->size[1]==0.0f || pob->size[2]==0.0f) {
1746                                                 zero_size= 1;
1747                                                 break;
1748                                         }
1749                                         pob= pob->parent;
1750                                 }
1751                         }
1752
1753                         if (zero_size) {
1754                                 unsigned int ml_count=0;
1755                                 while (ml) {
1756                                         ml_count++;
1757                                         ml= ml->next;
1758                                 }
1759                                 totelem -= ml_count;
1760                         }
1761                         else {
1762                         while (ml) {
1763                                 if (!(ml->flag & MB_HIDE)) {
1764                                         int i;
1765                                         float temp1[4][4], temp2[4][4], temp3[4][4];
1766                                         float (*mat)[4] = NULL, (*imat)[4] = NULL;
1767                                         float max_x, max_y, max_z, min_x, min_y, min_z;
1768
1769                                         max_x = max_y = max_z = -3.4e38;
1770                                         min_x = min_y = min_z =  3.4e38;
1771
1772                                         /* too big stiffness seems only ugly due to linear interpolation
1773                                          * no need to have possibility for too big stiffness */
1774                                         if (ml->s > 10.0f) ml->s = 10.0f;
1775                                         
1776                                         /* Rotation of MetaElem is stored in quat */
1777                                          quat_to_mat4(temp3, ml->quat);
1778
1779                                         /* Translation of MetaElem */
1780                                         unit_m4(temp2);
1781                                         temp2[3][0]= ml->x;
1782                                         temp2[3][1]= ml->y;
1783                                         temp2[3][2]= ml->z;
1784
1785                                         mult_m4_m4m4(temp1, temp2, temp3);
1786                                 
1787                                         /* make a copy because of duplicates */
1788                                         mainb[a]= new_pgn_element(sizeof(MetaElem));
1789                                         *(mainb[a])= *ml;
1790                                         mainb[a]->bb = new_pgn_element(sizeof(BoundBox));
1791                                 
1792                                         mat= new_pgn_element(4*4*sizeof(float));
1793                                         imat= new_pgn_element(4*4*sizeof(float));
1794                                         
1795                                         /* mat is the matrix to transform from mball into the basis-mball */
1796                                         invert_m4_m4(obinv, obmat);
1797                                         mult_m4_m4m4(temp2, obinv, bob->obmat);
1798                                         /* MetaBall transformation */
1799                                         mult_m4_m4m4(mat, temp2, temp1);
1800
1801                                         invert_m4_m4(imat, mat);                                
1802
1803                                         mainb[a]->rad2= ml->rad*ml->rad;
1804
1805                                         mainb[a]->mat= (float*) mat;
1806                                         mainb[a]->imat= (float*) imat;
1807
1808                                         /* untransformed Bounding Box of MetaElem */
1809                                         /* 0 */
1810                                         mainb[a]->bb->vec[0][0]= -ml->expx;
1811                                         mainb[a]->bb->vec[0][1]= -ml->expy;
1812                                         mainb[a]->bb->vec[0][2]= -ml->expz;
1813                                         /* 1 */
1814                                         mainb[a]->bb->vec[1][0]=  ml->expx;
1815                                         mainb[a]->bb->vec[1][1]= -ml->expy;
1816                                         mainb[a]->bb->vec[1][2]= -ml->expz;
1817                                         /* 2 */
1818                                         mainb[a]->bb->vec[2][0]=  ml->expx;
1819                                         mainb[a]->bb->vec[2][1]=  ml->expy;
1820                                         mainb[a]->bb->vec[2][2]= -ml->expz;
1821                                         /* 3 */
1822                                         mainb[a]->bb->vec[3][0]= -ml->expx;
1823                                         mainb[a]->bb->vec[3][1]=  ml->expy;
1824                                         mainb[a]->bb->vec[3][2]= -ml->expz;
1825                                         /* 4 */
1826                                         mainb[a]->bb->vec[4][0]= -ml->expx;
1827                                         mainb[a]->bb->vec[4][1]= -ml->expy;
1828                                         mainb[a]->bb->vec[4][2]=  ml->expz;
1829                                         /* 5 */
1830                                         mainb[a]->bb->vec[5][0]=  ml->expx;
1831                                         mainb[a]->bb->vec[5][1]= -ml->expy;
1832                                         mainb[a]->bb->vec[5][2]=  ml->expz;
1833                                         /* 6 */
1834                                         mainb[a]->bb->vec[6][0]=  ml->expx;
1835                                         mainb[a]->bb->vec[6][1]=  ml->expy;
1836                                         mainb[a]->bb->vec[6][2]=  ml->expz;
1837                                         /* 7 */
1838                                         mainb[a]->bb->vec[7][0]= -ml->expx;
1839                                         mainb[a]->bb->vec[7][1]=  ml->expy;
1840                                         mainb[a]->bb->vec[7][2]=  ml->expz;
1841
1842                                         /* transformation of Metalem bb */
1843                                         for (i=0; i<8; i++)
1844                                                 mul_m4_v3((float (*)[4])mat, mainb[a]->bb->vec[i]);
1845
1846                                         /* find max and min of transformed bb */
1847                                         for (i=0; i<8; i++) {
1848                                                 /* find maximums */
1849                                                 if (mainb[a]->bb->vec[i][0] > max_x) max_x = mainb[a]->bb->vec[i][0];
1850                                                 if (mainb[a]->bb->vec[i][1] > max_y) max_y = mainb[a]->bb->vec[i][1];
1851                                                 if (mainb[a]->bb->vec[i][2] > max_z) max_z = mainb[a]->bb->vec[i][2];
1852                                                 /* find  minimums */
1853                                                 if (mainb[a]->bb->vec[i][0] < min_x) min_x = mainb[a]->bb->vec[i][0];
1854                                                 if (mainb[a]->bb->vec[i][1] < min_y) min_y = mainb[a]->bb->vec[i][1];
1855                                                 if (mainb[a]->bb->vec[i][2] < min_z) min_z = mainb[a]->bb->vec[i][2];
1856                                         }
1857
1858                                         /* create "new" bb, only point 0 and 6, which are
1859                                          * necessary for octal tree filling */
1860                                         mainb[a]->bb->vec[0][0] = min_x - ml->rad;
1861                                         mainb[a]->bb->vec[0][1] = min_y - ml->rad;
1862                                         mainb[a]->bb->vec[0][2] = min_z - ml->rad;
1863
1864                                         mainb[a]->bb->vec[6][0] = max_x + ml->rad;
1865                                         mainb[a]->bb->vec[6][1] = max_y + ml->rad;
1866                                         mainb[a]->bb->vec[6][2] = max_z + ml->rad;
1867                                         
1868                                         a++;
1869                                 }
1870                                 ml= ml->next;
1871                         }
1872                         }
1873                 }
1874         }
1875
1876         
1877         /* totsize (= 'manhattan' radius) */
1878         totsize= 0.0;
1879         for (a=0; a<totelem; a++) {
1880                 
1881                 vec[0]= mainb[a]->x + mainb[a]->rad + mainb[a]->expx;
1882                 vec[1]= mainb[a]->y + mainb[a]->rad + mainb[a]->expy;
1883                 vec[2]= mainb[a]->z + mainb[a]->rad + mainb[a]->expz;
1884
1885                 calc_mballco(mainb[a], vec);
1886         
1887                 size= fabsf(vec[0]);
1888                 if ( size > totsize ) totsize= size;
1889                 size= fabsf(vec[1]);
1890                 if ( size > totsize ) totsize= size;
1891                 size= fabsf(vec[2]);
1892                 if ( size > totsize ) totsize= size;
1893
1894                 vec[0]= mainb[a]->x - mainb[a]->rad;
1895                 vec[1]= mainb[a]->y - mainb[a]->rad;
1896                 vec[2]= mainb[a]->z - mainb[a]->rad;
1897                                 
1898                 calc_mballco(mainb[a], vec);
1899         
1900                 size= fabsf(vec[0]);
1901                 if ( size > totsize ) totsize= size;
1902                 size= fabsf(vec[1]);
1903                 if ( size > totsize ) totsize= size;
1904                 size= fabsf(vec[2]);
1905                 if ( size > totsize ) totsize= size;
1906         }
1907
1908         for (a=0; a<totelem; a++) {
1909                 thresh += densfunc(mainb[a], 2.0f * totsize, 2.0f * totsize, 2.0f * totsize);
1910         }
1911
1912         return totsize;
1913 }
1914
1915 /* if MetaElem lies in node, then node includes MetaElem pointer (ml_p)
1916  * pointing at MetaElem (ml)
1917  */
1918 static void fill_metaball_octal_node(octal_node *node, MetaElem *ml, short i)
1919 {
1920         ml_pointer *ml_p;
1921
1922         ml_p= MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
1923         ml_p->ml= ml;
1924         BLI_addtail(&(node->nodes[i]->elems), ml_p);
1925         node->count++;
1926         
1927         if (ml->flag & MB_NEGATIVE) {
1928                 node->nodes[i]->neg++;
1929         }
1930         else {
1931                 node->nodes[i]->pos++;
1932         }
1933 }
1934
1935 /* Node is subdivided as is illustrated on the following figure:
1936  * 
1937  *      +------+------+
1938  *     /      /      /|
1939  *    +------+------+ |
1940  *   /      /      /| +
1941  *  +------+------+ |/|
1942  *  |      |      | + |
1943  *  |      |      |/| +
1944  *  +------+------+ |/
1945  *  |      |      | +
1946  *  |      |      |/
1947  *  +------+------+
1948  *  
1949  */
1950 static void subdivide_metaball_octal_node(octal_node *node, float size_x, float size_y, float size_z, short depth)
1951 {
1952         MetaElem *ml;
1953         ml_pointer *ml_p;
1954         float x, y, z;
1955         int a, i;
1956
1957         /* create new nodes */
1958         for (a=0;a<8;a++) {
1959                 node->nodes[a]= MEM_mallocN(sizeof(octal_node), "octal_node");
1960                 for (i=0;i<8;i++)
1961                         node->nodes[a]->nodes[i]= NULL;
1962                 node->nodes[a]->parent= node;
1963                 node->nodes[a]->elems.first= NULL;
1964                 node->nodes[a]->elems.last= NULL;
1965                 node->nodes[a]->count= 0;
1966                 node->nodes[a]->neg= 0;
1967                 node->nodes[a]->pos= 0;
1968         }
1969
1970         size_x /= 2;
1971         size_y /= 2;
1972         size_z /= 2;
1973         
1974         /* center of node */
1975         node->x = x = node->x_min + size_x;
1976         node->y = y = node->y_min + size_y;
1977         node->z = z = node->z_min + size_z;
1978
1979         /* setting up of border points of new nodes */
1980         node->nodes[0]->x_min = node->x_min;
1981         node->nodes[0]->y_min = node->y_min;
1982         node->nodes[0]->z_min = node->z_min;
1983         node->nodes[0]->x = node->nodes[0]->x_min + size_x/2;
1984         node->nodes[0]->y = node->nodes[0]->y_min + size_y/2;
1985         node->nodes[0]->z = node->nodes[0]->z_min + size_z/2;
1986         
1987         node->nodes[1]->x_min = x;
1988         node->nodes[1]->y_min = node->y_min;
1989         node->nodes[1]->z_min = node->z_min;
1990         node->nodes[1]->x = node->nodes[1]->x_min + size_x/2;
1991         node->nodes[1]->y = node->nodes[1]->y_min + size_y/2;
1992         node->nodes[1]->z = node->nodes[1]->z_min + size_z/2;
1993
1994         node->nodes[2]->x_min = x;
1995         node->nodes[2]->y_min = y;
1996         node->nodes[2]->z_min = node->z_min;
1997         node->nodes[2]->x = node->nodes[2]->x_min + size_x/2;
1998         node->nodes[2]->y = node->nodes[2]->y_min + size_y/2;
1999         node->nodes[2]->z = node->nodes[2]->z_min + size_z/2;
2000
2001         node->nodes[3]->x_min = node->x_min;
2002         node->nodes[3]->y_min = y;
2003         node->nodes[3]->z_min = node->z_min;
2004         node->nodes[3]->x = node->nodes[3]->x_min + size_x/2;
2005         node->nodes[3]->y = node->nodes[3]->y_min + size_y/2;
2006         node->nodes[3]->z = node->nodes[3]->z_min + size_z/2;
2007
2008         node->nodes[4]->x_min = node->x_min;
2009         node->nodes[4]->y_min = node->y_min;
2010         node->nodes[4]->z_min = z;
2011         node->nodes[4]->x = node->nodes[4]->x_min + size_x/2;
2012         node->nodes[4]->y = node->nodes[4]->y_min + size_y/2;
2013         node->nodes[4]->z = node->nodes[4]->z_min + size_z/2;
2014         
2015         node->nodes[5]->x_min = x;
2016         node->nodes[5]->y_min = node->y_min;
2017         node->nodes[5]->z_min = z;
2018         node->nodes[5]->x = node->nodes[5]->x_min + size_x/2;
2019         node->nodes[5]->y = node->nodes[5]->y_min + size_y/2;
2020         node->nodes[5]->z = node->nodes[5]->z_min + size_z/2;
2021
2022         node->nodes[6]->x_min = x;
2023         node->nodes[6]->y_min = y;
2024         node->nodes[6]->z_min = z;
2025         node->nodes[6]->x = node->nodes[6]->x_min + size_x/2;
2026         node->nodes[6]->y = node->nodes[6]->y_min + size_y/2;
2027         node->nodes[6]->z = node->nodes[6]->z_min + size_z/2;
2028
2029         node->nodes[7]->x_min = node->x_min;
2030         node->nodes[7]->y_min = y;
2031         node->nodes[7]->z_min = z;
2032         node->nodes[7]->x = node->nodes[7]->x_min + size_x/2;
2033         node->nodes[7]->y = node->nodes[7]->y_min + size_y/2;
2034         node->nodes[7]->z = node->nodes[7]->z_min + size_z/2;
2035
2036         ml_p= node->elems.first;
2037         
2038         /* setting up references of MetaElems for new nodes */
2039         while (ml_p) {
2040                 ml= ml_p->ml;
2041                 if (ml->bb->vec[0][2] < z) {
2042                         if (ml->bb->vec[0][1] < y) {
2043                                 /* vec[0][0] lies in first octant */
2044                                 if (ml->bb->vec[0][0] < x) {
2045                                         /* ml belongs to the (0)1st node */
2046                                         fill_metaball_octal_node(node, ml, 0);
2047
2048                                         /* ml belongs to the (3)4th node */
2049                                         if (ml->bb->vec[6][1] >= y) {
2050                                                 fill_metaball_octal_node(node, ml, 3);
2051
2052                                                 /* ml belongs to the (7)8th node */
2053                                                 if (ml->bb->vec[6][2] >= z) {
2054                                                         fill_metaball_octal_node(node, ml, 7);
2055                                                 }
2056                                         }
2057         
2058                                         /* ml belongs to the (1)2nd node */
2059                                         if (ml->bb->vec[6][0] >= x) {
2060                                                 fill_metaball_octal_node(node, ml, 1);
2061
2062                                                 /* ml belongs to the (5)6th node */
2063                                                 if (ml->bb->vec[6][2] >= z) {
2064                                                         fill_metaball_octal_node(node, ml, 5);
2065                                                 }
2066                                         }
2067
2068                                         /* ml belongs to the (2)3th node */
2069                                         if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
2070                                                 fill_metaball_octal_node(node, ml, 2);
2071                                                 
2072                                                 /* ml belong to the (6)7th node */
2073                                                 if (ml->bb->vec[6][2] >= z) {
2074                                                         fill_metaball_octal_node(node, ml, 6);
2075                                                 }
2076                                                 
2077                                         }
2078                         
2079                                         /* ml belongs to the (4)5th node too */ 
2080                                         if (ml->bb->vec[6][2] >= z) {
2081                                                 fill_metaball_octal_node(node, ml, 4);
2082                                         }
2083
2084                                         
2085                                         
2086                                 }
2087                                 /* vec[0][0] is in the (1)second octant */
2088                                 else {
2089                                         /* ml belong to the (1)2nd node */
2090                                         fill_metaball_octal_node(node, ml, 1);
2091
2092                                         /* ml belongs to the (2)3th node */
2093                                         if (ml->bb->vec[6][1] >= y) {
2094                                                 fill_metaball_octal_node(node, ml, 2);
2095
2096                                                 /* ml belongs to the (6)7th node */
2097                                                 if (ml->bb->vec[6][2] >= z) {
2098                                                         fill_metaball_octal_node(node, ml, 6);
2099                                                 }
2100                                                 
2101                                         }
2102                                         
2103                                         /* ml belongs to the (5)6th node */
2104                                         if (ml->bb->vec[6][2] >= z) {
2105                                                 fill_metaball_octal_node(node, ml, 5);
2106                                         }
2107                                 }
2108                         }
2109                         else {
2110                                 /* vec[0][0] is in the (3)4th octant */
2111                                 if (ml->bb->vec[0][0] < x) {
2112                                         /* ml belongs to the (3)4nd node */
2113                                         fill_metaball_octal_node(node, ml, 3);
2114                                         
2115                                         /* ml belongs to the (7)8th node */
2116                                         if (ml->bb->vec[6][2] >= z) {
2117                                                 fill_metaball_octal_node(node, ml, 7);
2118                                         }
2119                                 
2120
2121                                         /* ml belongs to the (2)3th node */
2122                                         if (ml->bb->vec[6][0] >= x) {
2123                                                 fill_metaball_octal_node(node, ml, 2);
2124                                         
2125                                                 /* ml belongs to the (6)7th node */
2126                                                 if (ml->bb->vec[6][2] >= z) {
2127                                                         fill_metaball_octal_node(node, ml, 6);
2128                                                 }
2129                                         }
2130                                 }
2131
2132                         }
2133
2134                         /* vec[0][0] is in the (2)3th octant */
2135                         if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
2136                                 /* ml belongs to the (2)3th node */
2137                                 fill_metaball_octal_node(node, ml, 2);
2138                                 
2139                                 /* ml belongs to the (6)7th node */
2140                                 if (ml->bb->vec[6][2] >= z) {
2141                                         fill_metaball_octal_node(node, ml, 6);
2142                                 }
2143                         }
2144                 }
2145                 else {
2146                         if (ml->bb->vec[0][1] < y) {
2147                                 /* vec[0][0] lies in (4)5th octant */
2148                                 if (ml->bb->vec[0][0] < x) {
2149                                         /* ml belongs to the (4)5th node */
2150                                         fill_metaball_octal_node(node, ml, 4);
2151
2152                                         if (ml->bb->vec[6][0] >= x) {
2153                                                 fill_metaball_octal_node(node, ml, 5);
2154                                         }
2155
2156                                         if (ml->bb->vec[6][1] >= y) {
2157                                                 fill_metaball_octal_node(node, ml, 7);
2158                                         }
2159                                         
2160                                         if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
2161                                                 fill_metaball_octal_node(node, ml, 6);
2162                                         }
2163                                 }
2164                                 /* vec[0][0] lies in (5)6th octant */
2165                                 else {
2166                                         fill_metaball_octal_node(node, ml, 5);
2167
2168                                         if (ml->bb->vec[6][1] >= y) {
2169                                                 fill_metaball_octal_node(node, ml, 6);
2170                                         }
2171                                 }
2172                         }
2173                         else {
2174                                 /* vec[0][0] lies in (7)8th octant */
2175                                 if (ml->bb->vec[0][0] < x) {
2176                                         fill_metaball_octal_node(node, ml, 7);
2177
2178                                         if (ml->bb->vec[6][0] >= x) {
2179                                                 fill_metaball_octal_node(node, ml, 6);
2180                                         }
2181                                 }
2182
2183                         }
2184                         
2185                         /* vec[0][0] lies in (6)7th octant */
2186                         if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
2187                                 fill_metaball_octal_node(node, ml, 6);
2188                         }
2189                 }
2190                 ml_p= ml_p->next;
2191         }
2192
2193         /* free references of MetaElems for curent node (it is not needed anymore) */
2194         BLI_freelistN(&node->elems);
2195
2196         depth--;
2197         
2198         if (depth>0) {
2199                 for (a=0;a<8;a++) {
2200                         if (node->nodes[a]->count > 0) /* if node is not empty, then it is subdivided */
2201                                 subdivide_metaball_octal_node(node->nodes[a], size_x, size_y, size_z, depth);
2202                 }
2203         }
2204 }
2205
2206 /* free all octal nodes recursively */
2207 static void free_metaball_octal_node(octal_node *node)
2208 {
2209         int a;
2210         for (a=0;a<8;a++) {
2211                 if (node->nodes[a]!=NULL) free_metaball_octal_node(node->nodes[a]);
2212         }
2213         BLI_freelistN(&node->elems);
2214         MEM_freeN(node);
2215 }
2216
2217 /* If scene include more then one MetaElem, then octree is used */
2218 static void init_metaball_octal_tree(int depth)
2219 {
2220         struct octal_node *node;
2221         ml_pointer *ml_p;
2222         float size[3];
2223         int a;
2224         
2225         metaball_tree= MEM_mallocN(sizeof(octal_tree), "metaball_octal_tree");
2226         metaball_tree->first= node= MEM_mallocN(sizeof(octal_node), "metaball_octal_node");
2227         /* maximal depth of octree */
2228         metaball_tree->depth= depth;
2229
2230         metaball_tree->neg= node->neg=0;
2231         metaball_tree->pos= node->pos=0;
2232         
2233         node->elems.first= NULL;
2234         node->elems.last= NULL;
2235         node->count=0;
2236
2237         for (a=0;a<8;a++)
2238                 node->nodes[a]=NULL;
2239
2240         node->x_min= node->y_min= node->z_min= FLT_MAX;
2241         node->x_max= node->y_max= node->z_max= -FLT_MAX;
2242
2243         /* size of octal tree scene */
2244         for (a=0;a<totelem;a++) {
2245                 if (mainb[a]->bb->vec[0][0] < node->x_min) node->x_min= mainb[a]->bb->vec[0][0];
2246                 if (mainb[a]->bb->vec[0][1] < node->y_min) node->y_min= mainb[a]->bb->vec[0][1];
2247                 if (mainb[a]->bb->vec[0][2] < node->z_min) node->z_min= mainb[a]->bb->vec[0][2];
2248                 
2249                 if (mainb[a]->bb->vec[6][0] > node->x_max) node->x_max= mainb[a]->bb->vec[6][0];
2250                 if (mainb[a]->bb->vec[6][1] > node->y_max) node->y_max= mainb[a]->bb->vec[6][1];
2251                 if (mainb[a]->bb->vec[6][2] > node->z_max) node->z_max= mainb[a]->bb->vec[6][2];
2252
2253                 ml_p= MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
2254                 ml_p->ml= mainb[a];
2255                 BLI_addtail(&node->elems, ml_p);
2256
2257                 if (mainb[a]->flag & MB_NEGATIVE) {
2258                         /* number of negative MetaElem in scene */
2259                         metaball_tree->neg++;
2260                 }
2261                 else {
2262                         /* number of positive MetaElem in scene */
2263                         metaball_tree->pos++;
2264                 }
2265         }
2266
2267         /* size of first node */        
2268         size[0]= node->x_max - node->x_min;
2269         size[1]= node->y_max - node->y_min;
2270         size[2]= node->z_max - node->z_min;
2271
2272         /* first node is subdivided recursively */
2273         subdivide_metaball_octal_node(node, size[0], size[1], size[2], metaball_tree->depth);
2274 }
2275
2276 void BKE_metaball_polygonize(Scene *scene, Object *ob, ListBase *dispbase)
2277 {
2278         PROCESS mbproc;
2279         MetaBall *mb;
2280         DispList *dl;
2281         int a, nr_cubes;
2282         float *ve, *no, totsize, width;
2283
2284         mb= ob->data;
2285
2286         if (totelem==0) return;
2287         if (!(G.rendering) && (mb->flag==MB_UPDATE_NEVER)) return;
2288         if (G.moving && mb->flag==MB_UPDATE_FAST) return;
2289
2290         curindex= totindex= 0;
2291         indices= NULL;
2292         thresh= mb->thresh;
2293
2294         /* total number of MetaElems (totelem) is precomputed in find_basis_mball() function */
2295         mainb= MEM_mallocN(sizeof(void *)*totelem, "mainb");
2296         
2297         /* initialize all mainb (MetaElems) */
2298         totsize= init_meta(scene, ob);
2299
2300         if (metaball_tree) {
2301                 free_metaball_octal_node(metaball_tree->first);
2302                 MEM_freeN(metaball_tree);
2303                 metaball_tree= NULL;
2304         }
2305
2306         /* if scene includes more then one MetaElem, then octal tree optimalisation is used */  
2307         if ((totelem > 1) && (totelem <= 64)) init_metaball_octal_tree(1);
2308         if ((totelem > 64) && (totelem <= 128)) init_metaball_octal_tree(2);
2309         if ((totelem > 128) && (totelem <= 512))        init_metaball_octal_tree(3);
2310         if ((totelem > 512) && (totelem <= 1024)) init_metaball_octal_tree(4);
2311         if (totelem > 1024) init_metaball_octal_tree(5);
2312
2313         /* don't polygonize metaballs with too high resolution (base mball to small)
2314          * note: Eps was 0.0001f but this was giving problems for blood animation for durian, using 0.00001f */
2315         if (metaball_tree) {
2316                 if (    ob->size[0] <= 0.00001f * (metaball_tree->first->x_max - metaball_tree->first->x_min) ||
2317                         ob->size[1] <= 0.00001f * (metaball_tree->first->y_max - metaball_tree->first->y_min) ||
2318                         ob->size[2] <= 0.00001f * (metaball_tree->first->z_max - metaball_tree->first->z_min))
2319                 {
2320                         new_pgn_element(-1); /* free values created by init_meta */
2321
2322                         MEM_freeN(mainb);
2323
2324                         /* free tree */
2325                         free_metaball_octal_node(metaball_tree->first);
2326                         MEM_freeN(metaball_tree);
2327                         metaball_tree= NULL;
2328
2329                         return;
2330                 }
2331         }
2332
2333         /* width is size per polygonize cube */
2334         if (G.rendering) width= mb->rendersize;
2335         else {
2336                 width= mb->wiresize;
2337                 if (G.moving && mb->flag==MB_UPDATE_HALFRES) width*= 2;
2338         }
2339         /* nr_cubes is just for safety, minimum is totsize */
2340         nr_cubes= (int)(0.5f+totsize/width);
2341
2342         /* init process */
2343         mbproc.function = metaball;
2344         mbproc.size = width;
2345         mbproc.bounds = nr_cubes;
2346         mbproc.cubes= NULL;
2347         mbproc.delta = width/(float)(RES*RES);
2348
2349         polygonize(&mbproc, mb);
2350         
2351         MEM_freeN(mainb);
2352
2353         /* free octal tree */
2354         if (totelem > 1) {
2355                 free_metaball_octal_node(metaball_tree->first);
2356                 MEM_freeN(metaball_tree);
2357                 metaball_tree= NULL;
2358         }
2359
2360         if (curindex) {
2361                 dl= MEM_callocN(sizeof(DispList), "mbaldisp");
2362                 BLI_addtail(dispbase, dl);
2363                 dl->type= DL_INDEX4;
2364                 dl->nr= mbproc.vertices.count;
2365                 dl->parts= curindex;
2366
2367                 dl->index= indices;
2368                 indices= NULL;
2369                 
2370                 a= mbproc.vertices.count;
2371                 dl->verts= ve= MEM_mallocN(sizeof(float)*3*a, "mballverts");
2372                 dl->nors= no= MEM_mallocN(sizeof(float)*3*a, "mballnors");
2373
2374                 for (a=0; a<mbproc.vertices.count; a++, no+=3, ve+=3) {
2375                         ve[0]= mbproc.vertices.ptr[a].position.x;
2376                         ve[1]= mbproc.vertices.ptr[a].position.y;
2377                         ve[2]= mbproc.vertices.ptr[a].position.z;
2378
2379                         no[0]= mbproc.vertices.ptr[a].normal.x;
2380                         no[1]= mbproc.vertices.ptr[a].normal.y;
2381                         no[2]= mbproc.vertices.ptr[a].normal.z;
2382                 }
2383         }
2384
2385         freepolygonize(&mbproc);
2386 }
2387
2388 /* basic vertex data functions */
2389 int BKE_metaball_minmax(MetaBall *mb, float min[3], float max[3])
2390 {
2391         MetaElem *ml;
2392
2393         INIT_MINMAX(min, max);
2394
2395         for (ml = mb->elems.first; ml; ml = ml->next) {
2396                 DO_MINMAX(&ml->x, min, max);
2397         }
2398
2399         return (mb->elems.first != NULL);
2400 }
2401
2402 int BKE_metaball_center_median(MetaBall *mb, float cent[3])
2403 {
2404         MetaElem *ml;
2405         int total= 0;
2406
2407         zero_v3(cent);
2408
2409         for (ml = mb->elems.first; ml; ml = ml->next) {
2410                 add_v3_v3(cent, &ml->x);
2411         }
2412
2413         if (total)
2414                 mul_v3_fl(cent, 1.0f/(float)total);
2415
2416         return (total != 0);
2417 }
2418
2419 int BKE_metaball_center_bounds(MetaBall *mb, float cent[3])
2420 {
2421         float min[3], max[3];
2422
2423         if (BKE_metaball_minmax(mb, min, max)) {
2424                 mid_v3_v3v3(cent, min, max);
2425                 return 1;
2426         }
2427
2428         return 0;
2429 }
2430
2431 void BKE_metaball_translate(MetaBall *mb, float offset[3])
2432 {
2433         MetaElem *ml;
2434
2435         for (ml = mb->elems.first; ml; ml = ml->next) {
2436                 add_v3_v3(&ml->x, offset);
2437         }
2438 }