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