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