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