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