source/blender/blenkernel/intern/mball.c

   1 /*
   2  * ***** BEGIN GPL LICENSE BLOCK *****
   3  *
   4  * This program is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU General Public License
   6  * as published by the Free Software Foundation; either version 2
   7  * of the License, or (at your option) any later version.
   8  *
   9  * This program is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12  * GNU General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License
  15  * along with this program; if not, write to the Free Software Foundation,
  16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  17  *
  18  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
  19  * All rights reserved.
  20  *
  21  * Contributor(s): Jiri Hnidek <jiri.hnidek@vslib.cz>.
  22  *
  23  * ***** END GPL LICENSE BLOCK *****
  24  *
  25  * MetaBalls are created from a single Object (with a name without number in it),
  26  * here the DispList and BoundBox also is located.
  27  * All objects with the same name (but with a number in it) are added to this.
  28  *
  29  * texture coordinates are patched within the displist
  30  */
  31
  32 /** \file blender/blenkernel/intern/mball.c
  33  *  \ingroup bke
  34  */
  35
  36 #include <stdio.h>
  37 #include <string.h>
  38 #include <math.h>
  39 #include <stdlib.h>
  40 #include <ctype.h>
  41 #include <float.h>
  42 #include "MEM_guardedalloc.h"
  43
  44 #include "DNA_material_types.h"
  45 #include "DNA_object_types.h"
  46 #include "DNA_meta_types.h"
  47 #include "DNA_scene_types.h"
  48
  49
  50 #include "BLI_blenlib.h"
  51 #include "BLI_math.h"
  52 #include "BLI_utildefines.h"
  53 #include "BLI_bpath.h"
  54
  55
  56 #include "BKE_global.h"
  57 #include "BKE_main.h"
  58
  59 /*  #include "BKE_object.h" */
  60 #include "BKE_animsys.h"
  61 #include "BKE_scene.h"
  62 #include "BKE_library.h"
  63 #include "BKE_displist.h"
  64 #include "BKE_mball.h"
  65 #include "BKE_object.h"
  66 #include "BKE_material.h"
  67
  68 /* Data types */
  69
  70 typedef struct vertex {         /* surface vertex */
  71         float co[3];  /* position and surface normal */
  72         float no[3];
  73 } VERTEX;
  74
  75 typedef struct vertices {       /* list of vertices in polygonization */
  76         int count, max;             /* # vertices, max # allowed */
  77         VERTEX *ptr;                /* dynamically allocated */
  78 } VERTICES;
  79
  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 */
  83         struct corner *next;
  84 } CORNER;
  85
  86 typedef struct cube {           /* partitioning cell (cube) */
  87         int i, j, k;                /* lattice location of cube */
  88         CORNER *corners[8];         /* eight corners */
  89 } CUBE;
  90
  91 typedef struct cubes {          /* linked list of cubes acting as stack */
  92         CUBE cube;                  /* a single cube */
  93         struct cubes *next;         /* remaining elements */
  94 } CUBES;
  95
  96 typedef struct centerlist {     /* list of cube locations */
  97         int i, j, k;                /* cube location */
  98         struct centerlist *next;    /* remaining elements */
  99 } CENTERLIST;
 100
 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 */
 105 } EDGELIST;
 106
 107 typedef struct intlist {        /* list of integers */
 108         int i;                      /* an integer */
 109         struct intlist *next;       /* remaining elements */
 110 } INTLIST;
 111
 112 typedef struct intlists {       /* list of list of integers */
 113         INTLIST *list;              /* a list of integers */
 114         struct intlists *next;      /* remaining elements */
 115 } INTLISTS;
 116
 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 */
 128 } PROCESS;
 129
 130 /* dividing scene using octal tree makes polygonisation faster */
 131 typedef struct ml_pointer {
 132         struct ml_pointer *next, *prev;
 133         struct MetaElem *ml;
 134 } ml_pointer;
 135
 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 */
 145 } octal_node;
 146
 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 */
 151 } octal_tree;
 152
 153 struct pgn_elements {
 154         struct pgn_elements *next, *prev;
 155         char *data;
 156 };
 157
 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);
 164
 165 /* Global variables */
 166 static struct {
 167         float thresh;
 168         int totelem;
 169         MetaElem **mainb;
 170         octal_tree *metaball_tree;
 171 } G_mb = {0};
 172
 173 /* Functions */
 174
 175 void BKE_mball_unlink(MetaBall *mb)
 176 {
 177         int a;
 178
 179         for (a = 0; a < mb->totcol; a++) {
 180                 if (mb->mat[a]) mb->mat[a]->id.us--;
 181                 mb->mat[a] = NULL;
 182         }
 183 }
 184
 185
 186 /* do not free mball itself */
 187 void BKE_mball_free(MetaBall *mb)
 188 {
 189         BKE_mball_unlink(mb);
 190
 191         if (mb->adt) {
 192                 BKE_free_animdata((ID *)mb);
 193                 mb->adt = NULL;
 194         }
 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);
 199 }
 200
 201 MetaBall *BKE_mball_add(const char *name)
 202 {
 203         MetaBall *mb;
 204
 205         mb = BKE_libblock_alloc(&G.main->mball, ID_MB, name);
 206
 207         mb->size[0] = mb->size[1] = mb->size[2] = 1.0;
 208         mb->texflag = MB_AUTOSPACE;
 209
 210         mb->wiresize = 0.4f;
 211         mb->rendersize = 0.2f;
 212         mb->thresh = 0.6f;
 213
 214         return mb;
 215 }
 216
 217 MetaBall *BKE_mball_copy(MetaBall *mb)
 218 {
 219         MetaBall *mbn;
 220         int a;
 221
 222         mbn = BKE_libblock_copy(&mb->id);
 223
 224         BLI_duplicatelist(&mbn->elems, &mb->elems);
 225
 226         mbn->mat = MEM_dupallocN(mb->mat);
 227         for (a = 0; a < mbn->totcol; a++) {
 228                 id_us_plus((ID *)mbn->mat[a]);
 229         }
 230         mbn->bb = MEM_dupallocN(mb->bb);
 231
 232         mbn->editelems = NULL;
 233         mbn->lastelem = NULL;
 234
 235         return mbn;
 236 }
 237
 238 static void extern_local_mball(MetaBall *mb)
 239 {
 240         if (mb->mat) {
 241                 extern_local_matarar(mb->mat, mb->totcol);
 242         }
 243 }
 244
 245 void BKE_mball_make_local(MetaBall *mb)
 246 {
 247         Main *bmain = G.main;
 248         Object *ob;
 249         int is_local = FALSE, is_lib = FALSE;
 250
 251         /* - only lib users: do nothing
 252          * - only local users: set flag
 253          * - mixed: make copy
 254          */
 255
 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);
 260
 261                 return;
 262         }
 263
 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;
 268                 }
 269         }
 270
 271         if (is_local && is_lib == FALSE) {
 272                 id_clear_lib_data(bmain, &mb->id);
 273                 extern_local_mball(mb);
 274         }
 275         else if (is_local && is_lib) {
 276                 MetaBall *mb_new = BKE_mball_copy(mb);
 277                 mb_new->id.us = 0;
 278
 279                 /* Remap paths of new ID using old library as base. */
 280                 BKE_id_lib_local_paths(bmain, mb->id.lib, &mb_new->id);
 281
 282                 for (ob = G.main->object.first; ob; ob = ob->id.next) {
 283                         if (ob->data == mb) {
 284                                 if (ob->id.lib == NULL) {
 285                                         ob->data = mb_new;
 286                                         mb_new->id.us++;
 287                                         mb->id.us--;
 288                                 }
 289                         }
 290                 }
 291         }
 292 }
 293
 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)
 297 {
 298         MetaElem *ml = MEM_callocN(sizeof(MetaElem), "metaelem");
 299
 300         unit_qt(ml->quat);
 301
 302         ml->rad = 2.0;
 303         ml->s = 2.0;
 304         ml->flag = MB_SCALE_RAD;
 305
 306         switch (type) {
 307                 case MB_BALL:
 308                         ml->type = MB_BALL;
 309                         ml->expx = ml->expy = ml->expz = 1.0;
 310
 311                         break;
 312                 case MB_TUBE:
 313                         ml->type = MB_TUBE;
 314                         ml->expx = ml->expy = ml->expz = 1.0;
 315
 316                         break;
 317                 case MB_PLANE:
 318                         ml->type = MB_PLANE;
 319                         ml->expx = ml->expy = ml->expz = 1.0;
 320
 321                         break;
 322                 case MB_ELIPSOID:
 323                         ml->type = MB_ELIPSOID;
 324                         ml->expx = 1.2f;
 325                         ml->expy = 0.8f;
 326                         ml->expz = 1.0;
 327
 328                         break;
 329                 case MB_CUBE:
 330                         ml->type = MB_CUBE;
 331                         ml->expx = ml->expy = ml->expz = 1.0;
 332
 333                         break;
 334                 default:
 335                         break;
 336         }
 337
 338         BLI_addtail(&mb->elems, ml);
 339
 340         return ml;
 341 }
 342 /** Compute bounding box of all MetaElems/MetaBalls.
 343  *
 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.
 347  */
 348 void BKE_mball_texspace_calc(Object *ob)
 349 {
 350         DispList *dl;
 351         BoundBox *bb;
 352         float *data, min[3], max[3] /*, loc[3], size[3] */;
 353         int tot, do_it = FALSE;
 354
 355         if (ob->bb == NULL) ob->bb = MEM_callocN(sizeof(BoundBox), "mb boundbox");
 356         bb = ob->bb;
 357
 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;
 362
 363         dl = ob->disp.first;
 364         while (dl) {
 365                 tot = dl->nr;
 366                 if (tot) do_it = TRUE;
 367                 data = dl->verts;
 368                 while (tot--) {
 369                         /* Also weird... but longer. From utildefines. */
 370                         minmax_v3v3_v3(min, max, data);
 371                         data += 3;
 372                 }
 373                 dl = dl->next;
 374         }
 375
 376         if (!do_it) {
 377                 min[0] = min[1] = min[2] = -1.0f;
 378                 max[0] = max[1] = max[2] = 1.0f;
 379         }
 380 #if 0
 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;
 384
 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;
 388 #endif
 389         BKE_boundbox_init_from_minmax(bb, min, max);
 390 }
 391
 392 float *BKE_mball_make_orco(Object *ob, ListBase *dispbase)
 393 {
 394         BoundBox *bb;
 395         DispList *dl;
 396         float *data, *orco, *orcodata;
 397         float loc[3], size[3];
 398         int a;
 399
 400         /* restore size and loc */
 401         bb = ob->bb;
 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];
 408
 409         dl = dispbase->first;
 410         orcodata = MEM_mallocN(sizeof(float) * 3 * dl->nr, "MballOrco");
 411
 412         data = dl->verts;
 413         orco = orcodata;
 414         a = dl->nr;
 415         while (a--) {
 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];
 419
 420                 data += 3;
 421                 orco += 3;
 422         }
 423
 424         return orcodata;
 425 }
 426
 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.
 431  *
 432  * Main idiot thing here is that the system returns find_basis_mball() objects which fail a is_basis_mball() test.
 433  *
 434  * Not only that but the depsgraph and their areas depend on this behavior!, so making small fixes here isn't worth it.
 435  * - Campbell
 436  */
 437
 438
 439 /** \brief Test, if Object *ob is basic MetaBall.
 440  *
 441  * It test last character of Object ID name. If last character
 442  * is digit it return 0, else it return 1.
 443  */
 444 int BKE_mball_is_basis(Object *ob)
 445 {
 446         int len;
 447
 448         /* just a quick test */
 449         len = strlen(ob->id.name);
 450         if (isdigit(ob->id.name[len - 1]) ) return 0;
 451         return 1;
 452 }
 453
 454 /* return nonzero if ob1 is a basis mball for ob */
 455 int BKE_mball_is_basis_for(Object *ob1, Object *ob2)
 456 {
 457         int basis1nr, basis2nr;
 458         char basis1name[MAX_ID_NAME], basis2name[MAX_ID_NAME];
 459
 460         BLI_split_name_num(basis1name, &basis1nr, ob1->id.name + 2, '.');
 461         BLI_split_name_num(basis2name, &basis2nr, ob2->id.name + 2, '.');
 462
 463         if (!strcmp(basis1name, basis2name)) return BKE_mball_is_basis(ob1);
 464         else return 0;
 465 }
 466
 467 /* \brief copy some properties from object to other metaball object with same base name
 468  *
 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)
 474 {
 475         Scene *sce_iter = scene;
 476         Base *base;
 477         Object *ob;
 478         MetaBall *active_mball = (MetaBall *)active_object->data;
 479         int basisnr, obnr;
 480         char basisname[MAX_ID_NAME], obname[MAX_ID_NAME];
 481
 482         BLI_split_name_num(basisname, &basisnr, active_object->id.name + 2, '.');
 483
 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))
 486                 return;
 487
 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, '.');
 492
 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;
 497
 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;
 503                                 }
 504                         }
 505                 }
 506         }
 507 }
 508
 509 /** \brief This function finds basic MetaBall.
 510  *
 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
 514  * blended.
 515  *
 516  * warning!, is_basis_mball() can fail on returned object, see long note above.
 517  */
 518 Object *BKE_mball_basis_find(Scene *scene, Object *basis)
 519 {
 520         Scene *sce_iter = scene;
 521         Base *base;
 522         Object *ob, *bob = basis;
 523         MetaElem *ml = NULL;
 524         int basisnr, obnr;
 525         char basisname[MAX_ID_NAME], obname[MAX_ID_NAME];
 526
 527         BLI_split_name_num(basisname, &basisnr, basis->id.name + 2, '.');
 528         G_mb.totelem = 0;
 529
 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))
 532                 return NULL;
 533
 534         while (BKE_scene_base_iter_next(&sce_iter, 1, &base, &ob)) {
 535
 536                 if (ob->type == OB_MBALL) {
 537                         if (ob == bob) {
 538                                 MetaBall *mb = ob->data;
 539
 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;
 545                         }
 546                         else {
 547                                 BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
 548
 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;
 553
 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;
 559
 560                                         if (obnr < basisnr) {
 561                                                 if (!(ob->flag & OB_FROMDUPLI)) {
 562                                                         basis = ob;
 563                                                         basisnr = obnr;
 564                                                 }
 565                                         }
 566                                 }
 567                         }
 568
 569                         for ( ; ml; ml = ml->next) {
 570                                 if (!(ml->flag & MB_HIDE)) {
 571                                         G_mb.totelem++;
 572                                 }
 573                         }
 574                 }
 575         }
 576
 577         return basis;
 578 }
 579
 580
 581 /* ******************** ARITH ************************* */
 582
 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
 588  *
 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. */
 593
 594 #define RES 12 /* # converge iterations    */
 595
 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     */
 610
 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) */
 613
 614 #define HASHBIT     (5)
 615 #define HASHSIZE    (size_t)(1 << (3 * HASHBIT))   /*! < hash table size (32768) */
 616
 617 #define HASH(i, j, k) ((((( (i) & 31) << 5) | ( (j) & 31)) << 5) | ( (k) & 31) )
 618
 619 #define MB_BIT(i, bit) (((i) >> (bit)) & 1)
 620 #define FLIP(i, bit) ((i) ^ 1 << (bit)) /* flip the given bit of i */
 621
 622
 623 /* **************** POLYGONIZATION ************************ */
 624
 625 static void calc_mballco(MetaElem *ml, float vec[3])
 626 {
 627         if (ml->mat) {
 628                 mul_m4_v3((float (*)[4])ml->mat, vec);
 629         }
 630 }
 631
 632 static float densfunc(MetaElem *ball, float x, float y, float z)
 633 {
 634         float dist2;
 635         float dvec[3] = {x, y, z};
 636
 637         mul_m4_v3((float (*)[4])ball->imat, dvec);
 638
 639         switch (ball->type) {
 640                 case MB_BALL:
 641                         /* do nothing */
 642                         break;
 643                 case MB_TUBEX:
 644                         if      (dvec[0] >  ball->len) dvec[0] -= ball->len;
 645                         else if (dvec[0] < -ball->len) dvec[0] += ball->len;
 646                         else                           dvec[0] = 0.0;
 647                         break;
 648                 case MB_TUBEY:
 649                         if      (dvec[1] >  ball->len) dvec[1] -= ball->len;
 650                         else if (dvec[1] < -ball->len) dvec[1] += ball->len;
 651                         else                           dvec[1] = 0.0;
 652                         break;
 653                 case MB_TUBEZ:
 654                         if      (dvec[2] >  ball->len) dvec[2] -= ball->len;
 655                         else if (dvec[2] < -ball->len) dvec[2] += ball->len;
 656                         else                           dvec[2] = 0.0;
 657                         break;
 658                 case MB_TUBE:
 659                         if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
 660                         else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
 661                         else                            dvec[0] = 0.0;
 662                         break;
 663                 case MB_PLANE:
 664                         if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
 665                         else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
 666                         else                            dvec[0] = 0.0;
 667                         if      (dvec[1] >  ball->expy) dvec[1] -= ball->expy;
 668                         else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
 669                         else                            dvec[1] = 0.0;
 670                         break;
 671                 case MB_ELIPSOID:
 672                         dvec[0] /= ball->expx;
 673                         dvec[1] /= ball->expy;
 674                         dvec[2] /= ball->expz;
 675                         break;
 676                 case MB_CUBE:
 677                         if      (dvec[0] >  ball->expx) dvec[0] -= ball->expx;
 678                         else if (dvec[0] < -ball->expx) dvec[0] += ball->expx;
 679                         else                            dvec[0] = 0.0;
 680
 681                         if      (dvec[1] >  ball->expy) dvec[1] -= ball->expy;
 682                         else if (dvec[1] < -ball->expy) dvec[1] += ball->expy;
 683                         else                            dvec[1] = 0.0;
 684
 685                         if      (dvec[2] >  ball->expz) dvec[2] -= ball->expz;
 686                         else if (dvec[2] < -ball->expz) dvec[2] += ball->expz;
 687                         else                            dvec[2] = 0.0;
 688                         break;
 689         }
 690
 691         dist2 = 1.0f - (len_v3(dvec) / ball->rad2);
 692
 693         if ((ball->flag & MB_NEGATIVE) == 0) {
 694                 return (dist2 < 0.0f) ? -0.5f : (ball->s * dist2 * dist2 * dist2) - 0.5f;
 695         }
 696         else {
 697                 return (dist2 < 0.0f) ? 0.5f : 0.5f - (ball->s * dist2 * dist2 * dist2);
 698         }
 699 }
 700
 701 static octal_node *find_metaball_octal_node(octal_node *node, float x, float y, float z, short depth)
 702 {
 703         if (!depth) return node;
 704
 705         if (z < node->z) {
 706                 if (y < node->y) {
 707                         if (x < node->x) {
 708                                 if (node->nodes[0])
 709                                         return find_metaball_octal_node(node->nodes[0], x, y, z, depth--);
 710                                 else
 711                                         return node;
 712                         }
 713                         else {
 714                                 if (node->nodes[1])
 715                                         return find_metaball_octal_node(node->nodes[1], x, y, z, depth--);
 716                                 else
 717                                         return node;
 718                         }
 719                 }
 720                 else {
 721                         if (x < node->x) {
 722                                 if (node->nodes[3])
 723                                         return find_metaball_octal_node(node->nodes[3], x, y, z, depth--);
 724                                 else
 725                                         return node;
 726                         }
 727                         else {
 728                                 if (node->nodes[2])
 729                                         return find_metaball_octal_node(node->nodes[2], x, y, z, depth--);
 730                                 else
 731                                         return node;
 732                         }
 733                 }
 734         }
 735         else {
 736                 if (y < node->y) {
 737                         if (x < node->x) {
 738                                 if (node->nodes[4])
 739                                         return find_metaball_octal_node(node->nodes[4], x, y, z, depth--);
 740                                 else
 741                                         return node;
 742                         }
 743                         else {
 744                                 if (node->nodes[5])
 745                                         return find_metaball_octal_node(node->nodes[5], x, y, z, depth--);
 746                                 else
 747                                         return node;
 748                         }
 749                 }
 750                 else {
 751                         if (x < node->x) {
 752                                 if (node->nodes[7])
 753                                         return find_metaball_octal_node(node->nodes[7], x, y, z, depth--);
 754                                 else
 755                                         return node;
 756                         }
 757                         else {
 758                                 if (node->nodes[6])
 759                                         return find_metaball_octal_node(node->nodes[6], x, y, z, depth--);
 760                                 else
 761                                         return node;
 762                         }
 763                 }
 764         }
 765
 766         return node;
 767 }
 768
 769 static float metaball(float x, float y, float z)
 770 /*  float x, y, z; */
 771 {
 772         struct octal_node *node;
 773         struct ml_pointer *ml_p;
 774         float dens = 0;
 775         int a;
 776
 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);
 779                 if (node) {
 780                         for (ml_p = node->elems.first; ml_p; ml_p = ml_p->next) {
 781                                 dens += densfunc(ml_p->ml, x, y, z);
 782                         }
 783
 784                         dens += -0.5f * (G_mb.metaball_tree->pos - node->pos);
 785                         dens +=  0.5f * (G_mb.metaball_tree->neg - node->neg);
 786                 }
 787                 else {
 788                         for (a = 0; a < G_mb.totelem; a++) {
 789                                 dens += densfunc(G_mb.mainb[a], x, y, z);
 790                         }
 791                 }
 792         }
 793         else {
 794                 dens += densfunc(G_mb.mainb[0], x, y, z);
 795         }
 796
 797         return G_mb.thresh - dens;
 798 }
 799
 800 /* ******************************************** */
 801
 802 static int *indices = NULL;
 803 static int totindex, curindex;
 804
 805
 806 static void accum_mballfaces(int i1, int i2, int i3, int i4)
 807 {
 808         int *newi, *cur;
 809         /* static int i = 0; I would like to delete altogether, but I don't dare to, yet */
 810
 811         if (totindex == curindex) {
 812                 totindex += 256;
 813                 newi = MEM_mallocN(4 * sizeof(int) * totindex, "vertindex");
 814
 815                 if (indices) {
 816                         memcpy(newi, indices, 4 * sizeof(int) * (totindex - 256));
 817                         MEM_freeN(indices);
 818                 }
 819                 indices = newi;
 820         }
 821
 822         cur = indices + 4 * curindex;
 823
 824         /* displists now support array drawing, we treat tri's as fake quad */
 825
 826         cur[0] = i1;
 827         cur[1] = i2;
 828         cur[2] = i3;
 829         if (i4 == 0)
 830                 cur[3] = i3;
 831         else
 832                 cur[3] = i4;
 833
 834         curindex++;
 835
 836 }
 837
 838 /* ******************* MEMORY MANAGEMENT *********************** */
 839 static void *new_pgn_element(int size)
 840 {
 841         /* during polygonize 1000s of elements are allocated
 842          * and never freed in between. Freeing only done at the end.
 843          */
 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};
 848         void *adr;
 849
 850         if (size > 10000 || size == 0) {
 851                 printf("incorrect use of new_pgn_element\n");
 852         }
 853         else if (size == -1) {
 854                 cur = lb.first;
 855                 while (cur) {
 856                         MEM_freeN(cur->data);
 857                         cur = cur->next;
 858                 }
 859                 BLI_freelistN(&lb);
 860
 861                 return NULL;
 862         }
 863
 864         size = 4 * ( (size + 3) / 4);
 865
 866         if (cur) {
 867                 if (size + offs < blocksize) {
 868                         adr = (void *) (cur->data + offs);
 869                         offs += size;
 870                         return adr;
 871                 }
 872         }
 873
 874         cur = MEM_callocN(sizeof(struct pgn_elements), "newpgn");
 875         cur->data = MEM_callocN(blocksize, "newpgn");
 876         BLI_addtail(&lb, cur);
 877
 878         offs = size;
 879         return cur->data;
 880 }
 881
 882 static void freepolygonize(PROCESS *p)
 883 {
 884         MEM_freeN(p->corners);
 885         MEM_freeN(p->edges);
 886         MEM_freeN(p->centers);
 887
 888         new_pgn_element(-1);
 889
 890         if (p->vertices.ptr) MEM_freeN(p->vertices.ptr);
 891 }
 892
 893 /**** Cubical Polygonization (optional) ****/
 894
 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  */
 907
 908 static INTLISTS *cubetable[256];
 909
 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
 913 };
 914 static int corner2[12] = {
 915         LBF, LTF, LTN, LTF, RBF, RTF, RTN, RTF, RBN, RBF, RTN, RTF
 916 };
 917 static int leftface[12] = {
 918         B,  L,  L,  F,  R,  T,  N,  R,  N,  B,  T,  F
 919 };
 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
 923 };
 924 /* face on right when going corner1 to corner2 */
 925
 926
 927 /* docube: triangulate the cube directly, without decomposition */
 928
 929 static void docube(CUBE *cube, PROCESS *p, MetaBall *mb)
 930 {
 931         INTLISTS *polys;
 932         CORNER *c1, *c2;
 933         int i, index = 0, count, indexar[8];
 934
 935         for (i = 0; i < 8; i++) if (cube->corners[i]->value > 0.0f) index += (1 << i);
 936
 937         for (polys = cubetable[index]; polys; polys = polys->next) {
 938                 INTLIST *edges;
 939
 940                 count = 0;
 941
 942                 for (edges = polys->list; edges; edges = edges->next) {
 943                         c1 = cube->corners[corner1[edges->i]];
 944                         c2 = cube->corners[corner2[edges->i]];
 945
 946                         indexar[count] = vertid(c1, c2, p, mb);
 947                         count++;
 948                 }
 949                 if (count > 2) {
 950                         switch (count) {
 951                                 case 3:
 952                                         accum_mballfaces(indexar[2], indexar[1], indexar[0], 0);
 953                                         break;
 954                                 case 4:
 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]);
 957                                         break;
 958                                 case 5:
 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]);
 961
 962                                         accum_mballfaces(indexar[4], indexar[3], indexar[0], 0);
 963                                         break;
 964                                 case 6:
 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]);
 968                                         }
 969                                         else {
 970                                                 accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]);
 971                                                 accum_mballfaces(indexar[5], indexar[4], indexar[3], indexar[0]);
 972                                         }
 973                                         break;
 974                                 case 7:
 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]);
 978                                         }
 979                                         else {
 980                                                 accum_mballfaces(indexar[3], indexar[2], indexar[1], indexar[0]);
 981                                                 accum_mballfaces(indexar[5], indexar[4], indexar[3], indexar[0]);
 982                                         }
 983
 984                                         accum_mballfaces(indexar[6], indexar[5], indexar[0], 0);
 985
 986                                         break;
 987                         }
 988                 }
 989         }
 990 }
 991
 992
 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 */
 996
 997 static void testface(int i, int j, int k, CUBE *old, int bit, int c1, int c2, int c3, int c4, PROCESS *p)
 998 {
 999         CUBE newc;
1000         CUBES *oldcubes = p->cubes;
1001         CORNER *corn1, *corn2, *corn3, *corn4;
1002         int n, pos;
1003
1004         corn1 = old->corners[c1];
1005         corn2 = old->corners[c2];
1006         corn3 = old->corners[c3];
1007         corn4 = old->corners[c4];
1008
1009         pos = corn1->value > 0.0f ? 1 : 0;
1010
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;
1017
1018
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;
1022
1023         newc.i = i;
1024         newc.j = j;
1025         newc.k = k;
1026         for (n = 0; n < 8; n++) newc.corners[n] = NULL;
1027
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;
1032
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);
1041
1042         p->cubes->cube = newc;
1043 }
1044
1045 /* setcorner: return corner with the given lattice location
1046  * set (and cache) its function value */
1047
1048 static CORNER *setcorner(PROCESS *p, int i, int j, int k)
1049 {
1050         /* for speed, do corner value caching here */
1051         CORNER *c;
1052         int index;
1053
1054         /* does corner exist? */
1055         index = HASH(i, j, k);
1056         c = p->corners[index];
1057
1058         for (; c != NULL; c = c->next) {
1059                 if (c->i == i && c->j == j && c->k == k) {
1060                         return c;
1061                 }
1062         }
1063
1064         c = (CORNER *) new_pgn_element(sizeof(CORNER));
1065
1066         c->i = i;
1067         c->co[0] = ((float)i - 0.5f) * p->size;
1068         c->j = j;
1069         c->co[1] = ((float)j - 0.5f) * p->size;
1070         c->k = k;
1071         c->co[2] = ((float)k - 0.5f) * p->size;
1072         c->value = p->function(c->co[0], c->co[1], c->co[2]);
1073
1074         c->next = p->corners[index];
1075         p->corners[index] = c;
1076
1077         return c;
1078 }
1079
1080
1081 /* nextcwedge: return next clockwise edge from given edge around given face */
1082
1083 static int nextcwedge(int edge, int face)
1084 {
1085         switch (edge) {
1086                 case LB:
1087                         return (face == L) ? LF : BN;
1088                 case LT:
1089                         return (face == L) ? LN : TF;
1090                 case LN:
1091                         return (face == L) ? LB : TN;
1092                 case LF:
1093                         return (face == L) ? LT : BF;
1094                 case RB:
1095                         return (face == R) ? RN : BF;
1096                 case RT:
1097                         return (face == R) ? RF : TN;
1098                 case RN:
1099                         return (face == R) ? RT : BN;
1100                 case RF:
1101                         return (face == R) ? RB : TF;
1102                 case BN:
1103                         return (face == B) ? RB : LN;
1104                 case BF:
1105                         return (face == B) ? LB : RF;
1106                 case TN:
1107                         return (face == T) ? LT : RN;
1108                 case TF:
1109                         return (face == T) ? RT : LF;
1110         }
1111         return 0;
1112 }
1113
1114
1115 /* otherface: return face adjoining edge that is not the given face */
1116
1117 static int otherface(int edge, int face)
1118 {
1119         int other = leftface[edge];
1120         return face == other ? rightface[edge] : other;
1121 }
1122
1123
1124 /* makecubetable: create the 256 entry table for cubical polygonization */
1125
1126 static void makecubetable(void)
1127 {
1128         static int is_done = FALSE;
1129         int i, e, c, done[12], pos[8];
1130
1131         if (is_done) return;
1132         is_done = TRUE;
1133
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;
1142
1143                                 /* get face that is to right of edge from pos to neg corner: */
1144                                 int face = pos[corner1[e]] ? rightface[e] : leftface[e];
1145
1146                                 while (1) {
1147                                         edge = nextcwedge(edge, face);
1148                                         done[edge] = 1;
1149                                         if (pos[corner1[edge]] != pos[corner2[edge]]) {
1150                                                 INTLIST *tmp = ints;
1151
1152                                                 ints = (INTLIST *) MEM_callocN(sizeof(INTLIST), "mball_intlist");
1153                                                 ints->i = edge;
1154                                                 ints->next = tmp; /* add edge to head of list */
1155
1156                                                 if (edge == start) break;
1157                                                 face = otherface(edge, face);
1158                                         }
1159                                 }
1160                                 lists->list = ints; /* add ints to head of table entry */
1161                                 lists->next = cubetable[i];
1162                                 cubetable[i] = lists;
1163                         }
1164         }
1165 }
1166
1167 void BKE_mball_cubeTable_free(void)
1168 {
1169         int i;
1170         INTLISTS *lists, *nlists;
1171         INTLIST *ints, *nints;
1172
1173         for (i = 0; i < 256; i++) {
1174                 lists = cubetable[i];
1175                 while (lists) {
1176                         nlists = lists->next;
1177
1178                         ints = lists->list;
1179                         while (ints) {
1180                                 nints = ints->next;
1181                                 MEM_freeN(ints);
1182                                 ints = nints;
1183                         }
1184
1185                         MEM_freeN(lists);
1186                         lists = nlists;
1187                 }
1188                 cubetable[i] = NULL;
1189         }
1190 }
1191
1192 /**** Storage ****/
1193
1194 /* setcenter: set (i, j, k) entry of table[]
1195  * return 1 if already set; otherwise, set and return 0 */
1196
1197 static int setcenter(CENTERLIST *table[], const int i, const int j, const int k)
1198 {
1199         int index;
1200         CENTERLIST *newc, *l, *q;
1201
1202         index = HASH(i, j, k);
1203         q = table[index];
1204
1205         for (l = q; l != NULL; l = l->next) {
1206                 if (l->i == i && l->j == j && l->k == k) return 1;
1207         }
1208
1209         newc = (CENTERLIST *) new_pgn_element(sizeof(CENTERLIST));
1210         newc->i = i;
1211         newc->j = j;
1212         newc->k = k;
1213         newc->next = q;
1214         table[index] = newc;
1215
1216         return 0;
1217 }
1218
1219
1220 /* setedge: set vertex id for edge */
1221
1222 static void setedge(EDGELIST *table[],
1223                     int i1, int j1,
1224                     int k1, int i2,
1225                     int j2, int k2,
1226                     int vid)
1227 {
1228         unsigned int index;
1229         EDGELIST *newe;
1230
1231         if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
1232                 int t = i1;
1233                 i1 = i2;
1234                 i2 = t;
1235                 t = j1;
1236                 j1 = j2;
1237                 j2 = t;
1238                 t = k1;
1239                 k1 = k2;
1240                 k2 = t;
1241         }
1242         index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
1243         newe = (EDGELIST *) new_pgn_element(sizeof(EDGELIST));
1244         newe->i1 = i1;
1245         newe->j1 = j1;
1246         newe->k1 = k1;
1247         newe->i2 = i2;
1248         newe->j2 = j2;
1249         newe->k2 = k2;
1250         newe->vid = vid;
1251         newe->next = table[index];
1252         table[index] = newe;
1253 }
1254
1255
1256 /* getedge: return vertex id for edge; return -1 if not set */
1257
1258 static int getedge(EDGELIST *table[],
1259                    int i1, int j1, int k1,
1260                    int i2, int j2, int k2)
1261 {
1262         EDGELIST *q;
1263
1264         if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
1265                 int t = i1;
1266                 i1 = i2;
1267                 i2 = t;
1268                 t = j1;
1269                 j1 = j2;
1270                 j2 = t;
1271                 t = k1;
1272                 k1 = k2;
1273                 k2 = t;
1274         }
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)
1279                 {
1280                         return q->vid;
1281                 }
1282         }
1283         return -1;
1284 }
1285
1286
1287 /**** Vertices ****/
1288
1289 #undef R
1290
1291
1292
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 */
1296
1297 /* addtovertices: add v to sequence of vertices */
1298
1299 static void addtovertices(VERTICES *vertices, VERTEX v)
1300 {
1301         if (vertices->count == vertices->max) {
1302                 int i;
1303                 VERTEX *newv;
1304                 vertices->max = vertices->count == 0 ? 10 : 2 * vertices->count;
1305                 newv = (VERTEX *) MEM_callocN(vertices->max * sizeof(VERTEX), "addtovertices");
1306
1307                 for (i = 0; i < vertices->count; i++) newv[i] = vertices->ptr[i];
1308
1309                 if (vertices->ptr != NULL) MEM_freeN(vertices->ptr);
1310                 vertices->ptr = newv;
1311         }
1312         vertices->ptr[vertices->count++] = v;
1313 }
1314
1315 /* vnormal: compute unit length surface normal at point */
1316
1317 static void vnormal(const float point[3], PROCESS *p, float r_no[3])
1318 {
1319         float delta = 0.2f * p->delta;
1320         float f = p->function(point[0], point[1], point[2]);
1321
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);
1326
1327         if (0) {
1328                 float tvec[3];
1329
1330                 delta *= 2.0f;
1331
1332                 f = p->function(point[0], point[1], point[2]);
1333
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;
1337
1338                 if (normalize_v3(tvec) != 0.0f) {
1339                         add_v3_v3(r_no, tvec);
1340                         normalize_v3(r_no);
1341                 }
1342         }
1343 }
1344
1345
1346 static int vertid(const CORNER *c1, const CORNER *c2, PROCESS *p, MetaBall *mb)
1347 {
1348         VERTEX v;
1349         int vid = getedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
1350
1351         if (vid != -1) {
1352                 return vid;  /* previously computed */
1353         }
1354
1355         converge(c1->co, c2->co, c1->value, c2->value, p->function, v.co, mb, 1); /* position */
1356         vnormal(v.co, p, v.no);
1357
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);
1361
1362         return vid;
1363 }
1364
1365
1366
1367
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)
1372 {
1373         int i = 0;
1374         float pos[3], neg[3];
1375         float positive = 0.0f, negative = 0.0f;
1376         float dvec[3];
1377
1378         if (v1 < 0) {
1379                 copy_v3_v3(pos, p2);
1380                 copy_v3_v3(neg, p1);
1381                 positive = v2;
1382                 negative = v1;
1383         }
1384         else {
1385                 copy_v3_v3(pos, p1);
1386                 copy_v3_v3(neg, p2);
1387                 positive = v1;
1388                 negative = v2;
1389         }
1390
1391         sub_v3_v3v3(dvec, pos, neg);
1392
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);
1398                         p[1] = neg[1];
1399                         p[2] = neg[2];
1400                         return;
1401                 }
1402                 if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
1403                         p[0] = neg[0];
1404                         p[1] = neg[1] - negative * dvec[1] / (positive - negative);
1405                         p[2] = neg[2];
1406                         return;
1407                 }
1408                 if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
1409                         p[0] = neg[0];
1410                         p[1] = neg[1];
1411                         p[2] = neg[2] - negative * dvec[2] / (positive - negative);
1412                         return;
1413                 }
1414         }
1415
1416         if ((dvec[1] == 0.0f) && (dvec[2] == 0.0f)) {
1417                 p[1] = neg[1];
1418                 p[2] = neg[2];
1419                 while (1) {
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];
1423                 }
1424         }
1425
1426         if ((dvec[0] == 0.0f) && (dvec[2] == 0.0f)) {
1427                 p[0] = neg[0];
1428                 p[2] = neg[2];
1429                 while (1) {
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];
1433                 }
1434         }
1435
1436         if ((dvec[0] == 0.0f) && (dvec[1] == 0.0f)) {
1437                 p[0] = neg[0];
1438                 p[1] = neg[1];
1439                 while (1) {
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];
1443                 }
1444         }
1445
1446         /* This is necessary to find start point */
1447         while (1) {
1448                 mid_v3_v3v3(&p[0], pos, neg);
1449
1450                 if (i++ == RES) {
1451                         return;
1452                 }
1453
1454                 if ((function(p[0], p[1], p[2])) > 0.0f) {
1455                         copy_v3_v3(pos, &p[0]);
1456                 }
1457                 else {
1458                         copy_v3_v3(neg, &p[0]);
1459                 }
1460         }
1461 }
1462
1463 /* ************************************** */
1464 static void add_cube(PROCESS *mbproc, int i, int j, int k, int count)
1465 {
1466         CUBES *ncube;
1467         int n;
1468         int a, b, c;
1469
1470         /* hmmm, not only one, but eight cube will be added on the stack
1471          * ... */
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;
1481
1482                                         ncube->cube.i = a;
1483                                         ncube->cube.j = b;
1484                                         ncube->cube.k = c;
1485
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));
1489                                 }
1490                         }
1491 }
1492
1493
1494 static void find_first_points(PROCESS *mbproc, MetaBall *mb, int a)
1495 {
1496         MetaElem *ml;
1497         float f = 0.0f;
1498
1499         ml = G_mb.mainb[a];
1500         f = 1.0 - (mb->thresh / ml->s);
1501
1502         /* Skip, when Stiffness of MetaElement is too small ... MetaElement can't be
1503          * visible alone ... but still can influence others MetaElements :-) */
1504         if (f > 0.0f) {
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*/;
1509                 float workp[3];
1510                 float dvec[3];
1511                 float tmp_v, workp_v, max_len, len, nx, ny, nz, MAXN;
1512
1513                 calc_mballco(ml, in);
1514                 in_v = mbproc->function(in[0], in[1], in[2]);
1515
1516                 for (i = 0; i < 3; i++) {
1517                         switch (ml->type) {
1518                                 case MB_BALL:
1519                                         OUT[0] = out[0] = IN[0] + index[i] * ml->rad;
1520                                         break;
1521                                 case MB_TUBE:
1522                                 case MB_PLANE:
1523                                 case MB_ELIPSOID:
1524                                 case MB_CUBE:
1525                                         OUT[0] = out[0] = IN[0] + index[i] * (ml->expx + ml->rad);
1526                                         break;
1527                         }
1528
1529                         for (j = 0; j < 3; j++) {
1530                                 switch (ml->type) {
1531                                         case MB_BALL:
1532                                                 OUT[1] = out[1] = IN[1] + index[j] * ml->rad;
1533                                                 break;
1534                                         case MB_TUBE:
1535                                         case MB_PLANE:
1536                                         case MB_ELIPSOID:
1537                                         case MB_CUBE:
1538                                                 OUT[1] = out[1] = IN[1] + index[j] * (ml->expy + ml->rad);
1539                                                 break;
1540                                 }
1541
1542                                 for (k = 0; k < 3; k++) {
1543                                         out[0] = OUT[0];
1544                                         out[1] = OUT[1];
1545                                         switch (ml->type) {
1546                                                 case MB_BALL:
1547                                                 case MB_TUBE:
1548                                                 case MB_PLANE:
1549                                                         out[2] = IN[2] + index[k] * ml->rad;
1550                                                         break;
1551                                                 case MB_ELIPSOID:
1552                                                 case MB_CUBE:
1553                                                         out[2] = IN[2] + index[k] * (ml->expz + ml->rad);
1554                                                         break;
1555                                         }
1556
1557                                         calc_mballco(ml, out);
1558
1559                                         /*out_v = mbproc->function(out[0], out[1], out[2]);*/ /*UNUSED*/
1560
1561                                         /* find "first points" on Implicit Surface of MetaElemnt ml */
1562                                         copy_v3_v3(workp, in);
1563                                         workp_v = in_v;
1564                                         max_len = len_v3v3(out, in);
1565
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);
1569
1570                                         MAXN = MAX3(nx, ny, nz);
1571                                         if (MAXN != 0.0f) {
1572                                                 dvec[0] = (out[0] - in[0]) / MAXN;
1573                                                 dvec[1] = (out[1] - in[1]) / MAXN;
1574                                                 dvec[2] = (out[2] - in[2]) / MAXN;
1575
1576                                                 len = 0.0;
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)) {
1585
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);
1590
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);
1595                                                                 }
1596                                                                 else {
1597                                                                         add_cube(mbproc, c_i, c_j, c_k, 2);
1598                                                                 }
1599                                                         }
1600                                                         len = len_v3v3(workp, in);
1601                                                         workp_v = tmp_v;
1602
1603                                                 }
1604                                         }
1605                                 }
1606                         }
1607                 }
1608         }
1609 }
1610
1611 static void polygonize(PROCESS *mbproc, MetaBall *mb)
1612 {
1613         CUBE c;
1614         int a;
1615
1616         mbproc->vertices.count = mbproc->vertices.max = 0;
1617         mbproc->vertices.ptr = NULL;
1618
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");
1623         makecubetable();
1624
1625         for (a = 0; a < G_mb.totelem; a++) {
1626
1627                 /* try to find 8 points on the surface for each MetaElem */
1628                 find_first_points(mbproc, mb, a);
1629         }
1630
1631         /* polygonize all MetaElems of current MetaBall */
1632         while (mbproc->cubes != NULL) { /* process active cubes till none left */
1633                 c = mbproc->cubes->cube;
1634
1635                 /* polygonize the cube directly: */
1636                 docube(&c, mbproc, mb);
1637
1638                 /* pop current cube from stack */
1639                 mbproc->cubes = mbproc->cubes->next;
1640
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);
1648         }
1649 }
1650
1651 static float init_meta(Scene *scene, Object *ob)    /* return totsize */
1652 {
1653         Scene *sce_iter = scene;
1654         Base *base;
1655         Object *bob;
1656         MetaBall *mb;
1657         MetaElem *ml;
1658         float size, totsize, obinv[4][4], obmat[4][4], vec[3];
1659         //float max = 0.0f;
1660         int a, obnr, zero_size = 0;
1661         char obname[MAX_ID_NAME];
1662
1663         copy_m4_m4(obmat, ob->obmat);   /* to cope with duplicators from BKE_scene_base_iter_next */
1664         invert_m4_m4(obinv, ob->obmat);
1665         a = 0;
1666
1667         BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
1668
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)) {
1672
1673                 if (bob->type == OB_MBALL) {
1674                         zero_size = 0;
1675                         ml = NULL;
1676
1677                         if (bob == ob && (base->flag & OB_FROMDUPLI) == 0) {
1678                                 mb = ob->data;
1679
1680                                 if (mb->editelems) ml = mb->editelems->first;
1681                                 else ml = mb->elems.first;
1682                         }
1683                         else {
1684                                 char name[MAX_ID_NAME];
1685                                 int nr;
1686
1687                                 BLI_split_name_num(name, &nr, bob->id.name + 2, '.');
1688                                 if (strcmp(obname, name) == 0) {
1689                                         mb = bob->data;
1690
1691                                         if (mb->editelems) ml = mb->editelems->first;
1692                                         else ml = mb->elems.first;
1693                                 }
1694                         }
1695
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) {
1699                                 zero_size = 1;
1700                         }
1701                         else if (bob->parent) {
1702                                 struct Object *pob = bob->parent;
1703                                 while (pob) {
1704                                         if (pob->size[0] == 0.0f || pob->size[1] == 0.0f || pob->size[2] == 0.0f) {
1705                                                 zero_size = 1;
1706                                                 break;
1707                                         }
1708                                         pob = pob->parent;
1709                                 }
1710                         }
1711
1712                         if (zero_size) {
1713                                 unsigned int ml_count = 0;
1714                                 while (ml) {
1715                                         ml_count++;
1716                                         ml = ml->next;
1717                                 }
1718                                 G_mb.totelem -= ml_count;
1719                         }
1720                         else {
1721                                 while (ml) {
1722                                         if (!(ml->flag & MB_HIDE)) {
1723                                                 int i;
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;
1727
1728                                                 max_x = max_y = max_z = -3.4e38;
1729                                                 min_x = min_y = min_z =  3.4e38;
1730
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;
1734
1735                                                 /* Rotation of MetaElem is stored in quat */
1736                                                 quat_to_mat4(temp3, ml->quat);
1737
1738                                                 /* Translation of MetaElem */
1739                                                 unit_m4(temp2);
1740                                                 temp2[3][0] = ml->x;
1741                                                 temp2[3][1] = ml->y;
1742                                                 temp2[3][2] = ml->z;
1743
1744                                                 mult_m4_m4m4(temp1, temp2, temp3);
1745
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));
1750
1751                                                 mat = new_pgn_element(4 * 4 * sizeof(float));
1752                                                 imat = new_pgn_element(4 * 4 * sizeof(float));
1753
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);
1759
1760                                                 invert_m4_m4(imat, mat);
1761
1762                                                 G_mb.mainb[a]->rad2 = ml->rad * ml->rad;
1763
1764                                                 G_mb.mainb[a]->mat = (float *) mat;
1765                                                 G_mb.mainb[a]->imat = (float *) imat;
1766
1767                                                 /* untransformed Bounding Box of MetaElem */
1768                                                 /* 0 */
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;
1772                                                 /* 1 */
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;
1776                                                 /* 2 */
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;
1780                                                 /* 3 */
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;
1784                                                 /* 4 */
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;
1788                                                 /* 5 */
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;
1792                                                 /* 6 */
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;
1796                                                 /* 7 */
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;
1800
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]);
1804
1805                                                 /* find max and min of transformed bb */
1806                                                 for (i = 0; i < 8; i++) {
1807                                                         /* find maximums */
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];
1811                                                         /* find  minimums */
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];
1815                                                 }
1816
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;
1822
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;
1826
1827                                                 a++;
1828                                         }
1829                                         ml = ml->next;
1830                                 }
1831                         }
1832                 }
1833         }
1834
1835
1836         /* totsize (= 'manhattan' radius) */
1837         totsize = 0.0;
1838         for (a = 0; a < G_mb.totelem; a++) {
1839
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;
1843
1844                 calc_mballco(G_mb.mainb[a], vec);
1845
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;
1852
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;
1856
1857                 calc_mballco(G_mb.mainb[a], vec);
1858
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;
1865         }
1866
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);
1869         }
1870
1871         return totsize;
1872 }
1873
1874 /* if MetaElem lies in node, then node includes MetaElem pointer (ml_p)
1875  * pointing at MetaElem (ml)
1876  */
1877 static void fill_metaball_octal_node(octal_node *node, MetaElem *ml, short i)
1878 {
1879         ml_pointer *ml_p;
1880
1881         ml_p = MEM_mallocN(sizeof(ml_pointer), "ml_pointer");
1882         ml_p->ml = ml;
1883         BLI_addtail(&(node->nodes[i]->elems), ml_p);
1884         node->count++;
1885
1886         if ((ml->flag & MB_NEGATIVE) == 0) {
1887                 node->nodes[i]->pos++;
1888         }
1889         else {
1890                 node->nodes[i]->neg++;
1891         }
1892 }
1893
1894 /* Node is subdivided as is illustrated on the following figure:
1895  *
1896  *      +------+------+
1897  *     /      /      /|
1898  *    +------+------+ |
1899  *   /      /      /| +
1900  *  +------+------+ |/|
1901  *  |      |      | + |
1902  *  |      |      |/| +
1903  *  +------+------+ |/
1904  *  |      |      | +
1905  *  |      |      |/
1906  *  +------+------+
1907  *
1908  */
1909 static void subdivide_metaball_octal_node(octal_node *node, float size_x, float size_y, float size_z, short depth)
1910 {
1911         MetaElem *ml;
1912         ml_pointer *ml_p;
1913         float x, y, z;
1914         int a, i;
1915
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;
1927         }
1928
1929         size_x /= 2;
1930         size_y /= 2;
1931         size_z /= 2;
1932
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;
1937
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;
1945
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;
1952
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;
1959
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;
1966
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;
1973
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;
1980
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;
1987
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;
1994
1995         ml_p = node->elems.first;
1996
1997         /* setting up references of MetaElems for new nodes */
1998         while (ml_p) {
1999                 ml = ml_p->ml;
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);
2006
2007                                         /* ml belongs to the (3)4th node */
2008                                         if (ml->bb->vec[6][1] >= y) {
2009                                                 fill_metaball_octal_node(node, ml, 3);
2010
2011                                                 /* ml belongs to the (7)8th node */
2012                                                 if (ml->bb->vec[6][2] >= z) {
2013                                                         fill_metaball_octal_node(node, ml, 7);
2014                                                 }
2015                                         }
2016
2017                                         /* ml belongs to the (1)2nd node */
2018                                         if (ml->bb->vec[6][0] >= x) {
2019                                                 fill_metaball_octal_node(node, ml, 1);
2020
2021                                                 /* ml belongs to the (5)6th node */
2022                                                 if (ml->bb->vec[6][2] >= z) {
2023                                                         fill_metaball_octal_node(node, ml, 5);
2024                                                 }
2025                                         }
2026
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);
2030
2031                                                 /* ml belong to the (6)7th node */
2032                                                 if (ml->bb->vec[6][2] >= z) {
2033                                                         fill_metaball_octal_node(node, ml, 6);
2034                                                 }
2035
2036                                         }
2037
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);
2041                                         }
2042
2043
2044
2045                                 }
2046                                 /* vec[0][0] is in the (1)second octant */
2047                                 else {
2048                                         /* ml belong to the (1)2nd node */
2049                                         fill_metaball_octal_node(node, ml, 1);
2050
2051                                         /* ml belongs to the (2)3th node */
2052                                         if (ml->bb->vec[6][1] >= y) {
2053                                                 fill_metaball_octal_node(node, ml, 2);
2054
2055                                                 /* ml belongs to the (6)7th node */
2056                                                 if (ml->bb->vec[6][2] >= z) {
2057                                                         fill_metaball_octal_node(node, ml, 6);
2058                                                 }
2059
2060                                         }
2061
2062                                         /* ml belongs to the (5)6th node */
2063                                         if (ml->bb->vec[6][2] >= z) {
2064                                                 fill_metaball_octal_node(node, ml, 5);
2065                                         }
2066                                 }
2067                         }
2068                         else {
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);
2073
2074                                         /* ml belongs to the (7)8th node */
2075                                         if (ml->bb->vec[6][2] >= z) {
2076                                                 fill_metaball_octal_node(node, ml, 7);
2077                                         }
2078
2079
2080                                         /* ml belongs to the (2)3th node */
2081                                         if (ml->bb->vec[6][0] >= x) {
2082                                                 fill_metaball_octal_node(node, ml, 2);
2083
2084                                                 /* ml belongs to the (6)7th node */
2085                                                 if (ml->bb->vec[6][2] >= z) {
2086                                                         fill_metaball_octal_node(node, ml, 6);
2087                                                 }
2088                                         }
2089                                 }
2090
2091                         }
2092
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);
2097
2098                                 /* ml belongs to the (6)7th node */
2099                                 if (ml->bb->vec[6][2] >= z) {
2100                                         fill_metaball_octal_node(node, ml, 6);
2101                                 }
2102                         }
2103                 }
2104                 else {
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);
2110
2111                                         if (ml->bb->vec[6][0] >= x) {
2112                                                 fill_metaball_octal_node(node, ml, 5);
2113                                         }
2114
2115                                         if (ml->bb->vec[6][1] >= y) {
2116                                                 fill_metaball_octal_node(node, ml, 7);
2117                                         }
2118
2119                                         if ((ml->bb->vec[6][0] >= x) && (ml->bb->vec[6][1] >= y)) {
2120                                                 fill_metaball_octal_node(node, ml, 6);
2121                                         }
2122                                 }
2123                                 /* vec[0][0] lies in (5)6th octant */
2124                                 else {
2125                                         fill_metaball_octal_node(node, ml, 5);
2126
2127                                         if (ml->bb->vec[6][1] >= y) {
2128                                                 fill_metaball_octal_node(node, ml, 6);
2129                                         }
2130                                 }
2131                         }
2132                         else {
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);
2136
2137                                         if (ml->bb->vec[6][0] >= x) {
2138                                                 fill_metaball_octal_node(node, ml, 6);
2139                                         }
2140                                 }
2141
2142                         }
2143
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);
2147                         }
2148                 }
2149                 ml_p = ml_p->next;
2150         }
2151
2152         /* free references of MetaElems for curent node (it is not needed anymore) */
2153         BLI_freelistN(&node->elems);
2154
2155         depth--;
2156
2157         if (depth > 0) {
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);
2161                 }
2162         }
2163 }
2164
2165 /* free all octal nodes recursively */
2166 static void free_metaball_octal_node(octal_node *node)
2167 {
2168         int a;
2169         for (a = 0; a < 8; a++) {
2170                 if (node->nodes[a] != NULL) free_metaball_octal_node(node->nodes[a]);
2171         }
2172         BLI_freelistN(&node->elems);
2173         MEM_freeN(node);
2174 }
2175
2176 /* If scene include more then one MetaElem, then octree is used */
2177 static void init_metaball_octal_tree(int depth)
2178 {
2179         struct octal_node *node;
2180         ml_pointer *ml_p;
2181         float size[3];
2182         int a;
2183
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;
2188
2189         G_mb.metaball_tree->neg = node->neg = 0;
2190         G_mb.metaball_tree->pos = node->pos = 0;
2191
2192         node->elems.first = NULL;
2193         node->elems.last = NULL;
2194         node->count = 0;
2195
2196         for (a = 0; a < 8; a++)
2197                 node->nodes[a] = NULL;
2198
2199         node->x_min = node->y_min = node->z_min = FLT_MAX;
2200         node->x_max = node->y_max = node->z_max = -FLT_MAX;
2201
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];
2207
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];
2211
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);
2215
2216                 if ((G_mb.mainb[a]->flag & MB_NEGATIVE) == 0) {
2217                         /* number of positive MetaElem in scene */
2218                         G_mb.metaball_tree->pos++;
2219                 }
2220                 else {
2221                         /* number of negative MetaElem in scene */
2222                         G_mb.metaball_tree->neg++;
2223                 }
2224         }
2225
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;
2230
2231         /* first node is subdivided recursively */
2232         subdivide_metaball_octal_node(node, size[0], size[1], size[2], G_mb.metaball_tree->depth);
2233 }
2234
2235 void BKE_mball_polygonize(Scene *scene, Object *ob, ListBase *dispbase)
2236 {
2237         PROCESS mbproc;
2238         MetaBall *mb;
2239         DispList *dl;
2240         int a, nr_cubes;
2241         float *co, *no, totsize, width;
2242
2243         mb = ob->data;
2244
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;
2248
2249         curindex = totindex = 0;
2250         indices = NULL;
2251         G_mb.thresh = mb->thresh;
2252
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");
2255
2256         /* initialize all mainb (MetaElems) */
2257         totsize = init_meta(scene, ob);
2258
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;
2263         }
2264
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);
2271
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))
2278                 {
2279                         new_pgn_element(-1); /* free values created by init_meta */
2280
2281                         MEM_freeN(G_mb.mainb);
2282
2283                         /* free tree */
2284                         free_metaball_octal_node(G_mb.metaball_tree->first);
2285                         MEM_freeN(G_mb.metaball_tree);
2286                         G_mb.metaball_tree = NULL;
2287
2288                         return;
2289                 }
2290         }
2291
2292         /* width is size per polygonize cube */
2293         if (G.is_rendering) width = mb->rendersize;
2294         else {
2295                 width = mb->wiresize;
2296                 if (G.moving && mb->flag == MB_UPDATE_HALFRES) width *= 2;
2297         }
2298         /* nr_cubes is just for safety, minimum is totsize */
2299         nr_cubes = (int)(0.5f + totsize / width);
2300
2301         /* init process */
2302         mbproc.function = metaball;
2303         mbproc.size = width;
2304         mbproc.bounds = nr_cubes;
2305         mbproc.cubes = NULL;
2306         mbproc.delta = width / (float)(RES * RES);
2307
2308         polygonize(&mbproc, mb);
2309
2310         MEM_freeN(G_mb.mainb);
2311
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;
2317         }
2318
2319         if (curindex) {
2320                 VERTEX *ptr = mbproc.vertices.ptr;
2321
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;
2327
2328                 dl->index = indices;
2329                 indices = NULL;
2330
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");
2334
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);
2338                 }
2339         }
2340
2341         freepolygonize(&mbproc);
2342 }
2343
2344 /* basic vertex data functions */
2345 int BKE_mball_minmax(MetaBall *mb, float min[3], float max[3])
2346 {
2347         MetaElem *ml;
2348
2349         INIT_MINMAX(min, max);
2350
2351         for (ml = mb->elems.first; ml; ml = ml->next) {
2352                 minmax_v3v3_v3(min, max, &ml->x);
2353         }
2354
2355         return (mb->elems.first != NULL);
2356 }
2357
2358 int BKE_mball_center_median(MetaBall *mb, float r_cent[3])
2359 {
2360         MetaElem *ml;
2361         int total = 0;
2362
2363         zero_v3(r_cent);
2364
2365         for (ml = mb->elems.first; ml; ml = ml->next) {
2366                 add_v3_v3(r_cent, &ml->x);
2367         }
2368
2369         if (total) {
2370                 mul_v3_fl(r_cent, 1.0f / (float)total);
2371         }
2372
2373         return (total != 0);
2374 }
2375
2376 int BKE_mball_center_bounds(MetaBall *mb, float r_cent[3])
2377 {
2378         float min[3], max[3];
2379
2380         if (BKE_mball_minmax(mb, min, max)) {
2381                 mid_v3_v3v3(r_cent, min, max);
2382                 return 1;
2383         }
2384
2385         return 0;
2386 }
2387
2388 void BKE_mball_translate(MetaBall *mb, float offset[3])
2389 {
2390         MetaElem *ml;
2391
2392         for (ml = mb->elems.first; ml; ml = ml->next) {
2393                 add_v3_v3(&ml->x, offset);
2394         }
2395 }