2 * ***** BEGIN GPL LICENSE BLOCK *****
4 * This program is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU General Public License
6 * as published by the Free Software Foundation; either version 2
7 * of the License, or (at your option) any later version.
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the Free Software Foundation,
16 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
19 * All rights reserved.
21 * Contributor(s): Jiri Hnidek <jiri.hnidek@vslib.cz>.
23 * ***** END GPL LICENSE BLOCK *****
25 * MetaBalls are created from a single Object (with a name without number in it),
26 * here the DispList and BoundBox also is located.
27 * All objects with the same name (but with a number in it) are added to this.
29 * texture coordinates are patched within the displist
32 /** \file blender/blenkernel/intern/mball.c
42 #include "MEM_guardedalloc.h"
44 #include "DNA_material_types.h"
45 #include "DNA_object_types.h"
46 #include "DNA_meta_types.h"
47 #include "DNA_scene_types.h"
50 #include "BLI_blenlib.h"
52 #include "BLI_utildefines.h"
53 #include "BLI_bpath.h"
56 #include "BKE_global.h"
59 /* #include "BKE_object.h" */
60 #include "BKE_animsys.h"
61 #include "BKE_scene.h"
62 #include "BKE_library.h"
63 #include "BKE_displist.h"
64 #include "BKE_mball.h"
65 #include "BKE_object.h"
66 #include "BKE_material.h"
70 typedef struct vertex { /* surface vertex */
71 float co[3]; /* position and surface normal */
75 typedef struct vertices { /* list of vertices in polygonization */
76 int count, max; /* # vertices, max # allowed */
77 VERTEX *ptr; /* dynamically allocated */
80 typedef struct corner { /* corner of a cube */
81 int i, j, k; /* (i, j, k) is index within lattice */
82 float co[3], value; /* location and function value */
86 typedef struct cube { /* partitioning cell (cube) */
87 int i, j, k; /* lattice location of cube */
88 CORNER *corners[8]; /* eight corners */
91 typedef struct cubes { /* linked list of cubes acting as stack */
92 CUBE cube; /* a single cube */
93 struct cubes *next; /* remaining elements */
96 typedef struct centerlist { /* list of cube locations */
97 int i, j, k; /* cube location */
98 struct centerlist *next; /* remaining elements */
101 typedef struct edgelist { /* list of edges */
102 int i1, j1, k1, i2, j2, k2; /* edge corner ids */
103 int vid; /* vertex id */
104 struct edgelist *next; /* remaining elements */
107 typedef struct intlist { /* list of integers */
108 int i; /* an integer */
109 struct intlist *next; /* remaining elements */
112 typedef struct intlists { /* list of list of integers */
113 INTLIST *list; /* a list of integers */
114 struct intlists *next; /* remaining elements */
117 typedef struct process { /* parameters, function, storage */
118 /* what happens here? floats, I think. */
119 /* float (*function)(void); */ /* implicit surface function */
120 float (*function)(float, float, float);
121 float size, delta; /* cube size, normal delta */
122 int bounds; /* cube range within lattice */
123 CUBES *cubes; /* active cubes */
124 VERTICES vertices; /* surface vertices */
125 CENTERLIST **centers; /* cube center hash table */
126 CORNER **corners; /* corner value hash table */
127 EDGELIST **edges; /* edge and vertex id hash table */
130 /* dividing scene using octal tree makes polygonisation faster */
131 typedef struct ml_pointer {
132 struct ml_pointer *next, *prev;
136 typedef struct octal_node {
137 struct octal_node *nodes[8];/* children of current node */
138 struct octal_node *parent; /* parent of current node */
139 struct ListBase elems; /* ListBase of MetaElem pointers (ml_pointer) */
140 float x_min, y_min, z_min; /* 1st border point */
141 float x_max, y_max, z_max; /* 7th border point */
142 float x, y, z; /* center of node */
143 int pos, neg; /* number of positive and negative MetaElements in the node */
144 int count; /* number of MetaElems, which belongs to the node */
147 typedef struct octal_tree {
148 struct octal_node *first; /* first node */
149 int pos, neg; /* number of positive and negative MetaElements in the scene */
150 short depth; /* number of scene subdivision */
153 struct pgn_elements {
154 struct pgn_elements *next, *prev;
158 /* Forward declarations */
159 static int vertid(const CORNER *c1, const CORNER *c2, PROCESS *p, MetaBall *mb);
160 static int setcenter(CENTERLIST *table[], const int i, const int j, const int k);
161 static CORNER *setcorner(PROCESS *p, int i, int j, int k);
162 static void converge(const float p1[3], const float p2[3], float v1, float v2,
163 float (*function)(float, float, float), float p[3], MetaBall *mb, int f);
165 /* Global variables */
170 octal_tree *metaball_tree;
175 void BKE_mball_unlink(MetaBall *mb)
179 for (a = 0; a < mb->totcol; a++) {
180 if (mb->mat[a]) mb->mat[a]->id.us--;
186 /* do not free mball itself */
187 void BKE_mball_free(MetaBall *mb)
189 BKE_mball_unlink(mb);
192 BKE_free_animdata((ID *)mb);
195 if (mb->mat) MEM_freeN(mb->mat);
196 if (mb->bb) MEM_freeN(mb->bb);
197 BLI_freelistN(&mb->elems);
198 if (mb->disp.first) BKE_displist_free(&mb->disp);
201 MetaBall *BKE_mball_add(const char *name)
205 mb = BKE_libblock_alloc(&G.main->mball, ID_MB, name);
207 mb->size[0] = mb->size[1] = mb->size[2] = 1.0;
208 mb->texflag = MB_AUTOSPACE;
211 mb->rendersize = 0.2f;
217 MetaBall *BKE_mball_copy(MetaBall *mb)
222 mbn = BKE_libblock_copy(&mb->id);
224 BLI_duplicatelist(&mbn->elems, &mb->elems);
226 mbn->mat = MEM_dupallocN(mb->mat);
227 for (a = 0; a < mbn->totcol; a++) {
228 id_us_plus((ID *)mbn->mat[a]);
230 mbn->bb = MEM_dupallocN(mb->bb);
232 mbn->editelems = NULL;
233 mbn->lastelem = NULL;
238 static void extern_local_mball(MetaBall *mb)
241 extern_local_matarar(mb->mat, mb->totcol);
245 void BKE_mball_make_local(MetaBall *mb)
247 Main *bmain = G.main;
249 int is_local = FALSE, is_lib = FALSE;
251 /* - only lib users: do nothing
252 * - only local users: set flag
256 if (mb->id.lib == NULL) return;
257 if (mb->id.us == 1) {
258 id_clear_lib_data(bmain, &mb->id);
259 extern_local_mball(mb);
264 for (ob = G.main->object.first; ob && ELEM(0, is_lib, is_local); ob = ob->id.next) {
265 if (ob->data == mb) {
266 if (ob->id.lib) is_lib = TRUE;
267 else is_local = TRUE;
271 if (is_local && is_lib == FALSE) {
272 id_clear_lib_data(bmain, &mb->id);
273 extern_local_mball(mb);
275 else if (is_local && is_lib) {
276 MetaBall *mb_new = BKE_mball_copy(mb);
279 /* Remap paths of new ID using old library as base. */
280 BKE_id_lib_local_paths(bmain, mb->id.lib, &mb_new->id);
282 for (ob = G.main->object.first; ob; ob = ob->id.next) {
283 if (ob->data == mb) {
284 if (ob->id.lib == NULL) {
294 /* most simple meta-element adding function
295 * don't do context manipulation here (rna uses) */
296 MetaElem *BKE_mball_element_add(MetaBall *mb, const int type)
298 MetaElem *ml = MEM_callocN(sizeof(MetaElem), "metaelem");
304 ml->flag = MB_SCALE_RAD;
309 ml->expx = ml->expy = ml->expz = 1.0;
314 ml->expx = ml->expy = ml->expz = 1.0;
319 ml->expx = ml->expy = ml->expz = 1.0;
323 ml->type = MB_ELIPSOID;
331 ml->expx = ml->expy = ml->expz = 1.0;
338 BLI_addtail(&mb->elems, ml);
342 /** Compute bounding box of all MetaElems/MetaBalls.
344 * Bounding box is computed from polygonized surface. Object *ob is
345 * basic MetaBall (usually with name Meta). All other MetaBalls (with
346 * names Meta.001, Meta.002, etc) are included in this Bounding Box.
348 void BKE_mball_texspace_calc(Object *ob)
352 float *data, min[3], max[3] /*, loc[3], size[3] */;
353 int tot, do_it = FALSE;
355 if (ob->bb == NULL) ob->bb = MEM_callocN(sizeof(BoundBox), "mb boundbox");
358 /* Weird one, this. */
359 /* INIT_MINMAX(min, max); */
360 (min)[0] = (min)[1] = (min)[2] = 1.0e30f;
361 (max)[0] = (max)[1] = (max)[2] = -1.0e30f;
366 if (tot) do_it = TRUE;
369 /* Also weird... but longer. From utildefines. */
370 minmax_v3v3_v3(min, max, data);
377 min[0] = min[1] = min[2] = -1.0f;
378 max[0] = max[1] = max[2] = 1.0f;
381 loc[0] = (min[0] + max[0]) / 2.0f;
382 loc[1] = (min[1] + max[1]) / 2.0f;
383 loc[2] = (min[2] + max[2]) / 2.0f;
385 size[0] = (max[0] - min[0]) / 2.0f;
386 size[1] = (max[1] - min[1]) / 2.0f;
387 size[2] = (max[2] - min[2]) / 2.0f;
389 BKE_boundbox_init_from_minmax(bb, min, max);
392 float *BKE_mball_make_orco(Object *ob, ListBase *dispbase)
396 float *data, *orco, *orcodata;
397 float loc[3], size[3];
400 /* restore size and loc */
402 loc[0] = (bb->vec[0][0] + bb->vec[4][0]) / 2.0f;
403 size[0] = bb->vec[4][0] - loc[0];
404 loc[1] = (bb->vec[0][1] + bb->vec[2][1]) / 2.0f;
405 size[1] = bb->vec[2][1] - loc[1];
406 loc[2] = (bb->vec[0][2] + bb->vec[1][2]) / 2.0f;
407 size[2] = bb->vec[1][2] - loc[2];
409 dl = dispbase->first;
410 orcodata = MEM_mallocN(sizeof(float) * 3 * dl->nr, "MballOrco");
416 orco[0] = (data[0] - loc[0]) / size[0];
417 orco[1] = (data[1] - loc[1]) / size[1];
418 orco[2] = (data[2] - loc[2]) / size[2];
427 /* Note on mball basis stuff 2.5x (this is a can of worms)
428 * This really needs a rewrite/refactor its totally broken in anything other then basic cases
429 * Multiple Scenes + Set Scenes & mixing mball basis SHOULD work but fails to update the depsgraph on rename
430 * and linking into scenes or removal of basis mball. so take care when changing this code.
432 * Main idiot thing here is that the system returns find_basis_mball() objects which fail a is_basis_mball() test.
434 * Not only that but the depsgraph and their areas depend on this behavior!, so making small fixes here isn't worth it.
439 /** \brief Test, if Object *ob is basic MetaBall.
441 * It test last character of Object ID name. If last character
442 * is digit it return 0, else it return 1.
444 int BKE_mball_is_basis(Object *ob)
448 /* just a quick test */
449 len = strlen(ob->id.name);
450 if (isdigit(ob->id.name[len - 1]) ) return 0;
454 /* return nonzero if ob1 is a basis mball for ob */
455 int BKE_mball_is_basis_for(Object *ob1, Object *ob2)
457 int basis1nr, basis2nr;
458 char basis1name[MAX_ID_NAME], basis2name[MAX_ID_NAME];
460 BLI_split_name_num(basis1name, &basis1nr, ob1->id.name + 2, '.');
461 BLI_split_name_num(basis2name, &basis2nr, ob2->id.name + 2, '.');
463 if (!strcmp(basis1name, basis2name)) return BKE_mball_is_basis(ob1);
467 /* \brief copy some properties from object to other metaball object with same base name
469 * When some properties (wiresize, threshold, update flags) of metaball are changed, then this properties
470 * are copied to all metaballs in same "group" (metaballs with same base name: MBall,
471 * MBall.001, MBall.002, etc). The most important is to copy properties to the base metaball,
472 * because this metaball influence polygonisation of metaballs. */
473 void BKE_mball_properties_copy(Scene *scene, Object *active_object)
475 Scene *sce_iter = scene;
478 MetaBall *active_mball = (MetaBall *)active_object->data;
480 char basisname[MAX_ID_NAME], obname[MAX_ID_NAME];
482 BLI_split_name_num(basisname, &basisnr, active_object->id.name + 2, '.');
484 /* XXX recursion check, see scene.c, just too simple code this BKE_scene_base_iter_next() */
485 if (F_ERROR == BKE_scene_base_iter_next(&sce_iter, 0, NULL, NULL))
488 while (BKE_scene_base_iter_next(&sce_iter, 1, &base, &ob)) {
489 if (ob->type == OB_MBALL) {
490 if (ob != active_object) {
491 BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
493 /* Object ob has to be in same "group" ... it means, that it has to have
494 * same base of its name */
495 if (strcmp(obname, basisname) == 0) {
496 MetaBall *mb = ob->data;
498 /* Copy properties from selected/edited metaball */
499 mb->wiresize = active_mball->wiresize;
500 mb->rendersize = active_mball->rendersize;
501 mb->thresh = active_mball->thresh;
502 mb->flag = active_mball->flag;
509 /** \brief This function finds basic MetaBall.
511 * Basic MetaBall doesn't include any number at the end of
512 * its name. All MetaBalls with same base of name can be
513 * blended. MetaBalls with different basic name can't be
516 * warning!, is_basis_mball() can fail on returned object, see long note above.
518 Object *BKE_mball_basis_find(Scene *scene, Object *basis)
520 Scene *sce_iter = scene;
522 Object *ob, *bob = basis;
525 char basisname[MAX_ID_NAME], obname[MAX_ID_NAME];
527 BLI_split_name_num(basisname, &basisnr, basis->id.name + 2, '.');
530 /* XXX recursion check, see scene.c, just too simple code this BKE_scene_base_iter_next() */
531 if (F_ERROR == BKE_scene_base_iter_next(&sce_iter, 0, NULL, NULL))
534 while (BKE_scene_base_iter_next(&sce_iter, 1, &base, &ob)) {
536 if (ob->type == OB_MBALL) {
538 MetaBall *mb = ob->data;
540 /* if bob object is in edit mode, then dynamic list of all MetaElems
541 * is stored in editelems */
542 if (mb->editelems) ml = mb->editelems->first;
543 /* if bob object is in object mode */
544 else ml = mb->elems.first;
547 BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
549 /* object ob has to be in same "group" ... it means, that it has to have
550 * same base of its name */
551 if (strcmp(obname, basisname) == 0) {
552 MetaBall *mb = ob->data;
554 /* if object is in edit mode, then dynamic list of all MetaElems
555 * is stored in editelems */
556 if (mb->editelems) ml = mb->editelems->first;
557 /* if bob object is in object mode */
558 else ml = mb->elems.first;
560 if (obnr < basisnr) {
561 if (!(ob->flag & OB_FROMDUPLI)) {
569 for ( ; ml; ml = ml->next) {
570 if (!(ml->flag & MB_HIDE)) {
581 /* ******************** ARITH ************************* */
583 /* BASED AT CODE (but mostly rewritten) :
584 * C code from the article
585 * "An Implicit Surface Polygonizer"
586 * by Jules Bloomenthal, jbloom@beauty.gmu.edu
587 * in "Graphics Gems IV", Academic Press, 1994
589 * Authored by Jules Bloomenthal, Xerox PARC.
590 * Copyright (c) Xerox Corporation, 1991. All rights reserved.
591 * Permission is granted to reproduce, use and distribute this code for
592 * any and all purposes, provided that this notice appears in all copies. */
594 #define RES 12 /* # converge iterations */
596 #define L 0 /* left direction: -x, -i */
597 #define R 1 /* right direction: +x, +i */
598 #define B 2 /* bottom direction: -y, -j */
599 #define T 3 /* top direction: +y, +j */
600 #define N 4 /* near direction: -z, -k */
601 #define F 5 /* far direction: +z, +k */
602 #define LBN 0 /* left bottom near corner */
603 #define LBF 1 /* left bottom far corner */
604 #define LTN 2 /* left top near corner */
605 #define LTF 3 /* left top far corner */
606 #define RBN 4 /* right bottom near corner */
607 #define RBF 5 /* right bottom far corner */
608 #define RTN 6 /* right top near corner */
609 #define RTF 7 /* right top far corner */
611 /* the LBN corner of cube (i, j, k), corresponds with location
612 * (i-0.5)*size, (j-0.5)*size, (k-0.5)*size) */
615 #define HASHSIZE (size_t)(1 << (3 * HASHBIT)) /*! < hash table size (32768) */
617 #define HASH(i, j, k) ((((( (i) & 31) << 5) | ( (j) & 31)) << 5) | ( (k) & 31) )
619 #define MB_BIT(i, bit) (((i) >> (bit)) & 1)
620 #define FLIP(i, bit) ((i) ^ 1 << (bit)) /* flip the given bit of i */
623 /* **************** POLYGONIZATION ************************ */
625 static void calc_mballco(MetaElem *ml, float vec[3])
628 mul_m4_v3((float (*)[4])ml->mat, vec);
632 static float densfunc(MetaElem *ball, float x, float y, float z)
635 float dvec[3] = {x, y, z};
637 mul_m4_v3((float (*)[4])ball->imat, dvec);
639 switch (ball->type) {
644 if (dvec[0] > ball->len) dvec[0] -= ball->len;
645 else if (dvec[0] < -ball->len) dvec[0] += ball->len;
649 if (dvec[1] > ball->len) dvec[1] -= ball->len;
650 else if (dvec[1] < -ball->len) dvec[1] += ball->len;
654 if (dvec[2] > ball->len) dvec[2] -= ball->len;
655 else if (dvec[2] < -ball->len) dvec[2] += ball->len;
659 if (dvec[0] > ball->expx) dvec[0] -= ball->expx;
660 else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
664 if (dvec[0] > ball->expx) dvec[0] -= ball->expx;
665 else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
667 if (dvec[1] > ball->expy) dvec[1] -= ball->expy;
668 else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
672 dvec[0] /= ball->expx;
673 dvec[1] /= ball->expy;
674 dvec[2] /= ball->expz;
677 if (dvec[0] > ball->expx) dvec[0] -= ball->expx;
678 else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
681 if (dvec[1] > ball->expy) dvec[1] -= ball->expy;
682 else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
685 if (dvec[2] > ball->expz) dvec[2] -= ball->expz;
686 else if (dvec[2] < -ball->expz) dvec[2] += ball->expz;
691 dist2 = 1.0f - (len_v3(dvec) / ball->rad2);
693 if ((ball->flag & MB_NEGATIVE) == 0) {
694 return (dist2 < 0.0f) ? -0.5f : (ball->s * dist2 * dist2 * dist2) - 0.5f;
697 return (dist2 < 0.0f) ? 0.5f : 0.5f - (ball->s * dist2 * dist2 * dist2);
701 static octal_node *find_metaball_octal_node(octal_node *node, float x, float y, float z, short depth)
703 if (!depth) return node;
709 return find_metaball_octal_node(node->nodes[0], x, y, z, depth--);
715 return find_metaball_octal_node(node->nodes[1], x, y, z, depth--);
723 return find_metaball_octal_node(node->nodes[3], x, y, z, depth--);
729 return find_metaball_octal_node(node->nodes[2], x, y, z, depth--);
739 return find_metaball_octal_node(node->nodes[4], x, y, z, depth--);
745 return find_metaball_octal_node(node->nodes[5], x, y, z, depth--);
753 return find_metaball_octal_node(node->nodes[7], x, y, z, depth--);
759 return find_metaball_octal_node(node->nodes[6], x, y, z, depth--);
769 static float metaball(float x, float y, float z)
772 struct octal_node *node;
773 struct ml_pointer *ml_p;
777 if (G_mb.totelem > 1) {
778 node = find_metaball_octal_node(G_mb.metaball_tree->first, x, y, z, G_mb.metaball_tree->depth);
780 for (ml_p = node->elems.first; ml_p; ml_p = ml_p->next) {
781 dens += densfunc(ml_p->ml, x, y, z);
784 dens += -0.5f * (G_mb.metaball_tree->pos - node->pos);
785 dens += 0.5f * (G_mb.metaball_tree->neg - node->neg);
788 for (a = 0; a < G_mb.totelem; a++) {
789 dens += densfunc(G_mb.mainb[a], x, y, z);
794 dens += densfunc(G_mb.mainb[0], x, y, z);
797 return G_mb.thresh - dens;
800 /* ******************************************** */
802 static int *indices = NULL;
803 static int totindex, curindex;
806 static void accum_mballfaces(int i1, int i2, int i3, int i4)
809 /* static int i = 0; I would like to delete altogether, but I don't dare to, yet */
811 if (totindex == curindex) {
813 newi = MEM_mallocN(4 * sizeof(int) * totindex, "vertindex");
816 memcpy(newi, indices, 4 * sizeof(int) * (totindex - 256));
822 cur = indices + 4 * curindex;
824 /* displists now support array drawing, we treat tri's as fake quad */
838 /* ******************* MEMORY MANAGEMENT *********************** */
839 static void *new_pgn_element(int size)
841 /* during polygonize 1000s of elements are allocated
842 * and never freed in between. Freeing only done at the end.
844 int blocksize = 16384;
845 static int offs = 0; /* the current free address */
846 static struct pgn_elements *cur = NULL;
847 static ListBase lb = {NULL, NULL};
850 if (size > 10000 || size == 0) {
851 printf("incorrect use of new_pgn_element\n");
853 else if (size == -1) {
856 MEM_freeN(cur->data);
864 size = 4 * ( (size + 3) / 4);
867 if (size + offs < blocksize) {
868 adr = (void *) (cur->data + offs);
874 cur = MEM_callocN(sizeof(struct pgn_elements), "newpgn");
875 cur->data = MEM_callocN(blocksize, "newpgn");
876 BLI_addtail(&lb, cur);
882 static void freepolygonize(PROCESS *p)
884 MEM_freeN(p->corners);
886 MEM_freeN(p->centers);
890 if (p->vertices.ptr) MEM_freeN(p->vertices.ptr);
893 /**** Cubical Polygonization (optional) ****/
895 #define LB 0 /* left bottom edge */
896 #define LT 1 /* left top edge */
897 #define LN 2 /* left near edge */
898 #define LF 3 /* left far edge */
899 #define RB 4 /* right bottom edge */
900 #define RT 5 /* right top edge */
901 #define RN 6 /* right near edge */
902 #define RF 7 /* right far edge */
903 #define BN 8 /* bottom near edge */
904 #define BF 9 /* bottom far edge */
905 #define TN 10 /* top near edge */
906 #define TF 11 /* top far edge */
908 static INTLISTS *cubetable[256];
910 /* edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
911 static int corner1[12] = {
912 LBN, LTN, LBN, LBF, RBN, RTN, RBN, RBF, LBN, LBF, LTN, LTF
914 static int corner2[12] = {
915 LBF, LTF, LTN, LTF, RBF, RTF, RTN, RTF, RBN, RBF, RTN, RTF
917 static int leftface[12] = {
918 B, L, L, F, R, T, N, R, N, B, T, F
920 /* face on left when going corner1 to corner2 */
921 static int rightface[12] = {
922 L, T, N, L, B, R, R, F, B, F, N, T
924 /* face on right when going corner1 to corner2 */
927 /* docube: triangulate the cube directly, without decomposition */
929 static void docube(CUBE *cube, PROCESS *p, MetaBall *mb)
933 int i, index = 0, count, indexar[8];
935 for (i = 0; i < 8; i++) if (cube->corners[i]->value > 0.0f) index += (1 << i);
937 for (polys = cubetable[index]; polys; polys = polys->next) {
942 for (edges = polys->list; edges; edges = edges->next) {
943 c1 = cube->corners[corner1[edges->i]];
944 c2 = cube->corners[corner2[edges->i]];
946 indexar[count] = vertid(c1, c2, p, mb);
952 accum_mballfaces(indexar[2], indexar[1], indexar[0], 0);
955 if (indexar[0] == 0) accum_mballfaces(indexar[0], indexar[3], indexar[2], indexar[1]);
956 else accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]);
959 if (indexar[0] == 0) accum_mballfaces(indexar[0], indexar[3], indexar[2], indexar[1]);
960 else accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]);
962 accum_mballfaces(indexar[4], indexar[3], indexar[0], 0);
965 if (indexar[0] == 0) {
966 accum_mballfaces(indexar[0], indexar[3], indexar[2], indexar[1]);
967 accum_mballfaces(indexar[0], indexar[5], indexar[4], indexar[3]);
970 accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]);
971 accum_mballfaces(indexar[5], indexar[4], indexar[3], indexar[0]);
975 if (indexar[0] == 0) {
976 accum_mballfaces(indexar[0], indexar[3], indexar[2], indexar[1]);
977 accum_mballfaces(indexar[0], indexar[5], indexar[4], indexar[3]);
980 accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]);
981 accum_mballfaces(indexar[5], indexar[4], indexar[3], indexar[0]);
984 accum_mballfaces(indexar[6], indexar[5], indexar[0], 0);
993 /* testface: given cube at lattice (i, j, k), and four corners of face,
994 * if surface crosses face, compute other four corners of adjacent cube
995 * and add new cube to cube stack */
997 static void testface(int i, int j, int k, CUBE *old, int bit, int c1, int c2, int c3, int c4, PROCESS *p)
1000 CUBES *oldcubes = p->cubes;
1001 CORNER *corn1, *corn2, *corn3, *corn4;
1004 corn1 = old->corners[c1];
1005 corn2 = old->corners[c2];
1006 corn3 = old->corners[c3];
1007 corn4 = old->corners[c4];
1009 pos = corn1->value > 0.0f ? 1 : 0;
1011 /* test if no surface crossing */
1012 if ( (corn2->value > 0) == pos && (corn3->value > 0) == pos && (corn4->value > 0) == pos) return;
1013 /* test if cube out of bounds */
1014 /*if ( abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) return;*/
1015 /* test if already visited (always as last) */
1016 if (setcenter(p->centers, i, j, k)) return;
1019 /* create new cube and add cube to top of stack: */
1020 p->cubes = (CUBES *) new_pgn_element(sizeof(CUBES));
1021 p->cubes->next = oldcubes;
1026 for (n = 0; n < 8; n++) newc.corners[n] = NULL;
1028 newc.corners[FLIP(c1, bit)] = corn1;
1029 newc.corners[FLIP(c2, bit)] = corn2;
1030 newc.corners[FLIP(c3, bit)] = corn3;
1031 newc.corners[FLIP(c4, bit)] = corn4;
1033 if (newc.corners[0] == NULL) newc.corners[0] = setcorner(p, i, j, k);
1034 if (newc.corners[1] == NULL) newc.corners[1] = setcorner(p, i, j, k + 1);
1035 if (newc.corners[2] == NULL) newc.corners[2] = setcorner(p, i, j + 1, k);
1036 if (newc.corners[3] == NULL) newc.corners[3] = setcorner(p, i, j + 1, k + 1);
1037 if (newc.corners[4] == NULL) newc.corners[4] = setcorner(p, i + 1, j, k);
1038 if (newc.corners[5] == NULL) newc.corners[5] = setcorner(p, i + 1, j, k + 1);
1039 if (newc.corners[6] == NULL) newc.corners[6] = setcorner(p, i + 1, j + 1, k);
1040 if (newc.corners[7] == NULL) newc.corners[7] = setcorner(p, i + 1, j + 1, k + 1);
1042 p->cubes->cube = newc;
1045 /* setcorner: return corner with the given lattice location
1046 * set (and cache) its function value */
1048 static CORNER *setcorner(PROCESS *p, int i, int j, int k)
1050 /* for speed, do corner value caching here */
1054 /* does corner exist? */
1055 index = HASH(i, j, k);
1056 c = p->corners[index];
1058 for (; c != NULL; c = c->next) {
1059 if (c->i == i && c->j == j && c->k == k) {
1064 c = (CORNER *) new_pgn_element(sizeof(CORNER));
1067 c->co[0] = ((float)i - 0.5f) * p->size;
1069 c->co[1] = ((float)j - 0.5f) * p->size;
1071 c->co[2] = ((float)k - 0.5f) * p->size;
1072 c->value = p->function(c->co[0], c->co[1], c->co[2]);
1074 c->next = p->corners[index];
1075 p->corners[index] = c;
1081 /* nextcwedge: return next clockwise edge from given edge around given face */
1083 static int nextcwedge(int edge, int face)
1087 return (face == L) ? LF : BN;
1089 return (face == L) ? LN : TF;
1091 return (face == L) ? LB : TN;
1093 return (face == L) ? LT : BF;
1095 return (face == R) ? RN : BF;
1097 return (face == R) ? RF : TN;
1099 return (face == R) ? RT : BN;
1101 return (face == R) ? RB : TF;
1103 return (face == B) ? RB : LN;
1105 return (face == B) ? LB : RF;
1107 return (face == T) ? LT : RN;
1109 return (face == T) ? RT : LF;
1115 /* otherface: return face adjoining edge that is not the given face */
1117 static int otherface(int edge, int face)
1119 int other = leftface[edge];
1120 return face == other ? rightface[edge] : other;
1124 /* makecubetable: create the 256 entry table for cubical polygonization */
1126 static void makecubetable(void)
1128 static int is_done = FALSE;
1129 int i, e, c, done[12], pos[8];
1131 if (is_done) return;
1134 for (i = 0; i < 256; i++) {
1135 for (e = 0; e < 12; e++) done[e] = 0;
1136 for (c = 0; c < 8; c++) pos[c] = MB_BIT(i, c);
1137 for (e = 0; e < 12; e++)
1138 if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
1139 INTLIST *ints = NULL;
1140 INTLISTS *lists = (INTLISTS *) MEM_callocN(sizeof(INTLISTS), "mball_intlist");
1141 int start = e, edge = e;
1143 /* get face that is to right of edge from pos to neg corner: */
1144 int face = pos[corner1[e]] ? rightface[e] : leftface[e];
1147 edge = nextcwedge(edge, face);
1149 if (pos[corner1[edge]] != pos[corner2[edge]]) {
1150 INTLIST *tmp = ints;
1152 ints = (INTLIST *) MEM_callocN(sizeof(INTLIST), "mball_intlist");
1154 ints->next = tmp; /* add edge to head of list */
1156 if (edge == start) break;
1157 face = otherface(edge, face);
1160 lists->list = ints; /* add ints to head of table entry */
1161 lists->next = cubetable[i];
1162 cubetable[i] = lists;
1167 void BKE_mball_cubeTable_free(void)
1170 INTLISTS *lists, *nlists;
1171 INTLIST *ints, *nints;
1173 for (i = 0; i < 256; i++) {
1174 lists = cubetable[i];
1176 nlists = lists->next;
1188 cubetable[i] = NULL;
1194 /* setcenter: set (i, j, k) entry of table[]
1195 * return 1 if already set; otherwise, set and return 0 */
1197 static int setcenter(CENTERLIST *table[], const int i, const int j, const int k)
1200 CENTERLIST *newc, *l, *q;
1202 index = HASH(i, j, k);
1205 for (l = q; l != NULL; l = l->next) {
1206 if (l->i == i && l->j == j && l->k == k) return 1;
1209 newc = (CENTERLIST *) new_pgn_element(sizeof(CENTERLIST));
1214 table[index] = newc;
1220 /* setedge: set vertex id for edge */
1222 static void setedge(EDGELIST *table[],
1231 if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
1242 index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
1243 newe = (EDGELIST *) new_pgn_element(sizeof(EDGELIST));
1251 newe->next = table[index];
1252 table[index] = newe;
1256 /* getedge: return vertex id for edge; return -1 if not set */
1258 static int getedge(EDGELIST *table[],
1259 int i1, int j1, int k1,
1260 int i2, int j2, int k2)
1264 if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
1275 q = table[HASH(i1, j1, k1) + HASH(i2, j2, k2)];
1276 for (; q != NULL; q = q->next) {
1277 if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 &&
1278 q->i2 == i2 && q->j2 == j2 && q->k2 == k2)
1287 /**** Vertices ****/
1293 /* vertid: return index for vertex on edge:
1294 * c1->value and c2->value are presumed of different sign
1295 * return saved index if any; else compute vertex and save */
1297 /* addtovertices: add v to sequence of vertices */
1299 static void addtovertices(VERTICES *vertices, VERTEX v)
1301 if (vertices->count == vertices->max) {
1304 vertices->max = vertices->count == 0 ? 10 : 2 * vertices->count;
1305 newv = (VERTEX *) MEM_callocN(vertices->max * sizeof(VERTEX), "addtovertices");
1307 for (i = 0; i < vertices->count; i++) newv[i] = vertices->ptr[i];
1309 if (vertices->ptr != NULL) MEM_freeN(vertices->ptr);
1310 vertices->ptr = newv;
1312 vertices->ptr[vertices->count++] = v;
1315 /* vnormal: compute unit length surface normal at point */
1317 static void vnormal(const float point[3], PROCESS *p, float r_no[3])
1319 float delta = 0.2f * p->delta;
1320 float f = p->function(point[0], point[1], point[2]);
1322 r_no[0] = p->function(point[0] + delta, point[1], point[2]) - f;
1323 r_no[1] = p->function(point[0], point[1] + delta, point[2]) - f;
1324 r_no[2] = p->function(point[0], point[1], point[2] + delta) - f;
1325 f = normalize_v3(r_no);
1332 f = p->function(point[0], point[1], point[2]);
1334 tvec[0] = p->function(point[0] + delta, point[1], point[2]) - f;
1335 tvec[1] = p->function(point[0], point[1] + delta, point[2]) - f;
1336 tvec[2] = p->function(point[0], point[1], point[2] + delta) - f;
1338 if (normalize_v3(tvec) != 0.0f) {
1339 add_v3_v3(r_no, tvec);
1346 static int vertid(const CORNER *c1, const CORNER *c2, PROCESS *p, MetaBall *mb)
1349 int vid = getedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
1352 return vid; /* previously computed */
1355 converge(c1->co, c2->co, c1->value, c2->value, p->function, v.co, mb, 1); /* position */
1356 vnormal(v.co, p, v.no);
1358 addtovertices(&p->vertices, v); /* save vertex */
1359 vid = p->vertices.count - 1;
1360 setedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
1368 /* converge: from two points of differing sign, converge to zero crossing */
1369 /* watch it: p1 and p2 are used to calculate */
1370 static void converge(const float p1[3], const float p2[3], float v1, float v2,
1371 float (*function)(float, float, float), float p[3], MetaBall *mb, int f)
1374 float pos[3], neg[3];
1375 float positive = 0.0f, negative = 0.0f;
1379 copy_v3_v3(pos, p2);
1380 copy_v3_v3(neg, p1);
1385 copy_v3_v3(pos, p1);
1386 copy_v3_v3(neg, p2);
1391 sub_v3_v3v3(dvec, pos, neg);
1393 /* Approximation by linear interpolation is faster then binary subdivision,
1394 * but it results sometimes (mb->thresh < 0.2) into the strange results */
1395 if ((mb->thresh > 0.2f) && (f == 1)) {
1396 if ((dvec[1] == 0.0f) && (dvec[2] == 0.0f)) {
1397 p[0] = neg[0] - negative * dvec[0] / (positive - negative);
1402 if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
1404 p[1] = neg[1] - negative * dvec[1] / (positive - negative);
1408 if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
1411 p[2] = neg[2] - negative * dvec[2] / (positive - negative);
1416 if ((dvec[1] == 0.0f) && (dvec[2] == 0.0f)) {
1420 if (i++ == RES) return;
1421 p[0] = 0.5f * (pos[0] + neg[0]);
1422 if ((function(p[0], p[1], p[2])) > 0.0f) pos[0] = p[0]; else neg[0] = p[0];
1426 if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
1430 if (i++ == RES) return;
1431 p[1] = 0.5f * (pos[1] + neg[1]);
1432 if ((function(p[0], p[1], p[2])) > 0.0f) pos[1] = p[1]; else neg[1] = p[1];
1436 if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
1440 if (i++ == RES) return;
1441 p[2] = 0.5f * (pos[2] + neg[2]);
1442 if ((function(p[0], p[1], p[2])) > 0.0f) pos[2] = p[2]; else neg[2] = p[2];
1446 /* This is necessary to find start point */
1448 mid_v3_v3v3(&p[0], pos, neg);
1454 if ((function(p[0], p[1], p[2])) > 0.0f) {
1455 copy_v3_v3(pos, &p[0]);
1458 copy_v3_v3(neg, &p[0]);
1463 /* ************************************** */
1464 static void add_cube(PROCESS *mbproc, int i, int j, int k, int count)
1470 /* hmmm, not only one, but eight cube will be added on the stack
1472 for (a = i - 1; a < i + count; a++)
1473 for (b = j - 1; b < j + count; b++)
1474 for (c = k - 1; c < k + count; c++) {
1475 /* test if cube has been found before */
1476 if (setcenter(mbproc->centers, a, b, c) == 0) {
1477 /* push cube on stack: */
1478 ncube = (CUBES *) new_pgn_element(sizeof(CUBES));
1479 ncube->next = mbproc->cubes;
1480 mbproc->cubes = ncube;
1486 /* set corners of initial cube: */
1487 for (n = 0; n < 8; n++)
1488 ncube->cube.corners[n] = setcorner(mbproc, a + MB_BIT(n, 2), b + MB_BIT(n, 1), c + MB_BIT(n, 0));
1494 static void find_first_points(PROCESS *mbproc, MetaBall *mb, int a)
1500 f = 1.0 - (mb->thresh / ml->s);
1502 /* Skip, when Stiffness of MetaElement is too small ... MetaElement can't be
1503 * visible alone ... but still can influence others MetaElements :-) */
1505 float IN[3] = {0.0f}, OUT[3] = {0.0f}, in[3] = {0.0f}, out[3];
1506 int i, j, k, c_i, c_j, c_k;
1507 int index[3] = {1, 0, -1};
1508 float in_v /*, out_v*/;
1511 float tmp_v, workp_v, max_len, len, nx, ny, nz, MAXN;
1513 calc_mballco(ml, in);
1514 in_v = mbproc->function(in[0], in[1], in[2]);
1516 for (i = 0; i < 3; i++) {
1519 OUT[0] = out[0] = IN[0] + index[i] * ml->rad;
1525 OUT[0] = out[0] = IN[0] + index[i] * (ml->expx + ml->rad);
1529 for (j = 0; j < 3; j++) {
1532 OUT[1] = out[1] = IN[1] + index[j] * ml->rad;
1538 OUT[1] = out[1] = IN[1] + index[j] * (ml->expy + ml->rad);
1542 for (k = 0; k < 3; k++) {
1549 out[2] = IN[2] + index[k] * ml->rad;
1553 out[2] = IN[2] + index[k] * (ml->expz + ml->rad);
1557 calc_mballco(ml, out);
1559 /*out_v = mbproc->function(out[0], out[1], out[2]);*/ /*UNUSED*/
1561 /* find "first points" on Implicit Surface of MetaElemnt ml */
1562 copy_v3_v3(workp, in);
1564 max_len = len_v3v3(out, in);
1566 nx = abs((out[0] - in[0]) / mbproc->size);
1567 ny = abs((out[1] - in[1]) / mbproc->size);
1568 nz = abs((out[2] - in[2]) / mbproc->size);
1570 MAXN = MAX3(nx, ny, nz);
1572 dvec[0] = (out[0] - in[0]) / MAXN;
1573 dvec[1] = (out[1] - in[1]) / MAXN;
1574 dvec[2] = (out[2] - in[2]) / MAXN;
1577 while (len <= max_len) {
1578 workp[0] += dvec[0];
1579 workp[1] += dvec[1];
1580 workp[2] += dvec[2];
1581 /* compute value of implicite function */
1582 tmp_v = mbproc->function(workp[0], workp[1], workp[2]);
1583 /* add cube to the stack, when value of implicite function crosses zero value */
1584 if ((tmp_v < 0.0f && workp_v >= 0.0f) || (tmp_v > 0.0f && workp_v <= 0.0f)) {
1586 /* indexes of CUBE, which includes "first point" */
1587 c_i = (int)floor(workp[0] / mbproc->size);
1588 c_j = (int)floor(workp[1] / mbproc->size);
1589 c_k = (int)floor(workp[2] / mbproc->size);
1591 /* add CUBE (with indexes c_i, c_j, c_k) to the stack,
1592 * this cube includes found point of Implicit Surface */
1593 if ((ml->flag & MB_NEGATIVE) == 0) {
1594 add_cube(mbproc, c_i, c_j, c_k, 1);
1597 add_cube(mbproc, c_i, c_j, c_k, 2);
1600 len = len_v3v3(workp, in);
1611 static void polygonize(PROCESS *mbproc, MetaBall *mb)
1616 mbproc->vertices.count = mbproc->vertices.max = 0;
1617 mbproc->vertices.ptr = NULL;
1619 /* allocate hash tables and build cube polygon table: */
1620 mbproc->centers = MEM_callocN(HASHSIZE * sizeof(CENTERLIST *), "mbproc->centers");
1621 mbproc->corners = MEM_callocN(HASHSIZE * sizeof(CORNER *), "mbproc->corners");
1622 mbproc->edges = MEM_callocN(2 * HASHSIZE * sizeof(EDGELIST *), "mbproc->edges");
1625 for (a = 0; a < G_mb.totelem; a++) {
1627 /* try to find 8 points on the surface for each MetaElem */
1628 find_first_points(mbproc, mb, a);
1631 /* polygonize all MetaElems of current MetaBall */
1632 while (mbproc->cubes != NULL) { /* process active cubes till none left */
1633 c = mbproc->cubes->cube;
1635 /* polygonize the cube directly: */
1636 docube(&c, mbproc, mb);
1638 /* pop current cube from stack */
1639 mbproc->cubes = mbproc->cubes->next;
1641 /* test six face directions, maybe add to stack: */
1642 testface(c.i - 1, c.j, c.k, &c, 2, LBN, LBF, LTN, LTF, mbproc);
1643 testface(c.i + 1, c.j, c.k, &c, 2, RBN, RBF, RTN, RTF, mbproc);
1644 testface(c.i, c.j - 1, c.k, &c, 1, LBN, LBF, RBN, RBF, mbproc);
1645 testface(c.i, c.j + 1, c.k, &c, 1, LTN, LTF, RTN, RTF, mbproc);
1646 testface(c.i, c.j, c.k - 1, &c, 0, LBN, LTN, RBN, RTN, mbproc);
1647 testface(c.i, c.j, c.k + 1, &c, 0, LBF, LTF, RBF, RTF, mbproc);
1651 static float init_meta(Scene *scene, Object *ob) /* return totsize */
1653 Scene *sce_iter = scene;
1658 float size, totsize, obinv[4][4], obmat[4][4], vec[3];
1660 int a, obnr, zero_size = 0;
1661 char obname[MAX_ID_NAME];
1663 copy_m4_m4(obmat, ob->obmat); /* to cope with duplicators from BKE_scene_base_iter_next */
1664 invert_m4_m4(obinv, ob->obmat);
1667 BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
1669 /* make main array */
1670 BKE_scene_base_iter_next(&sce_iter, 0, NULL, NULL);
1671 while (BKE_scene_base_iter_next(&sce_iter, 1, &base, &bob)) {
1673 if (bob->type == OB_MBALL) {
1677 if (bob == ob && (base->flag & OB_FROMDUPLI) == 0) {
1680 if (mb->editelems) ml = mb->editelems->first;
1681 else ml = mb->elems.first;
1684 char name[MAX_ID_NAME];
1687 BLI_split_name_num(name, &nr, bob->id.name + 2, '.');
1688 if (strcmp(obname, name) == 0) {
1691 if (mb->editelems) ml = mb->editelems->first;
1692 else ml = mb->elems.first;
1696 /* when metaball object has zero scale, then MetaElem to this MetaBall
1697 * will not be put to mainb array */
1698 if (bob->size[0] == 0.0f || bob->size[1] == 0.0f || bob->size[2] == 0.0f) {
1701 else if (bob->parent) {
1702 struct Object *pob = bob->parent;
1704 if (pob->size[0] == 0.0f || pob->size[1] == 0.0f || pob->size[2] == 0.0f) {
1713 unsigned int ml_count = 0;
1718 G_mb.totelem -= ml_count;
1722 if (!(ml->flag & MB_HIDE)) {
1724 float temp1[4][4], temp2[4][4], temp3[4][4];
1725 float (*mat)[4] = NULL, (*imat)[4] = NULL;
1726 float max_x, max_y, max_z, min_x, min_y, min_z;
1728 max_x = max_y = max_z = -3.4e38;
1729 min_x = min_y = min_z = 3.4e38;
1731 /* too big stiffness seems only ugly due to linear interpolation
1732 * no need to have possibility for too big stiffness */
1733 if (ml->s > 10.0f) ml->s = 10.0f;
1735 /* Rotation of MetaElem is stored in quat */
1736 quat_to_mat4(temp3, ml->quat);
1738 /* Translation of MetaElem */
1740 temp2[3][0] = ml->x;
1741 temp2[3][1] = ml->y;
1742 temp2[3][2] = ml->z;
1744 mult_m4_m4m4(temp1, temp2, temp3);
1746 /* make a copy because of duplicates */
1747 G_mb.mainb[a] = new_pgn_element(sizeof(MetaElem));
1748 *(G_mb.mainb[a]) = *ml;
1749 G_mb.mainb[a]->bb = new_pgn_element(sizeof(BoundBox));
1751 mat = new_pgn_element(4 * 4 * sizeof(float));
1752 imat = new_pgn_element(4 * 4 * sizeof(float));
1754 /* mat is the matrix to transform from mball into the basis-mball */
1755 invert_m4_m4(obinv, obmat);
1756 mult_m4_m4m4(temp2, obinv, bob->obmat);
1757 /* MetaBall transformation */
1758 mult_m4_m4m4(mat, temp2, temp1);
1760 invert_m4_m4(imat, mat);
1762 G_mb.mainb[a]->rad2 = ml->rad * ml->rad;
1764 G_mb.mainb[a]->mat = (float *) mat;
1765 G_mb.mainb[a]->imat = (float *) imat;
1767 /* untransformed Bounding Box of MetaElem */
1769 G_mb.mainb[a]->bb->vec[0][0] = -ml->expx;
1770 G_mb.mainb[a]->bb->vec[0][1] = -ml->expy;
1771 G_mb.mainb[a]->bb->vec[0][2] = -ml->expz;
1773 G_mb.mainb[a]->bb->vec[1][0] = ml->expx;
1774 G_mb.mainb[a]->bb->vec[1][1] = -ml->expy;
1775 G_mb.mainb[a]->bb->vec[1][2] = -ml->expz;
1777 G_mb.mainb[a]->bb->vec[2][0] = ml->expx;
1778 G_mb.mainb[a]->bb->vec[2][1] = ml->expy;
1779 G_mb.mainb[a]->bb->vec[2][2] = -ml->expz;
1781 G_mb.mainb[a]->bb->vec[3][0] = -ml->expx;
1782 G_mb.mainb[a]->bb->vec[3][1] = ml->expy;
1783 G_mb.mainb[a]->bb->vec[3][2] = -ml->expz;
1785 G_mb.mainb[a]->bb->vec[4][0] = -ml->expx;
1786 G_mb.mainb[a]->bb->vec[4][1] = -ml->expy;
1787 G_mb.mainb[a]->bb->vec[4][2] = ml->expz;
1789 G_mb.mainb[a]->bb->vec[5][0] = ml->expx;
1790 G_mb.mainb[a]->bb->vec[5][1] = -ml->expy;
1791 G_mb.mainb[a]->bb->vec[5][2] = ml->expz;
1793 G_mb.mainb[a]->bb->vec[6][0] = ml->expx;
1794 G_mb.mainb[a]->bb->vec[6][1] = ml->expy;
1795 G_mb.mainb[a]->bb->vec[6][2] = ml->expz;
1797 G_mb.mainb[a]->bb->vec[7][0] = -ml->expx;
1798 G_mb.mainb[a]->bb->vec[7][1] = ml->expy;
1799 G_mb.mainb[a]->bb->vec[7][2] = ml->expz;
1801 /* transformation of Metalem bb */
1802 for (i = 0; i < 8; i++)
1803 mul_m4_v3((float (*)[4])mat, G_mb.mainb[a]->bb->vec[i]);
1805 /* find max and min of transformed bb */
1806 for (i = 0; i < 8; i++) {
1808 if (G_mb.mainb[a]->bb->vec[i][0] > max_x) max_x = G_mb.mainb[a]->bb->vec[i][0];
1809 if (G_mb.mainb[a]->bb->vec[i][1] > max_y) max_y = G_mb.mainb[a]->bb->vec[i][1];
1810 if (G_mb.mainb[a]->bb->vec[i][2] > max_z) max_z = G_mb.mainb[a]->bb->vec[i][2];
1812 if (G_mb.mainb[a]->bb->vec[i][0] < min_x) min_x = G_mb.mainb[a]->bb->vec[i][0];
1813 if (G_mb.mainb[a]->bb->vec[i][1] < min_y) min_y = G_mb.mainb[a]->bb->vec[i][1];
1814 if (G_mb.mainb[a]->bb->vec[i][2] < min_z) min_z = G_mb.mainb[a]->bb->vec[i][2];
1817 /* create "new" bb, only point 0 and 6, which are
1818 * necessary for octal tree filling */
1819 G_mb.mainb[a]->bb->vec[0][0] = min_x - ml->rad;
1820 G_mb.mainb[a]->bb->vec[0][1] = min_y - ml->rad;
1821 G_mb.mainb[a]->bb->vec[0][2] = min_z - ml->rad;
1823 G_mb.mainb[a]->bb->vec[6][0] = max_x + ml->rad;
1824 G_mb.mainb[a]->bb->vec[6][1] = max_y + ml->rad;
1825 G_mb.mainb[a]->bb->vec[6][2] = max_z + ml->rad;
1836 /* totsize (= 'manhattan' radius) */
1838 for (a = 0; a < G_mb.totelem; a++) {
1840 vec[0] = G_mb.mainb[a]->x + G_mb.mainb[a]->rad + G_mb.mainb[a]->expx;
1841 vec[1] = G_mb.mainb[a]->y + G_mb.mainb[a]->rad + G_mb.mainb[a]->expy;
1842 vec[2] = G_mb.mainb[a]->z + G_mb.mainb[a]->rad + G_mb.mainb[a]->expz;
1844 calc_mballco(G_mb.mainb[a], vec);
1846 size = fabsf(vec[0]);
1847 if (size > totsize) totsize = size;
1848 size = fabsf(vec[1]);
1849 if (size > totsize) totsize = size;
1850 size = fabsf(vec[2]);
1851 if (size > totsize) totsize = size;
1853 vec[0] = G_mb.mainb[a]->x - G_mb.mainb[a]->rad;
1854 vec[1] = G_mb.mainb[a]->y - G_mb.mainb[a]->rad;
1855 vec[2] = G_mb.mainb[a]->z - G_mb.mainb[a]->rad;
1857 calc_mballco(G_mb.mainb[a], vec);
1859 size = fabsf(vec[0]);
1860 if (size > totsize) totsize = size;
1861 size = fabsf(vec[1]);
1862 if (size > totsize) totsize = size;
1863 size = fabsf(vec[2]);
1864 if (size > totsize) totsize = size;
1867 for (a = 0; a < G_mb.totelem; a++) {
1868 G_mb.thresh += densfunc(G_mb.mainb[a], 2.0f * totsize, 2.0f * totsize, 2.0f * totsize);
1874 /* if MetaElem lies in node, then node includes MetaElem pointer (ml_p)
1875 * pointing at MetaElem (ml)
1877 static void fill_metaball_octal_node(octal_node *node, MetaElem *ml, short i)
1881 ml_p = MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
1883 BLI_addtail(&(node->nodes[i]->elems), ml_p);
1886 if ((ml->flag & MB_NEGATIVE) == 0) {
1887 node->nodes[i]->pos++;
1890 node->nodes[i]->neg++;
1894 /* Node is subdivided as is illustrated on the following figure:
1900 * +------+------+ |/|
1903 * +------+------+ |/
1909 static void subdivide_metaball_octal_node(octal_node *node, float size_x, float size_y, float size_z, short depth)
1916 /* create new nodes */
1917 for (a = 0; a < 8; a++) {
1918 node->nodes[a] = MEM_mallocN(sizeof(octal_node), "octal_node");
1919 for (i = 0; i < 8; i++)
1920 node->nodes[a]->nodes[i] = NULL;
1921 node->nodes[a]->parent = node;
1922 node->nodes[a]->elems.first = NULL;
1923 node->nodes[a]->elems.last = NULL;
1924 node->nodes[a]->count = 0;
1925 node->nodes[a]->neg = 0;
1926 node->nodes[a]->pos = 0;
1933 /* center of node */
1934 node->x = x = node->x_min + size_x;
1935 node->y = y = node->y_min + size_y;
1936 node->z = z = node->z_min + size_z;
1938 /* setting up of border points of new nodes */
1939 node->nodes[0]->x_min = node->x_min;
1940 node->nodes[0]->y_min = node->y_min;
1941 node->nodes[0]->z_min = node->z_min;
1942 node->nodes[0]->x = node->nodes[0]->x_min + size_x / 2;
1943 node->nodes[0]->y = node->nodes[0]->y_min + size_y / 2;
1944 node->nodes[0]->z = node->nodes[0]->z_min + size_z / 2;
1946 node->nodes[1]->x_min = x;
1947 node->nodes[1]->y_min = node->y_min;
1948 node->nodes[1]->z_min = node->z_min;
1949 node->nodes[1]->x = node->nodes[1]->x_min + size_x / 2;
1950 node->nodes[1]->y = node->nodes[1]->y_min + size_y / 2;
1951 node->nodes[1]->z = node->nodes[1]->z_min + size_z / 2;
1953 node->nodes[2]->x_min = x;
1954 node->nodes[2]->y_min = y;
1955 node->nodes[2]->z_min = node->z_min;
1956 node->nodes[2]->x = node->nodes[2]->x_min + size_x / 2;
1957 node->nodes[2]->y = node->nodes[2]->y_min + size_y / 2;
1958 node->nodes[2]->z = node->nodes[2]->z_min + size_z / 2;
1960 node->nodes[3]->x_min = node->x_min;
1961 node->nodes[3]->y_min = y;
1962 node->nodes[3]->z_min = node->z_min;
1963 node->nodes[3]->x = node->nodes[3]->x_min + size_x / 2;
1964 node->nodes[3]->y = node->nodes[3]->y_min + size_y / 2;
1965 node->nodes[3]->z = node->nodes[3]->z_min + size_z / 2;
1967 node->nodes[4]->x_min = node->x_min;
1968 node->nodes[4]->y_min = node->y_min;
1969 node->nodes[4]->z_min = z;
1970 node->nodes[4]->x = node->nodes[4]->x_min + size_x / 2;
1971 node->nodes[4]->y = node->nodes[4]->y_min + size_y / 2;
1972 node->nodes[4]->z = node->nodes[4]->z_min + size_z / 2;
1974 node->nodes[5]->x_min = x;
1975 node->nodes[5]->y_min = node->y_min;
1976 node->nodes[5]->z_min = z;
1977 node->nodes[5]->x = node->nodes[5]->x_min + size_x / 2;
1978 node->nodes[5]->y = node->nodes[5]->y_min + size_y / 2;
1979 node->nodes[5]->z = node->nodes[5]->z_min + size_z / 2;
1981 node->nodes[6]->x_min = x;
1982 node->nodes[6]->y_min = y;
1983 node->nodes[6]->z_min = z;
1984 node->nodes[6]->x = node->nodes[6]->x_min + size_x / 2;
1985 node->nodes[6]->y = node->nodes[6]->y_min + size_y / 2;
1986 node->nodes[6]->z = node->nodes[6]->z_min + size_z / 2;
1988 node->nodes[7]->x_min = node->x_min;
1989 node->nodes[7]->y_min = y;
1990 node->nodes[7]->z_min = z;
1991 node->nodes[7]->x = node->nodes[7]->x_min + size_x / 2;
1992 node->nodes[7]->y = node->nodes[7]->y_min + size_y / 2;
1993 node->nodes[7]->z = node->nodes[7]->z_min + size_z / 2;
1995 ml_p = node->elems.first;
1997 /* setting up references of MetaElems for new nodes */
2000 if (ml->bb->vec[0][2] < z) {
2001 if (ml->bb->vec[0][1] < y) {
2002 /* vec[0][0] lies in first octant */
2003 if (ml->bb->vec[0][0] < x) {
2004 /* ml belongs to the (0)1st node */
2005 fill_metaball_octal_node(node, ml, 0);
2007 /* ml belongs to the (3)4th node */
2008 if (ml->bb->vec[6][1] >= y) {
2009 fill_metaball_octal_node(node, ml, 3);
2011 /* ml belongs to the (7)8th node */
2012 if (ml->bb->vec[6][2] >= z) {
2013 fill_metaball_octal_node(node, ml, 7);
2017 /* ml belongs to the (1)2nd node */
2018 if (ml->bb->vec[6][0] >= x) {
2019 fill_metaball_octal_node(node, ml, 1);
2021 /* ml belongs to the (5)6th node */
2022 if (ml->bb->vec[6][2] >= z) {
2023 fill_metaball_octal_node(node, ml, 5);
2027 /* ml belongs to the (2)3th node */
2028 if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
2029 fill_metaball_octal_node(node, ml, 2);
2031 /* ml belong to the (6)7th node */
2032 if (ml->bb->vec[6][2] >= z) {
2033 fill_metaball_octal_node(node, ml, 6);
2038 /* ml belongs to the (4)5th node too */
2039 if (ml->bb->vec[6][2] >= z) {
2040 fill_metaball_octal_node(node, ml, 4);
2046 /* vec[0][0] is in the (1)second octant */
2048 /* ml belong to the (1)2nd node */
2049 fill_metaball_octal_node(node, ml, 1);
2051 /* ml belongs to the (2)3th node */
2052 if (ml->bb->vec[6][1] >= y) {
2053 fill_metaball_octal_node(node, ml, 2);
2055 /* ml belongs to the (6)7th node */
2056 if (ml->bb->vec[6][2] >= z) {
2057 fill_metaball_octal_node(node, ml, 6);
2062 /* ml belongs to the (5)6th node */
2063 if (ml->bb->vec[6][2] >= z) {
2064 fill_metaball_octal_node(node, ml, 5);
2069 /* vec[0][0] is in the (3)4th octant */
2070 if (ml->bb->vec[0][0] < x) {
2071 /* ml belongs to the (3)4nd node */
2072 fill_metaball_octal_node(node, ml, 3);
2074 /* ml belongs to the (7)8th node */
2075 if (ml->bb->vec[6][2] >= z) {
2076 fill_metaball_octal_node(node, ml, 7);
2080 /* ml belongs to the (2)3th node */
2081 if (ml->bb->vec[6][0] >= x) {
2082 fill_metaball_octal_node(node, ml, 2);
2084 /* ml belongs to the (6)7th node */
2085 if (ml->bb->vec[6][2] >= z) {
2086 fill_metaball_octal_node(node, ml, 6);
2093 /* vec[0][0] is in the (2)3th octant */
2094 if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
2095 /* ml belongs to the (2)3th node */
2096 fill_metaball_octal_node(node, ml, 2);
2098 /* ml belongs to the (6)7th node */
2099 if (ml->bb->vec[6][2] >= z) {
2100 fill_metaball_octal_node(node, ml, 6);
2105 if (ml->bb->vec[0][1] < y) {
2106 /* vec[0][0] lies in (4)5th octant */
2107 if (ml->bb->vec[0][0] < x) {
2108 /* ml belongs to the (4)5th node */
2109 fill_metaball_octal_node(node, ml, 4);
2111 if (ml->bb->vec[6][0] >= x) {
2112 fill_metaball_octal_node(node, ml, 5);
2115 if (ml->bb->vec[6][1] >= y) {
2116 fill_metaball_octal_node(node, ml, 7);
2119 if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
2120 fill_metaball_octal_node(node, ml, 6);
2123 /* vec[0][0] lies in (5)6th octant */
2125 fill_metaball_octal_node(node, ml, 5);
2127 if (ml->bb->vec[6][1] >= y) {
2128 fill_metaball_octal_node(node, ml, 6);
2133 /* vec[0][0] lies in (7)8th octant */
2134 if (ml->bb->vec[0][0] < x) {
2135 fill_metaball_octal_node(node, ml, 7);
2137 if (ml->bb->vec[6][0] >= x) {
2138 fill_metaball_octal_node(node, ml, 6);
2144 /* vec[0][0] lies in (6)7th octant */
2145 if ((ml->bb->vec[0][0] >= x) && (ml->bb->vec[0][1] >= y)) {
2146 fill_metaball_octal_node(node, ml, 6);
2152 /* free references of MetaElems for curent node (it is not needed anymore) */
2153 BLI_freelistN(&node->elems);
2158 for (a = 0; a < 8; a++) {
2159 if (node->nodes[a]->count > 0) /* if node is not empty, then it is subdivided */
2160 subdivide_metaball_octal_node(node->nodes[a], size_x, size_y, size_z, depth);
2165 /* free all octal nodes recursively */
2166 static void free_metaball_octal_node(octal_node *node)
2169 for (a = 0; a < 8; a++) {
2170 if (node->nodes[a] != NULL) free_metaball_octal_node(node->nodes[a]);
2172 BLI_freelistN(&node->elems);
2176 /* If scene include more then one MetaElem, then octree is used */
2177 static void init_metaball_octal_tree(int depth)
2179 struct octal_node *node;
2184 G_mb.metaball_tree = MEM_mallocN(sizeof(octal_tree), "metaball_octal_tree");
2185 G_mb.metaball_tree->first = node = MEM_mallocN(sizeof(octal_node), "metaball_octal_node");
2186 /* maximal depth of octree */
2187 G_mb.metaball_tree->depth = depth;
2189 G_mb.metaball_tree->neg = node->neg = 0;
2190 G_mb.metaball_tree->pos = node->pos = 0;
2192 node->elems.first = NULL;
2193 node->elems.last = NULL;
2196 for (a = 0; a < 8; a++)
2197 node->nodes[a] = NULL;
2199 node->x_min = node->y_min = node->z_min = FLT_MAX;
2200 node->x_max = node->y_max = node->z_max = -FLT_MAX;
2202 /* size of octal tree scene */
2203 for (a = 0; a < G_mb.totelem; a++) {
2204 if (G_mb.mainb[a]->bb->vec[0][0] < node->x_min) node->x_min = G_mb.mainb[a]->bb->vec[0][0];
2205 if (G_mb.mainb[a]->bb->vec[0][1] < node->y_min) node->y_min = G_mb.mainb[a]->bb->vec[0][1];
2206 if (G_mb.mainb[a]->bb->vec[0][2] < node->z_min) node->z_min = G_mb.mainb[a]->bb->vec[0][2];
2208 if (G_mb.mainb[a]->bb->vec[6][0] > node->x_max) node->x_max = G_mb.mainb[a]->bb->vec[6][0];
2209 if (G_mb.mainb[a]->bb->vec[6][1] > node->y_max) node->y_max = G_mb.mainb[a]->bb->vec[6][1];
2210 if (G_mb.mainb[a]->bb->vec[6][2] > node->z_max) node->z_max = G_mb.mainb[a]->bb->vec[6][2];
2212 ml_p = MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
2213 ml_p->ml = G_mb.mainb[a];
2214 BLI_addtail(&node->elems, ml_p);
2216 if ((G_mb.mainb[a]->flag & MB_NEGATIVE) == 0) {
2217 /* number of positive MetaElem in scene */
2218 G_mb.metaball_tree->pos++;
2221 /* number of negative MetaElem in scene */
2222 G_mb.metaball_tree->neg++;
2226 /* size of first node */
2227 size[0] = node->x_max - node->x_min;
2228 size[1] = node->y_max - node->y_min;
2229 size[2] = node->z_max - node->z_min;
2231 /* first node is subdivided recursively */
2232 subdivide_metaball_octal_node(node, size[0], size[1], size[2], G_mb.metaball_tree->depth);
2235 void BKE_mball_polygonize(Scene *scene, Object *ob, ListBase *dispbase)
2241 float *co, *no, totsize, width;
2245 if (G_mb.totelem == 0) return;
2246 if ((G.is_rendering == FALSE) && (mb->flag == MB_UPDATE_NEVER)) return;
2247 if (G.moving && mb->flag == MB_UPDATE_FAST) return;
2249 curindex = totindex = 0;
2251 G_mb.thresh = mb->thresh;
2253 /* total number of MetaElems (totelem) is precomputed in find_basis_mball() function */
2254 G_mb.mainb = MEM_mallocN(sizeof(void *) * G_mb.totelem, "mainb");
2256 /* initialize all mainb (MetaElems) */
2257 totsize = init_meta(scene, ob);
2259 if (G_mb.metaball_tree) {
2260 free_metaball_octal_node(G_mb.metaball_tree->first);
2261 MEM_freeN(G_mb.metaball_tree);
2262 G_mb.metaball_tree = NULL;
2265 /* if scene includes more then one MetaElem, then octal tree optimization is used */
2266 if ((G_mb.totelem > 1) && (G_mb.totelem <= 64)) init_metaball_octal_tree(1);
2267 if ((G_mb.totelem > 64) && (G_mb.totelem <= 128)) init_metaball_octal_tree(2);
2268 if ((G_mb.totelem > 128) && (G_mb.totelem <= 512)) init_metaball_octal_tree(3);
2269 if ((G_mb.totelem > 512) && (G_mb.totelem <= 1024)) init_metaball_octal_tree(4);
2270 if (G_mb.totelem > 1024) init_metaball_octal_tree(5);
2272 /* don't polygonize metaballs with too high resolution (base mball to small)
2273 * note: Eps was 0.0001f but this was giving problems for blood animation for durian, using 0.00001f */
2274 if (G_mb.metaball_tree) {
2275 if (ob->size[0] <= 0.00001f * (G_mb.metaball_tree->first->x_max - G_mb.metaball_tree->first->x_min) ||
2276 ob->size[1] <= 0.00001f * (G_mb.metaball_tree->first->y_max - G_mb.metaball_tree->first->y_min) ||
2277 ob->size[2] <= 0.00001f * (G_mb.metaball_tree->first->z_max - G_mb.metaball_tree->first->z_min))
2279 new_pgn_element(-1); /* free values created by init_meta */
2281 MEM_freeN(G_mb.mainb);
2284 free_metaball_octal_node(G_mb.metaball_tree->first);
2285 MEM_freeN(G_mb.metaball_tree);
2286 G_mb.metaball_tree = NULL;
2292 /* width is size per polygonize cube */
2293 if (G.is_rendering) width = mb->rendersize;
2295 width = mb->wiresize;
2296 if (G.moving && mb->flag == MB_UPDATE_HALFRES) width *= 2;
2298 /* nr_cubes is just for safety, minimum is totsize */
2299 nr_cubes = (int)(0.5f + totsize / width);
2302 mbproc.function = metaball;
2303 mbproc.size = width;
2304 mbproc.bounds = nr_cubes;
2305 mbproc.cubes = NULL;
2306 mbproc.delta = width / (float)(RES * RES);
2308 polygonize(&mbproc, mb);
2310 MEM_freeN(G_mb.mainb);
2312 /* free octal tree */
2313 if (G_mb.totelem > 1) {
2314 free_metaball_octal_node(G_mb.metaball_tree->first);
2315 MEM_freeN(G_mb.metaball_tree);
2316 G_mb.metaball_tree = NULL;
2320 VERTEX *ptr = mbproc.vertices.ptr;
2322 dl = MEM_callocN(sizeof(DispList), "mbaldisp");
2323 BLI_addtail(dispbase, dl);
2324 dl->type = DL_INDEX4;
2325 dl->nr = mbproc.vertices.count;
2326 dl->parts = curindex;
2328 dl->index = indices;
2331 a = mbproc.vertices.count;
2332 dl->verts = co = MEM_mallocN(sizeof(float) * 3 * a, "mballverts");
2333 dl->nors = no = MEM_mallocN(sizeof(float) * 3 * a, "mballnors");
2335 for (a = 0; a < mbproc.vertices.count; ptr++, a++, no += 3, co += 3) {
2336 copy_v3_v3(co, ptr->co);
2337 copy_v3_v3(no, ptr->no);
2341 freepolygonize(&mbproc);
2344 /* basic vertex data functions */
2345 int BKE_mball_minmax(MetaBall *mb, float min[3], float max[3])
2349 INIT_MINMAX(min, max);
2351 for (ml = mb->elems.first; ml; ml = ml->next) {
2352 minmax_v3v3_v3(min, max, &ml->x);
2355 return (mb->elems.first != NULL);
2358 int BKE_mball_center_median(MetaBall *mb, float r_cent[3])
2365 for (ml = mb->elems.first; ml; ml = ml->next) {
2366 add_v3_v3(r_cent, &ml->x);
2370 mul_v3_fl(r_cent, 1.0f / (float)total);
2373 return (total != 0);
2376 int BKE_mball_center_bounds(MetaBall *mb, float r_cent[3])
2378 float min[3], max[3];
2380 if (BKE_mball_minmax(mb, min, max)) {
2381 mid_v3_v3v3(r_cent, min, max);
2388 void BKE_mball_translate(MetaBall *mb, float offset[3])
2392 for (ml = mb->elems.first; ml; ml = ml->next) {
2393 add_v3_v3(&ml->x, offset);