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