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