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