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