Merge branch 'blender2.7'
[blender.git] / source / blender / blenkernel / intern / mball_tessellate.c
1 /*
2  * This program is free software; you can redistribute it and/or
3  * modify it under the terms of the GNU General Public License
4  * as published by the Free Software Foundation; either version 2
5  * of the License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program; if not, write to the Free Software Foundation,
14  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
15  *
16  * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV.
17  * All rights reserved.
18  */
19
20 /** \file
21  * \ingroup bke
22  */
23
24 #include <stdio.h>
25 #include <string.h>
26 #include <math.h>
27 #include <stdlib.h>
28 #include <ctype.h>
29 #include <float.h>
30
31 #include "MEM_guardedalloc.h"
32
33 #include "DNA_object_types.h"
34 #include "DNA_meta_types.h"
35 #include "DNA_scene_types.h"
36
37 #include "BLI_listbase.h"
38 #include "BLI_math.h"
39 #include "BLI_string_utils.h"
40 #include "BLI_utildefines.h"
41 #include "BLI_memarena.h"
42
43 #include "BKE_global.h"
44
45 #include "BKE_displist.h"
46 #include "BKE_mball_tessellate.h" /* own include */
47 #include "BKE_scene.h"
48
49 #include "DEG_depsgraph.h"
50 #include "DEG_depsgraph_query.h"
51
52 #include "BLI_strict_flags.h"
53
54 /* experimental (faster) normal calculation */
55 // #define USE_ACCUM_NORMAL
56
57 /* Data types */
58
59 typedef struct corner { /* corner of a cube */
60   int i, j, k;          /* (i, j, k) is index within lattice */
61   float co[3], value;   /* location and function value */
62   struct corner *next;
63 } CORNER;
64
65 typedef struct cube { /* partitioning cell (cube) */
66   int i, j, k;        /* lattice location of cube */
67   CORNER *corners[8]; /* eight corners */
68 } CUBE;
69
70 typedef struct cubes { /* linked list of cubes acting as stack */
71   CUBE cube;           /* a single cube */
72   struct cubes *next;  /* remaining elements */
73 } CUBES;
74
75 typedef struct centerlist { /* list of cube locations */
76   int i, j, k;              /* cube location */
77   struct centerlist *next;  /* remaining elements */
78 } CENTERLIST;
79
80 typedef struct edgelist {     /* list of edges */
81   int i1, j1, k1, i2, j2, k2; /* edge corner ids */
82   int vid;                    /* vertex id */
83   struct edgelist *next;      /* remaining elements */
84 } EDGELIST;
85
86 typedef struct intlist { /* list of integers */
87   int i;                 /* an integer */
88   struct intlist *next;  /* remaining elements */
89 } INTLIST;
90
91 typedef struct intlists { /* list of list of integers */
92   INTLIST *list;          /* a list of integers */
93   struct intlists *next;  /* remaining elements */
94 } INTLISTS;
95
96 typedef struct Box { /* an AABB with pointer to metalelem */
97   float min[3], max[3];
98   const MetaElem *ml;
99 } Box;
100
101 typedef struct MetaballBVHNode { /* BVH node */
102   Box bb[2];                     /* AABB of children */
103   struct MetaballBVHNode *child[2];
104 } MetaballBVHNode;
105
106 typedef struct process {     /* parameters, storage */
107   float thresh, size;        /* mball threshold, single cube size */
108   float delta;               /* small delta for calculating normals */
109   unsigned int converge_res; /* converge procedure resolution (more = slower) */
110
111   MetaElem **mainb;          /* array of all metaelems */
112   unsigned int totelem, mem; /* number of metaelems */
113
114   MetaballBVHNode metaball_bvh; /* The simplest bvh */
115   Box allbb;                    /* Bounding box of all metaelems */
116
117   MetaballBVHNode **bvh_queue; /* Queue used during bvh traversal */
118   unsigned int bvh_queue_size;
119
120   CUBES *cubes;         /* stack of cubes waiting for polygonization */
121   CENTERLIST **centers; /* cube center hash table */
122   CORNER **corners;     /* corner value hash table */
123   EDGELIST **edges;     /* edge and vertex id hash table */
124
125   int (*indices)[4];     /* output indices */
126   unsigned int totindex; /* size of memory allocated for indices */
127   unsigned int curindex; /* number of currently added indices */
128
129   float (*co)[3], (*no)[3]; /* surface vertices - positions and normals */
130   unsigned int totvertex;   /* memory size */
131   unsigned int curvertex;   /* currently added vertices */
132
133   /* memory allocation from common pool */
134   MemArena *pgn_elements;
135 } PROCESS;
136
137 /* Forward declarations */
138 static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2);
139 static void add_cube(PROCESS *process, int i, int j, int k);
140 static void make_face(PROCESS *process, int i1, int i2, int i3, int i4);
141 static void converge(PROCESS *process, const CORNER *c1, const CORNER *c2, float r_p[3]);
142
143 /* ******************* SIMPLE BVH ********************* */
144
145 static void make_box_union(const BoundBox *a, const Box *b, Box *r_out)
146 {
147   r_out->min[0] = min_ff(a->vec[0][0], b->min[0]);
148   r_out->min[1] = min_ff(a->vec[0][1], b->min[1]);
149   r_out->min[2] = min_ff(a->vec[0][2], b->min[2]);
150
151   r_out->max[0] = max_ff(a->vec[6][0], b->max[0]);
152   r_out->max[1] = max_ff(a->vec[6][1], b->max[1]);
153   r_out->max[2] = max_ff(a->vec[6][2], b->max[2]);
154 }
155
156 static void make_box_from_metaelem(Box *r, const MetaElem *ml)
157 {
158   copy_v3_v3(r->max, ml->bb->vec[6]);
159   copy_v3_v3(r->min, ml->bb->vec[0]);
160   r->ml = ml;
161 }
162
163 /**
164  * Partitions part of mainb array [start, end) along axis s. Returns i,
165  * where centroids of elements in the [start, i) segment lie "on the right side" of div,
166  * and elements in the [i, end) segment lie "on the left"
167  */
168 static unsigned int partition_mainb(
169     MetaElem **mainb, unsigned int start, unsigned int end, unsigned int s, float div)
170 {
171   unsigned int i = start, j = end - 1;
172   div *= 2.0f;
173
174   while (1) {
175     while (i < j && div > (mainb[i]->bb->vec[6][s] + mainb[i]->bb->vec[0][s]))
176       i++;
177     while (j > i && div < (mainb[j]->bb->vec[6][s] + mainb[j]->bb->vec[0][s]))
178       j--;
179
180     if (i >= j)
181       break;
182
183     SWAP(MetaElem *, mainb[i], mainb[j]);
184     i++;
185     j--;
186   }
187
188   if (i == start) {
189     i++;
190   }
191
192   return i;
193 }
194
195 /**
196  * Recursively builds a BVH, dividing elements along the middle of the longest axis of allbox.
197  */
198 static void build_bvh_spatial(PROCESS *process,
199                               MetaballBVHNode *node,
200                               unsigned int start,
201                               unsigned int end,
202                               const Box *allbox)
203 {
204   unsigned int part, j, s;
205   float dim[3], div;
206
207   /* Maximum bvh queue size is number of nodes which are made, equals calls to this function. */
208   process->bvh_queue_size++;
209
210   dim[0] = allbox->max[0] - allbox->min[0];
211   dim[1] = allbox->max[1] - allbox->min[1];
212   dim[2] = allbox->max[2] - allbox->min[2];
213
214   s = 0;
215   if (dim[1] > dim[0] && dim[1] > dim[2])
216     s = 1;
217   else if (dim[2] > dim[1] && dim[2] > dim[0])
218     s = 2;
219
220   div = allbox->min[s] + (dim[s] / 2.0f);
221
222   part = partition_mainb(process->mainb, start, end, s, div);
223
224   make_box_from_metaelem(&node->bb[0], process->mainb[start]);
225   node->child[0] = NULL;
226
227   if (part > start + 1) {
228     for (j = start; j < part; j++) {
229       make_box_union(process->mainb[j]->bb, &node->bb[0], &node->bb[0]);
230     }
231
232     node->child[0] = BLI_memarena_alloc(process->pgn_elements, sizeof(MetaballBVHNode));
233     build_bvh_spatial(process, node->child[0], start, part, &node->bb[0]);
234   }
235
236   node->child[1] = NULL;
237   if (part < end) {
238     make_box_from_metaelem(&node->bb[1], process->mainb[part]);
239
240     if (part < end - 1) {
241       for (j = part; j < end; j++) {
242         make_box_union(process->mainb[j]->bb, &node->bb[1], &node->bb[1]);
243       }
244
245       node->child[1] = BLI_memarena_alloc(process->pgn_elements, sizeof(MetaballBVHNode));
246       build_bvh_spatial(process, node->child[1], part, end, &node->bb[1]);
247     }
248   }
249   else {
250     INIT_MINMAX(node->bb[1].min, node->bb[1].max);
251   }
252 }
253
254 /* ******************** ARITH ************************* */
255
256 /**
257  * BASED AT CODE (but mostly rewritten) :
258  * C code from the article
259  * "An Implicit Surface Polygonizer"
260  * by Jules Bloomenthal, jbloom@beauty.gmu.edu
261  * in "Graphics Gems IV", Academic Press, 1994
262  *
263  * Authored by Jules Bloomenthal, Xerox PARC.
264  * Copyright (c) Xerox Corporation, 1991.  All rights reserved.
265  * Permission is granted to reproduce, use and distribute this code for
266  * any and all purposes, provided that this notice appears in all copies.
267  */
268
269 #define L 0   /* left direction:   -x, -i */
270 #define R 1   /* right direction:  +x, +i */
271 #define B 2   /* bottom direction: -y, -j */
272 #define T 3   /* top direction:    +y, +j */
273 #define N 4   /* near direction:   -z, -k */
274 #define F 5   /* far direction:    +z, +k */
275 #define LBN 0 /* left bottom near corner  */
276 #define LBF 1 /* left bottom far corner   */
277 #define LTN 2 /* left top near corner     */
278 #define LTF 3 /* left top far corner      */
279 #define RBN 4 /* right bottom near corner */
280 #define RBF 5 /* right bottom far corner  */
281 #define RTN 6 /* right top near corner    */
282 #define RTF 7 /* right top far corner     */
283
284 /**
285  * the LBN corner of cube (i, j, k), corresponds with location
286  * (i-0.5)*size, (j-0.5)*size, (k-0.5)*size)
287  */
288
289 #define HASHBIT (5)
290 #define HASHSIZE (size_t)(1 << (3 * HASHBIT)) /*! < hash table size (32768) */
291
292 #define HASH(i, j, k) ((((((i)&31) << 5) | ((j)&31)) << 5) | ((k)&31))
293
294 #define MB_BIT(i, bit) (((i) >> (bit)) & 1)
295 // #define FLIP(i, bit) ((i) ^ 1 << (bit)) /* flip the given bit of i */
296
297 /* ******************** DENSITY COPMPUTATION ********************* */
298
299 /**
300  * Computes density from given metaball at given position.
301  * Metaball equation is: ``(1 - r^2 / R^2)^3 * s``
302  *
303  * r = distance from center
304  * R = metaball radius
305  * s - metaball stiffness
306  */
307 static float densfunc(const MetaElem *ball, float x, float y, float z)
308 {
309   float dist2;
310   float dvec[3] = {x, y, z};
311
312   mul_m4_v3((float(*)[4])ball->imat, dvec);
313
314   switch (ball->type) {
315     case MB_BALL:
316       /* do nothing */
317       break;
318     case MB_CUBE:
319       if (dvec[2] > ball->expz)
320         dvec[2] -= ball->expz;
321       else if (dvec[2] < -ball->expz)
322         dvec[2] += ball->expz;
323       else
324         dvec[2] = 0.0;
325       ATTR_FALLTHROUGH;
326     case MB_PLANE:
327       if (dvec[1] > ball->expy)
328         dvec[1] -= ball->expy;
329       else if (dvec[1] < -ball->expy)
330         dvec[1] += ball->expy;
331       else
332         dvec[1] = 0.0;
333       ATTR_FALLTHROUGH;
334     case MB_TUBE:
335       if (dvec[0] > ball->expx)
336         dvec[0] -= ball->expx;
337       else if (dvec[0] < -ball->expx)
338         dvec[0] += ball->expx;
339       else
340         dvec[0] = 0.0;
341       break;
342     case MB_ELIPSOID:
343       dvec[0] /= ball->expx;
344       dvec[1] /= ball->expy;
345       dvec[2] /= ball->expz;
346       break;
347
348     /* *** deprecated, could be removed?, do-versioned at least *** */
349     case MB_TUBEX:
350       if (dvec[0] > ball->len)
351         dvec[0] -= ball->len;
352       else if (dvec[0] < -ball->len)
353         dvec[0] += ball->len;
354       else
355         dvec[0] = 0.0;
356       break;
357     case MB_TUBEY:
358       if (dvec[1] > ball->len)
359         dvec[1] -= ball->len;
360       else if (dvec[1] < -ball->len)
361         dvec[1] += ball->len;
362       else
363         dvec[1] = 0.0;
364       break;
365     case MB_TUBEZ:
366       if (dvec[2] > ball->len)
367         dvec[2] -= ball->len;
368       else if (dvec[2] < -ball->len)
369         dvec[2] += ball->len;
370       else
371         dvec[2] = 0.0;
372       break;
373       /* *** end deprecated *** */
374   }
375
376   /* ball->rad2 is inverse of squared rad */
377   dist2 = 1.0f - (len_squared_v3(dvec) * ball->rad2);
378
379   /* ball->s is negative if metaball is negative */
380   return (dist2 < 0.0f) ? 0.0f : (ball->s * dist2 * dist2 * dist2);
381 }
382
383 /**
384  * Computes density at given position form all metaballs which contain this point in their box.
385  * Traverses BVH using a queue.
386  */
387 static float metaball(PROCESS *process, float x, float y, float z)
388 {
389   int i;
390   float dens = 0.0f;
391   unsigned int front = 0, back = 0;
392   MetaballBVHNode *node;
393
394   process->bvh_queue[front++] = &process->metaball_bvh;
395
396   while (front != back) {
397     node = process->bvh_queue[back++];
398
399     for (i = 0; i < 2; i++) {
400       if ((node->bb[i].min[0] <= x) && (node->bb[i].max[0] >= x) && (node->bb[i].min[1] <= y) &&
401           (node->bb[i].max[1] >= y) && (node->bb[i].min[2] <= z) && (node->bb[i].max[2] >= z)) {
402         if (node->child[i])
403           process->bvh_queue[front++] = node->child[i];
404         else
405           dens += densfunc(node->bb[i].ml, x, y, z);
406       }
407     }
408   }
409
410   return process->thresh - dens;
411 }
412
413 /**
414  * Adds face to indices, expands memory if needed.
415  */
416 static void make_face(PROCESS *process, int i1, int i2, int i3, int i4)
417 {
418   int *cur;
419
420 #ifdef USE_ACCUM_NORMAL
421   float n[3];
422 #endif
423
424   if (UNLIKELY(process->totindex == process->curindex)) {
425     process->totindex += 4096;
426     process->indices = MEM_reallocN(process->indices, sizeof(int[4]) * process->totindex);
427   }
428
429   cur = process->indices[process->curindex++];
430
431   /* displists now support array drawing, we treat tri's as fake quad */
432
433   cur[0] = i1;
434   cur[1] = i2;
435   cur[2] = i3;
436   cur[3] = i4;
437
438 #ifdef USE_ACCUM_NORMAL
439   if (i4 == i3) {
440     normal_tri_v3(n, process->co[i1], process->co[i2], process->co[i3]);
441     accumulate_vertex_normals_v3(process->no[i1],
442                                  process->no[i2],
443                                  process->no[i3],
444                                  NULL,
445                                  n,
446                                  process->co[i1],
447                                  process->co[i2],
448                                  process->co[i3],
449                                  NULL);
450   }
451   else {
452     normal_quad_v3(n, process->co[i1], process->co[i2], process->co[i3], process->co[i4]);
453     accumulate_vertex_normals_v3(process->no[i1],
454                                  process->no[i2],
455                                  process->no[i3],
456                                  process->no[i4],
457                                  n,
458                                  process->co[i1],
459                                  process->co[i2],
460                                  process->co[i3],
461                                  process->co[i4]);
462   }
463 #endif
464 }
465
466 /* Frees allocated memory */
467 static void freepolygonize(PROCESS *process)
468 {
469   if (process->corners)
470     MEM_freeN(process->corners);
471   if (process->edges)
472     MEM_freeN(process->edges);
473   if (process->centers)
474     MEM_freeN(process->centers);
475   if (process->mainb)
476     MEM_freeN(process->mainb);
477   if (process->bvh_queue)
478     MEM_freeN(process->bvh_queue);
479   if (process->pgn_elements)
480     BLI_memarena_free(process->pgn_elements);
481 }
482
483 /* **************** POLYGONIZATION ************************ */
484
485 /**** Cubical Polygonization (optional) ****/
486
487 #define LB 0  /* left bottom edge */
488 #define LT 1  /* left top edge */
489 #define LN 2  /* left near edge */
490 #define LF 3  /* left far edge */
491 #define RB 4  /* right bottom edge */
492 #define RT 5  /* right top edge */
493 #define RN 6  /* right near edge */
494 #define RF 7  /* right far edge */
495 #define BN 8  /* bottom near edge */
496 #define BF 9  /* bottom far edge */
497 #define TN 10 /* top near edge */
498 #define TF 11 /* top far edge */
499
500 static INTLISTS *cubetable[256];
501 static char faces[256];
502
503 /* edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
504 static int corner1[12] = {
505     LBN,
506     LTN,
507     LBN,
508     LBF,
509     RBN,
510     RTN,
511     RBN,
512     RBF,
513     LBN,
514     LBF,
515     LTN,
516     LTF,
517 };
518 static int corner2[12] = {
519     LBF,
520     LTF,
521     LTN,
522     LTF,
523     RBF,
524     RTF,
525     RTN,
526     RTF,
527     RBN,
528     RBF,
529     RTN,
530     RTF,
531 };
532 static int leftface[12] = {
533     B,
534     L,
535     L,
536     F,
537     R,
538     T,
539     N,
540     R,
541     N,
542     B,
543     T,
544     F,
545 };
546 /* face on left when going corner1 to corner2 */
547 static int rightface[12] = {
548     L,
549     T,
550     N,
551     L,
552     B,
553     R,
554     R,
555     F,
556     B,
557     F,
558     N,
559     T,
560 };
561 /* face on right when going corner1 to corner2 */
562
563 /**
564  * triangulate the cube directly, without decomposition
565  */
566 static void docube(PROCESS *process, CUBE *cube)
567 {
568   INTLISTS *polys;
569   CORNER *c1, *c2;
570   int i, index = 0, count, indexar[8];
571
572   /* Determine which case cube falls into. */
573   for (i = 0; i < 8; i++) {
574     if (cube->corners[i]->value > 0.0f) {
575       index += (1 << i);
576     }
577   }
578
579   /* Using faces[] table, adds neighbouring cube if surface intersects face in this direction. */
580   if (MB_BIT(faces[index], 0))
581     add_cube(process, cube->i - 1, cube->j, cube->k);
582   if (MB_BIT(faces[index], 1))
583     add_cube(process, cube->i + 1, cube->j, cube->k);
584   if (MB_BIT(faces[index], 2))
585     add_cube(process, cube->i, cube->j - 1, cube->k);
586   if (MB_BIT(faces[index], 3))
587     add_cube(process, cube->i, cube->j + 1, cube->k);
588   if (MB_BIT(faces[index], 4))
589     add_cube(process, cube->i, cube->j, cube->k - 1);
590   if (MB_BIT(faces[index], 5))
591     add_cube(process, cube->i, cube->j, cube->k + 1);
592
593   /* Using cubetable[], determines polygons for output. */
594   for (polys = cubetable[index]; polys; polys = polys->next) {
595     INTLIST *edges;
596
597     count = 0;
598     /* Sets needed vertex id's lying on the edges. */
599     for (edges = polys->list; edges; edges = edges->next) {
600       c1 = cube->corners[corner1[edges->i]];
601       c2 = cube->corners[corner2[edges->i]];
602
603       indexar[count] = vertid(process, c1, c2);
604       count++;
605     }
606
607     /* Adds faces to output. */
608     if (count > 2) {
609       switch (count) {
610         case 3:
611           make_face(process, indexar[2], indexar[1], indexar[0], indexar[0]); /* triangle */
612           break;
613         case 4:
614           make_face(process, indexar[3], indexar[2], indexar[1], indexar[0]);
615           break;
616         case 5:
617           make_face(process, indexar[3], indexar[2], indexar[1], indexar[0]);
618           make_face(process, indexar[4], indexar[3], indexar[0], indexar[0]); /* triangle */
619           break;
620         case 6:
621           make_face(process, indexar[3], indexar[2], indexar[1], indexar[0]);
622           make_face(process, indexar[5], indexar[4], indexar[3], indexar[0]);
623           break;
624         case 7:
625           make_face(process, indexar[3], indexar[2], indexar[1], indexar[0]);
626           make_face(process, indexar[5], indexar[4], indexar[3], indexar[0]);
627           make_face(process, indexar[6], indexar[5], indexar[0], indexar[0]); /* triangle */
628           break;
629       }
630     }
631   }
632 }
633
634 /**
635  * return corner with the given lattice location
636  * set (and cache) its function value
637  */
638 static CORNER *setcorner(PROCESS *process, int i, int j, int k)
639 {
640   /* for speed, do corner value caching here */
641   CORNER *c;
642   int index;
643
644   /* does corner exist? */
645   index = HASH(i, j, k);
646   c = process->corners[index];
647
648   for (; c != NULL; c = c->next) {
649     if (c->i == i && c->j == j && c->k == k) {
650       return c;
651     }
652   }
653
654   c = BLI_memarena_alloc(process->pgn_elements, sizeof(CORNER));
655
656   c->i = i;
657   c->co[0] = ((float)i - 0.5f) * process->size;
658   c->j = j;
659   c->co[1] = ((float)j - 0.5f) * process->size;
660   c->k = k;
661   c->co[2] = ((float)k - 0.5f) * process->size;
662
663   c->value = metaball(process, c->co[0], c->co[1], c->co[2]);
664
665   c->next = process->corners[index];
666   process->corners[index] = c;
667
668   return c;
669 }
670
671 /**
672  * return next clockwise edge from given edge around given face
673  */
674 static int nextcwedge(int edge, int face)
675 {
676   switch (edge) {
677     case LB:
678       return (face == L) ? LF : BN;
679     case LT:
680       return (face == L) ? LN : TF;
681     case LN:
682       return (face == L) ? LB : TN;
683     case LF:
684       return (face == L) ? LT : BF;
685     case RB:
686       return (face == R) ? RN : BF;
687     case RT:
688       return (face == R) ? RF : TN;
689     case RN:
690       return (face == R) ? RT : BN;
691     case RF:
692       return (face == R) ? RB : TF;
693     case BN:
694       return (face == B) ? RB : LN;
695     case BF:
696       return (face == B) ? LB : RF;
697     case TN:
698       return (face == T) ? LT : RN;
699     case TF:
700       return (face == T) ? RT : LF;
701   }
702   return 0;
703 }
704
705 /**
706  * \return the face adjoining edge that is not the given face
707  */
708 static int otherface(int edge, int face)
709 {
710   int other = leftface[edge];
711   return face == other ? rightface[edge] : other;
712 }
713
714 /**
715  * create the 256 entry table for cubical polygonization
716  */
717 static void makecubetable(void)
718 {
719   static bool is_done = false;
720   int i, e, c, done[12], pos[8];
721
722   if (is_done)
723     return;
724   is_done = true;
725
726   for (i = 0; i < 256; i++) {
727     for (e = 0; e < 12; e++)
728       done[e] = 0;
729     for (c = 0; c < 8; c++)
730       pos[c] = MB_BIT(i, c);
731     for (e = 0; e < 12; e++) {
732       if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
733         INTLIST *ints = NULL;
734         INTLISTS *lists = MEM_callocN(sizeof(INTLISTS), "mball_intlist");
735         int start = e, edge = e;
736
737         /* get face that is to right of edge from pos to neg corner: */
738         int face = pos[corner1[e]] ? rightface[e] : leftface[e];
739
740         while (1) {
741           edge = nextcwedge(edge, face);
742           done[edge] = 1;
743           if (pos[corner1[edge]] != pos[corner2[edge]]) {
744             INTLIST *tmp = ints;
745
746             ints = MEM_callocN(sizeof(INTLIST), "mball_intlist");
747             ints->i = edge;
748             ints->next = tmp; /* add edge to head of list */
749
750             if (edge == start)
751               break;
752             face = otherface(edge, face);
753           }
754         }
755         lists->list = ints; /* add ints to head of table entry */
756         lists->next = cubetable[i];
757         cubetable[i] = lists;
758       }
759     }
760   }
761
762   for (i = 0; i < 256; i++) {
763     INTLISTS *polys;
764     faces[i] = 0;
765     for (polys = cubetable[i]; polys; polys = polys->next) {
766       INTLIST *edges;
767
768       for (edges = polys->list; edges; edges = edges->next) {
769         if (edges->i == LB || edges->i == LT || edges->i == LN || edges->i == LF)
770           faces[i] |= 1 << L;
771         if (edges->i == RB || edges->i == RT || edges->i == RN || edges->i == RF)
772           faces[i] |= 1 << R;
773         if (edges->i == LB || edges->i == RB || edges->i == BN || edges->i == BF)
774           faces[i] |= 1 << B;
775         if (edges->i == LT || edges->i == RT || edges->i == TN || edges->i == TF)
776           faces[i] |= 1 << T;
777         if (edges->i == LN || edges->i == RN || edges->i == BN || edges->i == TN)
778           faces[i] |= 1 << N;
779         if (edges->i == LF || edges->i == RF || edges->i == BF || edges->i == TF)
780           faces[i] |= 1 << F;
781       }
782     }
783   }
784 }
785
786 void BKE_mball_cubeTable_free(void)
787 {
788   int i;
789   INTLISTS *lists, *nlists;
790   INTLIST *ints, *nints;
791
792   for (i = 0; i < 256; i++) {
793     lists = cubetable[i];
794     while (lists) {
795       nlists = lists->next;
796
797       ints = lists->list;
798       while (ints) {
799         nints = ints->next;
800         MEM_freeN(ints);
801         ints = nints;
802       }
803
804       MEM_freeN(lists);
805       lists = nlists;
806     }
807     cubetable[i] = NULL;
808   }
809 }
810
811 /**** Storage ****/
812
813 /**
814  * Inserts cube at lattice i, j, k into hash table, marking it as "done"
815  */
816 static int setcenter(PROCESS *process, CENTERLIST *table[], const int i, const int j, const int k)
817 {
818   int index;
819   CENTERLIST *newc, *l, *q;
820
821   index = HASH(i, j, k);
822   q = table[index];
823
824   for (l = q; l != NULL; l = l->next) {
825     if (l->i == i && l->j == j && l->k == k)
826       return 1;
827   }
828
829   newc = BLI_memarena_alloc(process->pgn_elements, sizeof(CENTERLIST));
830   newc->i = i;
831   newc->j = j;
832   newc->k = k;
833   newc->next = q;
834   table[index] = newc;
835
836   return 0;
837 }
838
839 /**
840  * Sets vid of vertex lying on given edge.
841  */
842 static void setedge(PROCESS *process, int i1, int j1, int k1, int i2, int j2, int k2, int vid)
843 {
844   int index;
845   EDGELIST *newe;
846
847   if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
848     int t = i1;
849     i1 = i2;
850     i2 = t;
851     t = j1;
852     j1 = j2;
853     j2 = t;
854     t = k1;
855     k1 = k2;
856     k2 = t;
857   }
858   index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
859   newe = BLI_memarena_alloc(process->pgn_elements, sizeof(EDGELIST));
860
861   newe->i1 = i1;
862   newe->j1 = j1;
863   newe->k1 = k1;
864   newe->i2 = i2;
865   newe->j2 = j2;
866   newe->k2 = k2;
867   newe->vid = vid;
868   newe->next = process->edges[index];
869   process->edges[index] = newe;
870 }
871
872 /**
873  * \return vertex id for edge; return -1 if not set
874  */
875 static int getedge(EDGELIST *table[], int i1, int j1, int k1, int i2, int j2, int k2)
876 {
877   EDGELIST *q;
878
879   if (i1 > i2 || (i1 == i2 && (j1 > j2 || (j1 == j2 && k1 > k2)))) {
880     int t = i1;
881     i1 = i2;
882     i2 = t;
883     t = j1;
884     j1 = j2;
885     j2 = t;
886     t = k1;
887     k1 = k2;
888     k2 = t;
889   }
890   q = table[HASH(i1, j1, k1) + HASH(i2, j2, k2)];
891   for (; q != NULL; q = q->next) {
892     if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 && q->i2 == i2 && q->j2 == j2 && q->k2 == k2) {
893       return q->vid;
894     }
895   }
896   return -1;
897 }
898
899 /**
900  * Adds a vertex, expands memory if needed.
901  */
902 static void addtovertices(PROCESS *process, const float v[3], const float no[3])
903 {
904   if (process->curvertex == process->totvertex) {
905     process->totvertex += 4096;
906     process->co = MEM_reallocN(process->co, process->totvertex * sizeof(float[3]));
907     process->no = MEM_reallocN(process->no, process->totvertex * sizeof(float[3]));
908   }
909
910   copy_v3_v3(process->co[process->curvertex], v);
911   copy_v3_v3(process->no[process->curvertex], no);
912
913   process->curvertex++;
914 }
915
916 #ifndef USE_ACCUM_NORMAL
917 /**
918  * Computes normal from density field at given point.
919  *
920  * \note Doesn't do normalization!
921  */
922 static void vnormal(PROCESS *process, const float point[3], float r_no[3])
923 {
924   const float delta = process->delta;
925   const float f = metaball(process, point[0], point[1], point[2]);
926
927   r_no[0] = metaball(process, point[0] + delta, point[1], point[2]) - f;
928   r_no[1] = metaball(process, point[0], point[1] + delta, point[2]) - f;
929   r_no[2] = metaball(process, point[0], point[1], point[2] + delta) - f;
930 }
931 #endif /* USE_ACCUM_NORMAL */
932
933 /**
934  * \return the id of vertex between two corners.
935  *
936  * If it wasn't previously computed, does #converge() and adds vertex to process.
937  */
938 static int vertid(PROCESS *process, const CORNER *c1, const CORNER *c2)
939 {
940   float v[3], no[3];
941   int vid = getedge(process->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
942
943   if (vid != -1)
944     return vid; /* previously computed */
945
946   converge(process, c1, c2, v); /* position */
947
948 #ifdef USE_ACCUM_NORMAL
949   zero_v3(no);
950 #else
951   vnormal(process, v, no);
952 #endif
953
954   addtovertices(process, v, no); /* save vertex */
955   vid = (int)process->curvertex - 1;
956   setedge(process, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
957
958   return vid;
959 }
960
961 /**
962  * Given two corners, computes approximation of surface intersection point between them.
963  * In case of small threshold, do bisection.
964  */
965 static void converge(PROCESS *process, const CORNER *c1, const CORNER *c2, float r_p[3])
966 {
967   float tmp, dens;
968   unsigned int i;
969   float c1_value, c1_co[3];
970   float c2_value, c2_co[3];
971
972   if (c1->value < c2->value) {
973     c1_value = c2->value;
974     copy_v3_v3(c1_co, c2->co);
975     c2_value = c1->value;
976     copy_v3_v3(c2_co, c1->co);
977   }
978   else {
979     c1_value = c1->value;
980     copy_v3_v3(c1_co, c1->co);
981     c2_value = c2->value;
982     copy_v3_v3(c2_co, c2->co);
983   }
984
985   for (i = 0; i < process->converge_res; i++) {
986     interp_v3_v3v3(r_p, c1_co, c2_co, 0.5f);
987     dens = metaball(process, r_p[0], r_p[1], r_p[2]);
988
989     if (dens > 0.0f) {
990       c1_value = dens;
991       copy_v3_v3(c1_co, r_p);
992     }
993     else {
994       c2_value = dens;
995       copy_v3_v3(c2_co, r_p);
996     }
997   }
998
999   tmp = -c1_value / (c2_value - c1_value);
1000   interp_v3_v3v3(r_p, c1_co, c2_co, tmp);
1001 }
1002
1003 /**
1004  * Adds cube at given lattice position to cube stack of process.
1005  */
1006 static void add_cube(PROCESS *process, int i, int j, int k)
1007 {
1008   CUBES *ncube;
1009   int n;
1010
1011   /* test if cube has been found before */
1012   if (setcenter(process, process->centers, i, j, k) == 0) {
1013     /* push cube on stack: */
1014     ncube = BLI_memarena_alloc(process->pgn_elements, sizeof(CUBES));
1015     ncube->next = process->cubes;
1016     process->cubes = ncube;
1017
1018     ncube->cube.i = i;
1019     ncube->cube.j = j;
1020     ncube->cube.k = k;
1021
1022     /* set corners of initial cube: */
1023     for (n = 0; n < 8; n++)
1024       ncube->cube.corners[n] = setcorner(
1025           process, i + MB_BIT(n, 2), j + MB_BIT(n, 1), k + MB_BIT(n, 0));
1026   }
1027 }
1028
1029 static void next_lattice(int r[3], const float pos[3], const float size)
1030 {
1031   r[0] = (int)ceil((pos[0] / size) + 0.5f);
1032   r[1] = (int)ceil((pos[1] / size) + 0.5f);
1033   r[2] = (int)ceil((pos[2] / size) + 0.5f);
1034 }
1035 static void prev_lattice(int r[3], const float pos[3], const float size)
1036 {
1037   next_lattice(r, pos, size);
1038   r[0]--;
1039   r[1]--;
1040   r[2]--;
1041 }
1042 static void closest_latice(int r[3], const float pos[3], const float size)
1043 {
1044   r[0] = (int)floorf(pos[0] / size + 1.0f);
1045   r[1] = (int)floorf(pos[1] / size + 1.0f);
1046   r[2] = (int)floorf(pos[2] / size + 1.0f);
1047 }
1048
1049 /**
1050  * Find at most 26 cubes to start polygonization from.
1051  */
1052 static void find_first_points(PROCESS *process, const unsigned int em)
1053 {
1054   const MetaElem *ml;
1055   int center[3], lbn[3], rtf[3], it[3], dir[3], add[3];
1056   float tmp[3], a, b;
1057
1058   ml = process->mainb[em];
1059
1060   mid_v3_v3v3(tmp, ml->bb->vec[0], ml->bb->vec[6]);
1061   closest_latice(center, tmp, process->size);
1062   prev_lattice(lbn, ml->bb->vec[0], process->size);
1063   next_lattice(rtf, ml->bb->vec[6], process->size);
1064
1065   for (dir[0] = -1; dir[0] <= 1; dir[0]++) {
1066     for (dir[1] = -1; dir[1] <= 1; dir[1]++) {
1067       for (dir[2] = -1; dir[2] <= 1; dir[2]++) {
1068         if (dir[0] == 0 && dir[1] == 0 && dir[2] == 0) {
1069           continue;
1070         }
1071
1072         copy_v3_v3_int(it, center);
1073
1074         b = setcorner(process, it[0], it[1], it[2])->value;
1075         do {
1076           it[0] += dir[0];
1077           it[1] += dir[1];
1078           it[2] += dir[2];
1079           a = b;
1080           b = setcorner(process, it[0], it[1], it[2])->value;
1081
1082           if (a * b < 0.0f) {
1083             add[0] = it[0] - dir[0];
1084             add[1] = it[1] - dir[1];
1085             add[2] = it[2] - dir[2];
1086             DO_MIN(it, add);
1087             add_cube(process, add[0], add[1], add[2]);
1088             break;
1089           }
1090         } while ((it[0] > lbn[0]) && (it[1] > lbn[1]) && (it[2] > lbn[2]) && (it[0] < rtf[0]) &&
1091                  (it[1] < rtf[1]) && (it[2] < rtf[2]));
1092       }
1093     }
1094   }
1095 }
1096
1097 /**
1098  * The main polygonization proc.
1099  * Allocates memory, makes cubetable,
1100  * finds starting surface points
1101  * and processes cubes on the stack until none left.
1102  */
1103 static void polygonize(PROCESS *process)
1104 {
1105   CUBE c;
1106   unsigned int i;
1107
1108   process->centers = MEM_callocN(HASHSIZE * sizeof(CENTERLIST *), "mbproc->centers");
1109   process->corners = MEM_callocN(HASHSIZE * sizeof(CORNER *), "mbproc->corners");
1110   process->edges = MEM_callocN(2 * HASHSIZE * sizeof(EDGELIST *), "mbproc->edges");
1111   process->bvh_queue = MEM_callocN(sizeof(MetaballBVHNode *) * process->bvh_queue_size,
1112                                    "Metaball BVH Queue");
1113
1114   makecubetable();
1115
1116   for (i = 0; i < process->totelem; i++) {
1117     find_first_points(process, i);
1118   }
1119
1120   while (process->cubes != NULL) {
1121     c = process->cubes->cube;
1122     process->cubes = process->cubes->next;
1123
1124     docube(process, &c);
1125   }
1126 }
1127
1128 /**
1129  * Iterates over ALL objects in the scene and all of its sets, including
1130  * making all duplis(not only metas). Copies metas to mainb array.
1131  * Computes bounding boxes for building BVH. */
1132 static void init_meta(Depsgraph *depsgraph, PROCESS *process, Scene *scene, Object *ob)
1133 {
1134   Scene *sce_iter = scene;
1135   Base *base;
1136   Object *bob;
1137   MetaBall *mb;
1138   const MetaElem *ml;
1139   float obinv[4][4], obmat[4][4];
1140   unsigned int i;
1141   int obnr, zero_size = 0;
1142   char obname[MAX_ID_NAME];
1143   SceneBaseIter iter;
1144
1145   copy_m4_m4(obmat, ob->obmat); /* to cope with duplicators from BKE_scene_base_iter_next */
1146   invert_m4_m4(obinv, ob->obmat);
1147
1148   BLI_split_name_num(obname, &obnr, ob->id.name + 2, '.');
1149
1150   /* make main array */
1151   BKE_scene_base_iter_next(depsgraph, &iter, &sce_iter, 0, NULL, NULL);
1152   while (BKE_scene_base_iter_next(depsgraph, &iter, &sce_iter, 1, &base, &bob)) {
1153     if (bob->type == OB_MBALL) {
1154       zero_size = 0;
1155       ml = NULL;
1156
1157       if (bob == ob && (base->flag_legacy & OB_FROMDUPLI) == 0) {
1158         mb = ob->data;
1159
1160         if (mb->editelems)
1161           ml = mb->editelems->first;
1162         else
1163           ml = mb->elems.first;
1164       }
1165       else {
1166         char name[MAX_ID_NAME];
1167         int nr;
1168
1169         BLI_split_name_num(name, &nr, bob->id.name + 2, '.');
1170         if (STREQ(obname, name)) {
1171           mb = bob->data;
1172
1173           if (mb->editelems)
1174             ml = mb->editelems->first;
1175           else
1176             ml = mb->elems.first;
1177         }
1178       }
1179
1180       /* when metaball object has zero scale, then MetaElem to this MetaBall
1181        * will not be put to mainb array */
1182       if (has_zero_axis_m4(bob->obmat)) {
1183         zero_size = 1;
1184       }
1185       else if (bob->parent) {
1186         struct Object *pob = bob->parent;
1187         while (pob) {
1188           if (has_zero_axis_m4(pob->obmat)) {
1189             zero_size = 1;
1190             break;
1191           }
1192           pob = pob->parent;
1193         }
1194       }
1195
1196       if (zero_size) {
1197         while (ml) {
1198           ml = ml->next;
1199         }
1200       }
1201       else {
1202         while (ml) {
1203           if (!(ml->flag & MB_HIDE)) {
1204             float pos[4][4], rot[4][4];
1205             float expx, expy, expz;
1206             float tempmin[3], tempmax[3];
1207
1208             MetaElem *new_ml;
1209
1210             /* make a copy because of duplicates */
1211             new_ml = BLI_memarena_alloc(process->pgn_elements, sizeof(MetaElem));
1212             *(new_ml) = *ml;
1213             new_ml->bb = BLI_memarena_alloc(process->pgn_elements, sizeof(BoundBox));
1214             new_ml->mat = BLI_memarena_alloc(process->pgn_elements, 4 * 4 * sizeof(float));
1215             new_ml->imat = BLI_memarena_alloc(process->pgn_elements, 4 * 4 * sizeof(float));
1216
1217             /* too big stiffness seems only ugly due to linear interpolation
1218              * no need to have possibility for too big stiffness */
1219             if (ml->s > 10.0f)
1220               new_ml->s = 10.0f;
1221             else
1222               new_ml->s = ml->s;
1223
1224             /* if metaball is negative, set stiffness negative */
1225             if (new_ml->flag & MB_NEGATIVE)
1226               new_ml->s = -new_ml->s;
1227
1228             /* Translation of MetaElem */
1229             unit_m4(pos);
1230             pos[3][0] = ml->x;
1231             pos[3][1] = ml->y;
1232             pos[3][2] = ml->z;
1233
1234             /* Rotation of MetaElem is stored in quat */
1235             quat_to_mat4(rot, ml->quat);
1236
1237             /* basis object space -> world -> ml object space -> position -> rotation -> ml local space */
1238             mul_m4_series((float(*)[4])new_ml->mat, obinv, bob->obmat, pos, rot);
1239             /* ml local space -> basis object space */
1240             invert_m4_m4((float(*)[4])new_ml->imat, (float(*)[4])new_ml->mat);
1241
1242             /* rad2 is inverse of squared radius */
1243             new_ml->rad2 = 1 / (ml->rad * ml->rad);
1244
1245             /* initial dimensions = radius */
1246             expx = ml->rad;
1247             expy = ml->rad;
1248             expz = ml->rad;
1249
1250             switch (ml->type) {
1251               case MB_BALL:
1252                 break;
1253               case MB_CUBE: /* cube is "expanded" by expz, expy and expx */
1254                 expz += ml->expz;
1255                 ATTR_FALLTHROUGH;
1256               case MB_PLANE: /* plane is "expanded" by expy and expx */
1257                 expy += ml->expy;
1258                 ATTR_FALLTHROUGH;
1259               case MB_TUBE: /* tube is "expanded" by expx */
1260                 expx += ml->expx;
1261                 break;
1262               case MB_ELIPSOID: /* ellipsoid is "stretched" by exp* */
1263                 expx *= ml->expx;
1264                 expy *= ml->expy;
1265                 expz *= ml->expz;
1266                 break;
1267             }
1268
1269             /* untransformed Bounding Box of MetaElem */
1270             /* TODO, its possible the elem type has been changed and the exp* values can use a fallback */
1271             copy_v3_fl3(new_ml->bb->vec[0], -expx, -expy, -expz); /* 0 */
1272             copy_v3_fl3(new_ml->bb->vec[1], +expx, -expy, -expz); /* 1 */
1273             copy_v3_fl3(new_ml->bb->vec[2], +expx, +expy, -expz); /* 2 */
1274             copy_v3_fl3(new_ml->bb->vec[3], -expx, +expy, -expz); /* 3 */
1275             copy_v3_fl3(new_ml->bb->vec[4], -expx, -expy, +expz); /* 4 */
1276             copy_v3_fl3(new_ml->bb->vec[5], +expx, -expy, +expz); /* 5 */
1277             copy_v3_fl3(new_ml->bb->vec[6], +expx, +expy, +expz); /* 6 */
1278             copy_v3_fl3(new_ml->bb->vec[7], -expx, +expy, +expz); /* 7 */
1279
1280             /* transformation of Metalem bb */
1281             for (i = 0; i < 8; i++)
1282               mul_m4_v3((float(*)[4])new_ml->mat, new_ml->bb->vec[i]);
1283
1284             /* find max and min of transformed bb */
1285             INIT_MINMAX(tempmin, tempmax);
1286             for (i = 0; i < 8; i++) {
1287               DO_MINMAX(new_ml->bb->vec[i], tempmin, tempmax);
1288             }
1289
1290             /* set only point 0 and 6 - AABB of Metaelem */
1291             copy_v3_v3(new_ml->bb->vec[0], tempmin);
1292             copy_v3_v3(new_ml->bb->vec[6], tempmax);
1293
1294             /* add new_ml to mainb[] */
1295             if (UNLIKELY(process->totelem == process->mem)) {
1296               process->mem = process->mem * 2 + 10;
1297               process->mainb = MEM_reallocN(process->mainb, sizeof(MetaElem *) * process->mem);
1298             }
1299             process->mainb[process->totelem++] = new_ml;
1300           }
1301           ml = ml->next;
1302         }
1303       }
1304     }
1305   }
1306
1307   /* compute AABB of all Metaelems */
1308   if (process->totelem > 0) {
1309     copy_v3_v3(process->allbb.min, process->mainb[0]->bb->vec[0]);
1310     copy_v3_v3(process->allbb.max, process->mainb[0]->bb->vec[6]);
1311     for (i = 1; i < process->totelem; i++)
1312       make_box_union(process->mainb[i]->bb, &process->allbb, &process->allbb);
1313   }
1314 }
1315
1316 void BKE_mball_polygonize(Depsgraph *depsgraph, Scene *scene, Object *ob, ListBase *dispbase)
1317 {
1318   MetaBall *mb;
1319   DispList *dl;
1320   unsigned int a;
1321   PROCESS process = {0};
1322   bool is_render = DEG_get_mode(depsgraph) == DAG_EVAL_RENDER;
1323
1324   mb = ob->data;
1325
1326   process.thresh = mb->thresh;
1327
1328   if (process.thresh < 0.001f)
1329     process.converge_res = 16;
1330   else if (process.thresh < 0.01f)
1331     process.converge_res = 8;
1332   else if (process.thresh < 0.1f)
1333     process.converge_res = 4;
1334   else
1335     process.converge_res = 2;
1336
1337   if (is_render && (mb->flag == MB_UPDATE_NEVER))
1338     return;
1339   if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_FAST)
1340     return;
1341
1342   if (is_render) {
1343     process.size = mb->rendersize;
1344   }
1345   else {
1346     process.size = mb->wiresize;
1347     if ((G.moving & (G_TRANSFORM_OBJ | G_TRANSFORM_EDIT)) && mb->flag == MB_UPDATE_HALFRES) {
1348       process.size *= 2.0f;
1349     }
1350   }
1351
1352   process.delta = process.size * 0.001f;
1353
1354   process.pgn_elements = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, "Metaball memarena");
1355
1356   /* initialize all mainb (MetaElems) */
1357   init_meta(depsgraph, &process, scene, ob);
1358
1359   if (process.totelem > 0) {
1360     build_bvh_spatial(&process, &process.metaball_bvh, 0, process.totelem, &process.allbb);
1361
1362     /* don't polygonize metaballs with too high resolution (base mball to small)
1363      * note: Eps was 0.0001f but this was giving problems for blood animation for durian, using 0.00001f */
1364     if (ob->scale[0] > 0.00001f * (process.allbb.max[0] - process.allbb.min[0]) ||
1365         ob->scale[1] > 0.00001f * (process.allbb.max[1] - process.allbb.min[1]) ||
1366         ob->scale[2] > 0.00001f * (process.allbb.max[2] - process.allbb.min[2])) {
1367       polygonize(&process);
1368
1369       /* add resulting surface to displist */
1370       if (process.curindex) {
1371         dl = MEM_callocN(sizeof(DispList), "mballdisp");
1372         BLI_addtail(dispbase, dl);
1373         dl->type = DL_INDEX4;
1374         dl->nr = (int)process.curvertex;
1375         dl->parts = (int)process.curindex;
1376
1377         dl->index = (int *)process.indices;
1378
1379         for (a = 0; a < process.curvertex; a++) {
1380           normalize_v3(process.no[a]);
1381         }
1382
1383         dl->verts = (float *)process.co;
1384         dl->nors = (float *)process.no;
1385       }
1386     }
1387   }
1388
1389   freepolygonize(&process);
1390 }