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