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