Merge branch 'blender-v2.81-release'
[blender.git] / source / blender / modifiers / intern / MOD_skin.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
17 /** \file
18  * \ingroup modifiers
19  */
20
21 /* Implementation based in part off the paper "B-Mesh: A Fast Modeling
22  * System for Base Meshes of 3D Articulated Shapes" (Zhongping Ji,
23  * Ligang Liu, Yigang Wang)
24  *
25  * Note that to avoid confusion with Blender's BMesh data structure,
26  * this tool is renamed as the Skin modifier.
27  *
28  * The B-Mesh paper is current available here:
29  * http://www.math.zju.edu.cn/ligangliu/CAGD/Projects/BMesh/
30  *
31  * The main missing features in this code compared to the paper are:
32  *
33  * + No mesh evolution. The paper suggests iteratively subsurfing the
34  *   skin output and adapting the output to better conform with the
35  *   spheres of influence surrounding each vertex.
36  *
37  * + No mesh fairing. The paper suggests re-aligning output edges to
38  *   follow principal mesh curvatures.
39  *
40  * + No auxiliary balls. These would serve to influence mesh
41  *   evolution, which as noted above is not implemented.
42  *
43  * The code also adds some features not present in the paper:
44  *
45  * + Loops in the input edge graph.
46  *
47  * + Concave surfaces around branch nodes. The paper does not discuss
48  *   how to handle non-convex regions; this code adds a number of
49  *   cleanup operations to handle many (though not all) of these
50  *   cases.
51  */
52
53 #include "MEM_guardedalloc.h"
54
55 #include "BLI_utildefines.h"
56
57 #include "BLI_array.h"
58 #include "BLI_bitmap.h"
59 #include "BLI_heap_simple.h"
60 #include "BLI_math.h"
61 #include "BLI_math_geom.h"
62 #include "BLI_stack.h"
63
64 #include "DNA_mesh_types.h"
65 #include "DNA_meshdata_types.h"
66 #include "DNA_object_types.h"
67 #include "DNA_modifier_types.h"
68
69 #include "BKE_deform.h"
70 #include "BKE_library.h"
71 #include "BKE_mesh.h"
72 #include "BKE_mesh_mapping.h"
73 #include "BKE_modifier.h"
74
75 #include "MOD_modifiertypes.h"
76
77 #include "bmesh.h"
78
79 typedef struct {
80   float mat[3][3];
81   /* Vert that edge is pointing away from, no relation to
82    * MEdge.v1 */
83   int origin;
84 } EMat;
85
86 typedef enum {
87   CAP_START = 1,
88   CAP_END = 2,
89   SEAM_FRAME = 4,
90   FLIP_NORMAL = 8,
91 } SkinNodeFlag;
92
93 typedef struct Frame {
94   /* Index in the MVert array */
95   BMVert *verts[4];
96   /* Location of each corner */
97   float co[4][3];
98   /* Indicates which corners have been merged with another
99    * frame's corner (so they share an MVert index) */
100   struct {
101     /* Merge to target frame/corner (no merge if frame is null) */
102     struct Frame *frame;
103     int corner;
104     /* checked to avoid chaining.
105      * (merging when we're already been referenced), see T39775 */
106     uint is_target : 1;
107   } merge[4];
108
109   /* For hull frames, whether each vertex is detached or not */
110   bool inside_hull[4];
111   /* Whether any part of the frame (corner or edge) is detached */
112   bool detached;
113 } Frame;
114
115 #define MAX_SKIN_NODE_FRAMES 2
116 typedef struct {
117   Frame frames[MAX_SKIN_NODE_FRAMES];
118   int totframe;
119
120   SkinNodeFlag flag;
121
122   /* Used for hulling a loop seam */
123   int seam_edges[2];
124 } SkinNode;
125
126 typedef struct {
127   BMesh *bm;
128   SkinModifierData *smd;
129   int mat_nr;
130 } SkinOutput;
131
132 static void add_poly(SkinOutput *so, BMVert *v1, BMVert *v2, BMVert *v3, BMVert *v4);
133
134 /***************************** Convex Hull ****************************/
135
136 static bool is_quad_symmetric(BMVert *quad[4], const SkinModifierData *smd)
137 {
138   const float threshold = 0.0001f;
139   const float threshold_squared = threshold * threshold;
140   int axis;
141
142   for (axis = 0; axis < 3; axis++) {
143     if (smd->symmetry_axes & (1 << axis)) {
144       float a[3];
145
146       copy_v3_v3(a, quad[0]->co);
147       a[axis] = -a[axis];
148
149       if (len_squared_v3v3(a, quad[1]->co) < threshold_squared) {
150         copy_v3_v3(a, quad[2]->co);
151         a[axis] = -a[axis];
152         if (len_squared_v3v3(a, quad[3]->co) < threshold_squared) {
153           return 1;
154         }
155       }
156       else if (len_squared_v3v3(a, quad[3]->co) < threshold_squared) {
157         copy_v3_v3(a, quad[2]->co);
158         a[axis] = -a[axis];
159         if (len_squared_v3v3(a, quad[1]->co) < threshold_squared) {
160           return 1;
161         }
162       }
163     }
164   }
165
166   return 0;
167 }
168
169 /* Returns true if the quad crosses the plane of symmetry, false otherwise */
170 static bool quad_crosses_symmetry_plane(BMVert *quad[4], const SkinModifierData *smd)
171 {
172   int axis;
173
174   for (axis = 0; axis < 3; axis++) {
175     if (smd->symmetry_axes & (1 << axis)) {
176       bool left = false, right = false;
177       int i;
178
179       for (i = 0; i < 4; i++) {
180         if (quad[i]->co[axis] < 0.0f) {
181           left = true;
182         }
183         else if (quad[i]->co[axis] > 0.0f) {
184           right = true;
185         }
186
187         if (left && right) {
188           return true;
189         }
190       }
191     }
192   }
193
194   return false;
195 }
196
197 /* Returns true if the frame is filled by precisely two faces (and
198  * outputs those faces to fill_faces), otherwise returns false. */
199 static bool skin_frame_find_contained_faces(const Frame *frame, BMFace *fill_faces[2])
200 {
201   BMEdge *diag;
202
203   /* See if the frame is bisected by a diagonal edge */
204   diag = BM_edge_exists(frame->verts[0], frame->verts[2]);
205   if (!diag) {
206     diag = BM_edge_exists(frame->verts[1], frame->verts[3]);
207   }
208
209   if (diag) {
210     return BM_edge_face_pair(diag, &fill_faces[0], &fill_faces[1]);
211   }
212   else {
213     return false;
214   }
215 }
216
217 /* Returns true if hull is successfully built, false otherwise */
218 static bool build_hull(SkinOutput *so, Frame **frames, int totframe)
219 {
220 #ifdef WITH_BULLET
221   BMesh *bm = so->bm;
222   BMOperator op;
223   BMIter iter;
224   BMOIter oiter;
225   BMVert *v;
226   BMFace *f;
227   BMEdge *e;
228   int i, j;
229
230   BM_mesh_elem_hflag_disable_all(bm, BM_VERT, BM_ELEM_TAG, false);
231
232   for (i = 0; i < totframe; i++) {
233     for (j = 0; j < 4; j++) {
234       BM_elem_flag_enable(frames[i]->verts[j], BM_ELEM_TAG);
235     }
236   }
237
238   /* Deselect all faces so that only new hull output faces are
239    * selected after the operator is run */
240   BM_mesh_elem_hflag_disable_all(bm, BM_ALL_NOLOOP, BM_ELEM_SELECT, false);
241
242   BMO_op_initf(
243       bm, &op, (BMO_FLAG_DEFAULTS & ~BMO_FLAG_RESPECT_HIDE), "convex_hull input=%hv", BM_ELEM_TAG);
244   BMO_op_exec(bm, &op);
245
246   if (BMO_error_occurred(bm)) {
247     BMO_op_finish(bm, &op);
248     return false;
249   }
250
251   /* Apply face attributes to hull output */
252   BMO_ITER (f, &oiter, op.slots_out, "geom.out", BM_FACE) {
253     BM_face_normal_update(f);
254     if (so->smd->flag & MOD_SKIN_SMOOTH_SHADING) {
255       BM_elem_flag_enable(f, BM_ELEM_SMOOTH);
256     }
257     f->mat_nr = so->mat_nr;
258   }
259
260   /* Mark interior frames */
261   BMO_ITER (v, &oiter, op.slots_out, "geom_interior.out", BM_VERT) {
262     for (i = 0; i < totframe; i++) {
263       Frame *frame = frames[i];
264
265       if (!frame->detached) {
266         for (j = 0; j < 4; j++) {
267           if (frame->verts[j] == v) {
268             frame->inside_hull[j] = true;
269             frame->detached = true;
270             break;
271           }
272         }
273       }
274     }
275   }
276
277   /* Also mark frames as interior if an edge is not in the hull */
278   for (i = 0; i < totframe; i++) {
279     Frame *frame = frames[i];
280
281     if (!frame->detached && (!BM_edge_exists(frame->verts[0], frame->verts[1]) ||
282                              !BM_edge_exists(frame->verts[1], frame->verts[2]) ||
283                              !BM_edge_exists(frame->verts[2], frame->verts[3]) ||
284                              !BM_edge_exists(frame->verts[3], frame->verts[0]))) {
285       frame->detached = true;
286     }
287   }
288
289   /* Remove triangles that would fill the original frames -- skip if
290    * frame is partially detached */
291   BM_mesh_elem_hflag_disable_all(bm, BM_ALL_NOLOOP, BM_ELEM_TAG, false);
292   for (i = 0; i < totframe; i++) {
293     Frame *frame = frames[i];
294     if (!frame->detached) {
295       BMFace *fill_faces[2];
296
297       /* Check if the frame is filled by precisely two
298        * triangles. If so, delete the triangles and their shared
299        * edge. Otherwise, give up and mark the frame as
300        * detached. */
301       if (skin_frame_find_contained_faces(frame, fill_faces)) {
302         BM_elem_flag_enable(fill_faces[0], BM_ELEM_TAG);
303         BM_elem_flag_enable(fill_faces[1], BM_ELEM_TAG);
304       }
305       else {
306         frame->detached = true;
307       }
308     }
309   }
310
311   /* Check if removing triangles above will create wire triangles,
312    * mark them too */
313   BMO_ITER (e, &oiter, op.slots_out, "geom.out", BM_EDGE) {
314     bool is_wire = true;
315     BM_ITER_ELEM (f, &iter, e, BM_FACES_OF_EDGE) {
316       if (!BM_elem_flag_test(f, BM_ELEM_TAG)) {
317         is_wire = false;
318         break;
319       }
320     }
321     if (is_wire) {
322       BM_elem_flag_enable(e, BM_ELEM_TAG);
323     }
324   }
325
326   BMO_op_finish(bm, &op);
327
328   BM_mesh_delete_hflag_tagged(bm, BM_ELEM_TAG, BM_EDGE | BM_FACE);
329
330   return true;
331 #else
332   UNUSED_VARS(so, frames, totframe, skin_frame_find_contained_faces);
333   return false;
334 #endif
335 }
336
337 /* Returns the average frame side length (frames are rectangular, so
338  * just the average of two adjacent edge lengths) */
339 static float frame_len(const Frame *frame)
340 {
341   return (len_v3v3(frame->co[0], frame->co[1]) + len_v3v3(frame->co[1], frame->co[2])) * 0.5f;
342 }
343
344 static void merge_frame_corners(Frame **frames, int totframe)
345 {
346   float dist, side_a, side_b, thresh, mid[3];
347   int i, j, k, l;
348
349   for (i = 0; i < totframe; i++) {
350     side_a = frame_len(frames[i]);
351
352     /* For each corner of each frame... */
353     for (j = 0; j < 4; j++) {
354
355       /* Ensure the merge target is not itself a merge target */
356       if (frames[i]->merge[j].frame) {
357         continue;
358       }
359
360       for (k = i + 1; k < totframe; k++) {
361         BLI_assert(frames[i] != frames[k]);
362
363         side_b = frame_len(frames[k]);
364         thresh = min_ff(side_a, side_b) / 2.0f;
365
366         /* Compare with each corner of all other frames... */
367         for (l = 0; l < 4; l++) {
368           if (frames[k]->merge[l].frame || frames[k]->merge[l].is_target) {
369             continue;
370           }
371
372           /* Some additional concerns that could be checked
373            * further:
374            *
375            * + Vertex coords are being used for the
376            *   edge-length test, but are also being
377            *   modified, might cause symmetry problems.
378            *
379            * + A frame could be merged diagonally across
380            *   another, would generate a weird (bad) T
381            *   junction
382            */
383
384           /* Check if corners are near each other, where
385            * 'near' is based in the frames' minimum side
386            * length */
387           dist = len_v3v3(frames[i]->co[j], frames[k]->co[l]);
388           if (dist < thresh) {
389             mid_v3_v3v3(mid, frames[i]->co[j], frames[k]->co[l]);
390
391             copy_v3_v3(frames[i]->co[j], mid);
392             copy_v3_v3(frames[k]->co[l], mid);
393
394             frames[k]->merge[l].frame = frames[i];
395             frames[k]->merge[l].corner = j;
396             frames[i]->merge[j].is_target = true;
397
398             /* Can't merge another corner into the same
399              * frame corner, so move on to frame k+1 */
400             break;
401           }
402         }
403       }
404     }
405   }
406 }
407
408 static Frame **collect_hull_frames(
409     int v, SkinNode *frames, const MeshElemMap *emap, const MEdge *medge, int *tothullframe)
410 {
411   SkinNode *f;
412   Frame **hull_frames;
413   int nbr, i;
414
415   (*tothullframe) = emap[v].count;
416   hull_frames = MEM_calloc_arrayN(
417       (*tothullframe), sizeof(Frame *), "hull_from_frames.hull_frames");
418   i = 0;
419   for (nbr = 0; nbr < emap[v].count; nbr++) {
420     const MEdge *e = &medge[emap[v].indices[nbr]];
421     f = &frames[BKE_mesh_edge_other_vert(e, v)];
422     /* Can't have adjacent branch nodes yet */
423     if (f->totframe) {
424       hull_frames[i++] = &f->frames[0];
425     }
426     else {
427       (*tothullframe)--;
428     }
429   }
430
431   return hull_frames;
432 }
433
434 /**************************** Create Frames ***************************/
435
436 static void node_frames_init(SkinNode *nf, int totframe)
437 {
438   int i;
439
440   nf->totframe = totframe;
441   memset(nf->frames, 0, sizeof(nf->frames));
442
443   nf->flag = 0;
444   for (i = 0; i < 2; i++) {
445     nf->seam_edges[i] = -1;
446   }
447 }
448
449 static void create_frame(
450     Frame *frame, const float co[3], const float radius[2], float mat[3][3], float offset)
451 {
452   float rx[3], ry[3], rz[3];
453   int i;
454
455   mul_v3_v3fl(ry, mat[1], radius[0]);
456   mul_v3_v3fl(rz, mat[2], radius[1]);
457
458   add_v3_v3v3(frame->co[3], co, ry);
459   add_v3_v3v3(frame->co[3], frame->co[3], rz);
460
461   sub_v3_v3v3(frame->co[2], co, ry);
462   add_v3_v3v3(frame->co[2], frame->co[2], rz);
463
464   sub_v3_v3v3(frame->co[1], co, ry);
465   sub_v3_v3v3(frame->co[1], frame->co[1], rz);
466
467   add_v3_v3v3(frame->co[0], co, ry);
468   sub_v3_v3v3(frame->co[0], frame->co[0], rz);
469
470   mul_v3_v3fl(rx, mat[0], offset);
471   for (i = 0; i < 4; i++) {
472     add_v3_v3v3(frame->co[i], frame->co[i], rx);
473   }
474 }
475
476 static float half_v2(const float v[2])
477 {
478   return (v[0] + v[1]) * 0.5f;
479 }
480
481 static void end_node_frames(int v,
482                             SkinNode *skin_nodes,
483                             const MVert *mvert,
484                             const MVertSkin *nodes,
485                             const MeshElemMap *emap,
486                             EMat *emat)
487 {
488   const float *rad = nodes[v].radius;
489   float mat[3][3];
490
491   if (emap[v].count == 0) {
492     float avg = half_v2(rad);
493
494     /* For solitary nodes, just build a box (two frames) */
495     node_frames_init(&skin_nodes[v], 2);
496     skin_nodes[v].flag |= (CAP_START | CAP_END);
497
498     /* Hardcoded basis */
499     zero_m3(mat);
500     mat[0][2] = mat[1][0] = mat[2][1] = 1;
501
502     /* Caps */
503     create_frame(&skin_nodes[v].frames[0], mvert[v].co, rad, mat, avg);
504     create_frame(&skin_nodes[v].frames[1], mvert[v].co, rad, mat, -avg);
505   }
506   else {
507     /* For nodes with an incoming edge, create a single (capped) frame */
508     node_frames_init(&skin_nodes[v], 1);
509     skin_nodes[v].flag |= CAP_START;
510
511     /* Use incoming edge for orientation */
512     copy_m3_m3(mat, emat[emap[v].indices[0]].mat);
513     if (emat[emap[v].indices[0]].origin != v) {
514       negate_v3(mat[0]);
515     }
516
517     Frame *frame = &skin_nodes[v].frames[0];
518
519     /* End frame */
520     create_frame(frame, mvert[v].co, rad, mat, 0);
521
522     /* The caps might need to have their normals inverted. So check if they
523      * need to be flipped when creating faces. */
524     float normal[3];
525     normal_quad_v3(normal, frame->co[0], frame->co[1], frame->co[2], frame->co[3]);
526     if (dot_v3v3(mat[0], normal) < 0.0f) {
527       skin_nodes[v].flag |= FLIP_NORMAL;
528     }
529   }
530 }
531
532 /* Returns 1 for seam, 0 otherwise */
533 static int connection_node_mat(float mat[3][3], int v, const MeshElemMap *emap, EMat *emat)
534 {
535   float axis[3], angle, ine[3][3], oute[3][3];
536   EMat *e1, *e2;
537
538   e1 = &emat[emap[v].indices[0]];
539   e2 = &emat[emap[v].indices[1]];
540
541   if (e1->origin != v && e2->origin == v) {
542     copy_m3_m3(ine, e1->mat);
543     copy_m3_m3(oute, e2->mat);
544   }
545   else if (e1->origin == v && e2->origin != v) {
546     copy_m3_m3(ine, e2->mat);
547     copy_m3_m3(oute, e1->mat);
548   }
549   else {
550     return 1;
551   }
552
553   /* Get axis and angle to rotate frame by */
554   angle = angle_normalized_v3v3(ine[0], oute[0]) / 2.0f;
555   cross_v3_v3v3(axis, ine[0], oute[0]);
556   normalize_v3(axis);
557
558   /* Build frame matrix (don't care about X axis here) */
559   copy_v3_v3(mat[0], ine[0]);
560   rotate_normalized_v3_v3v3fl(mat[1], ine[1], axis, angle);
561   rotate_normalized_v3_v3v3fl(mat[2], ine[2], axis, angle);
562
563   return 0;
564 }
565
566 static void connection_node_frames(int v,
567                                    SkinNode *skin_nodes,
568                                    const MVert *mvert,
569                                    const MVertSkin *nodes,
570                                    const MeshElemMap *emap,
571                                    EMat *emat)
572 {
573   const float *rad = nodes[v].radius;
574   float mat[3][3];
575   EMat *e1, *e2;
576
577   if (connection_node_mat(mat, v, emap, emat)) {
578     float avg = half_v2(rad);
579
580     /* Get edges */
581     e1 = &emat[emap[v].indices[0]];
582     e2 = &emat[emap[v].indices[1]];
583
584     /* Handle seam separately to avoid twisting */
585     /* Create two frames, will be hulled to neighbors later */
586     node_frames_init(&skin_nodes[v], 2);
587     skin_nodes[v].flag |= SEAM_FRAME;
588
589     copy_m3_m3(mat, e1->mat);
590     if (e1->origin != v) {
591       negate_v3(mat[0]);
592     }
593     create_frame(&skin_nodes[v].frames[0], mvert[v].co, rad, mat, avg);
594     skin_nodes[v].seam_edges[0] = emap[v].indices[0];
595
596     copy_m3_m3(mat, e2->mat);
597     if (e2->origin != v) {
598       negate_v3(mat[0]);
599     }
600     create_frame(&skin_nodes[v].frames[1], mvert[v].co, rad, mat, avg);
601     skin_nodes[v].seam_edges[1] = emap[v].indices[1];
602
603     return;
604   }
605
606   /* Build regular frame */
607   node_frames_init(&skin_nodes[v], 1);
608   create_frame(&skin_nodes[v].frames[0], mvert[v].co, rad, mat, 0);
609 }
610
611 static SkinNode *build_frames(
612     const MVert *mvert, int totvert, const MVertSkin *nodes, const MeshElemMap *emap, EMat *emat)
613 {
614   SkinNode *skin_nodes;
615   int v;
616
617   skin_nodes = MEM_calloc_arrayN(totvert, sizeof(SkinNode), "build_frames.skin_nodes");
618
619   for (v = 0; v < totvert; v++) {
620     if (emap[v].count <= 1) {
621       end_node_frames(v, skin_nodes, mvert, nodes, emap, emat);
622     }
623     else if (emap[v].count == 2) {
624       connection_node_frames(v, skin_nodes, mvert, nodes, emap, emat);
625     }
626     else {
627       /* Branch node generates no frames */
628     }
629   }
630
631   return skin_nodes;
632 }
633
634 /**************************** Edge Matrices ***************************/
635
636 static void calc_edge_mat(float mat[3][3], const float a[3], const float b[3])
637 {
638   const float z_up[3] = {0, 0, 1};
639   float dot;
640
641   /* X = edge direction */
642   sub_v3_v3v3(mat[0], b, a);
643   normalize_v3(mat[0]);
644
645   dot = dot_v3v3(mat[0], z_up);
646   if (dot > -1 + FLT_EPSILON && dot < 1 - FLT_EPSILON) {
647     /* Y = Z cross x */
648     cross_v3_v3v3(mat[1], z_up, mat[0]);
649     normalize_v3(mat[1]);
650
651     /* Z = x cross y */
652     cross_v3_v3v3(mat[2], mat[0], mat[1]);
653     normalize_v3(mat[2]);
654   }
655   else {
656     mat[1][0] = 1;
657     mat[1][1] = 0;
658     mat[1][2] = 0;
659     mat[2][0] = 0;
660     mat[2][1] = 1;
661     mat[2][2] = 0;
662   }
663 }
664
665 typedef struct {
666   float mat[3][3];
667   int parent_v;
668   int e;
669 } EdgeStackElem;
670
671 static void build_emats_stack(BLI_Stack *stack,
672                               BLI_bitmap *visited_e,
673                               EMat *emat,
674                               const MeshElemMap *emap,
675                               const MEdge *medge,
676                               const MVertSkin *vs,
677                               const MVert *mvert)
678 {
679   EdgeStackElem stack_elem;
680   float axis[3], angle;
681   int i, e, v, parent_v, parent_is_branch;
682
683   BLI_stack_pop(stack, &stack_elem);
684   parent_v = stack_elem.parent_v;
685   e = stack_elem.e;
686
687   /* Skip if edge already visited */
688   if (BLI_BITMAP_TEST(visited_e, e)) {
689     return;
690   }
691
692   /* Mark edge as visited */
693   BLI_BITMAP_ENABLE(visited_e, e);
694
695   /* Process edge */
696
697   parent_is_branch = ((emap[parent_v].count > 2) || (vs[parent_v].flag & MVERT_SKIN_ROOT));
698
699   v = BKE_mesh_edge_other_vert(&medge[e], parent_v);
700   emat[e].origin = parent_v;
701
702   /* If parent is a branch node, start a new edge chain */
703   if (parent_is_branch) {
704     calc_edge_mat(emat[e].mat, mvert[parent_v].co, mvert[v].co);
705   }
706   else {
707     /* Build edge matrix guided by parent matrix */
708     sub_v3_v3v3(emat[e].mat[0], mvert[v].co, mvert[parent_v].co);
709     normalize_v3(emat[e].mat[0]);
710     angle = angle_normalized_v3v3(stack_elem.mat[0], emat[e].mat[0]);
711     cross_v3_v3v3(axis, stack_elem.mat[0], emat[e].mat[0]);
712     normalize_v3(axis);
713     rotate_normalized_v3_v3v3fl(emat[e].mat[1], stack_elem.mat[1], axis, angle);
714     rotate_normalized_v3_v3v3fl(emat[e].mat[2], stack_elem.mat[2], axis, angle);
715   }
716
717   /* Add neighbors to stack */
718   for (i = 0; i < emap[v].count; i++) {
719     /* Add neighbors to stack */
720     copy_m3_m3(stack_elem.mat, emat[e].mat);
721     stack_elem.e = emap[v].indices[i];
722     stack_elem.parent_v = v;
723     BLI_stack_push(stack, &stack_elem);
724   }
725 }
726
727 static EMat *build_edge_mats(const MVertSkin *vs,
728                              const MVert *mvert,
729                              int totvert,
730                              const MEdge *medge,
731                              const MeshElemMap *emap,
732                              int totedge,
733                              bool *has_valid_root)
734 {
735   BLI_Stack *stack;
736   EMat *emat;
737   EdgeStackElem stack_elem;
738   BLI_bitmap *visited_e;
739   int i, v;
740
741   stack = BLI_stack_new(sizeof(stack_elem), "build_edge_mats.stack");
742
743   visited_e = BLI_BITMAP_NEW(totedge, "build_edge_mats.visited_e");
744   emat = MEM_calloc_arrayN(totedge, sizeof(EMat), "build_edge_mats.emat");
745
746   /* Edge matrices are built from the root nodes, add all roots with
747    * children to the stack */
748   for (v = 0; v < totvert; v++) {
749     if (vs[v].flag & MVERT_SKIN_ROOT) {
750       if (emap[v].count >= 1) {
751         const MEdge *e = &medge[emap[v].indices[0]];
752         calc_edge_mat(stack_elem.mat, mvert[v].co, mvert[BKE_mesh_edge_other_vert(e, v)].co);
753         stack_elem.parent_v = v;
754
755         /* Add adjacent edges to stack */
756         for (i = 0; i < emap[v].count; i++) {
757           stack_elem.e = emap[v].indices[i];
758           BLI_stack_push(stack, &stack_elem);
759         }
760
761         *has_valid_root = true;
762       }
763     }
764   }
765
766   while (!BLI_stack_is_empty(stack)) {
767     build_emats_stack(stack, visited_e, emat, emap, medge, vs, mvert);
768   }
769
770   MEM_freeN(visited_e);
771   BLI_stack_free(stack);
772
773   return emat;
774 }
775
776 /************************** Input Subdivision *************************/
777
778 /* Returns number of edge subdivisions, taking into account the radius
779  * of the endpoints and the edge length. If both endpoints are branch
780  * nodes, at least two intermediate frames are required. (This avoids
781  * having any special cases for dealing with sharing a frame between
782  * two hulls.) */
783 static int calc_edge_subdivisions(const MVert *mvert,
784                                   const MVertSkin *nodes,
785                                   const MEdge *e,
786                                   const int *degree)
787 {
788   /* prevent memory errors [#38003] */
789 #define NUM_SUBDIVISIONS_MAX 128
790
791   const MVertSkin *evs[2] = {&nodes[e->v1], &nodes[e->v2]};
792   float avg_radius;
793   const bool v1_branch = degree[e->v1] > 2;
794   const bool v2_branch = degree[e->v2] > 2;
795   int num_subdivisions;
796
797   /* If either end is a branch node marked 'loose', don't subdivide
798    * the edge (or subdivide just twice if both are branches) */
799   if ((v1_branch && (evs[0]->flag & MVERT_SKIN_LOOSE)) ||
800       (v2_branch && (evs[1]->flag & MVERT_SKIN_LOOSE))) {
801     if (v1_branch && v2_branch) {
802       return 2;
803     }
804     else {
805       return 0;
806     }
807   }
808
809   avg_radius = half_v2(evs[0]->radius) + half_v2(evs[1]->radius);
810
811   if (avg_radius != 0.0f) {
812     /* possible (but unlikely) that we overflow INT_MAX */
813     float num_subdivisions_fl;
814     const float edge_len = len_v3v3(mvert[e->v1].co, mvert[e->v2].co);
815     num_subdivisions_fl = (edge_len / avg_radius);
816     if (num_subdivisions_fl < NUM_SUBDIVISIONS_MAX) {
817       num_subdivisions = (int)num_subdivisions_fl;
818     }
819     else {
820       num_subdivisions = NUM_SUBDIVISIONS_MAX;
821     }
822   }
823   else {
824     num_subdivisions = 0;
825   }
826
827   /* If both ends are branch nodes, two intermediate nodes are
828    * required */
829   if (num_subdivisions < 2 && v1_branch && v2_branch) {
830     num_subdivisions = 2;
831   }
832
833   return num_subdivisions;
834
835 #undef NUM_SUBDIVISIONS_MAX
836 }
837
838 /* Take a Mesh and subdivide its edges to keep skin nodes
839  * reasonably close. */
840 static Mesh *subdivide_base(Mesh *orig)
841 {
842   Mesh *result;
843   MVertSkin *orignode, *outnode;
844   MVert *origvert, *outvert;
845   MEdge *origedge, *outedge, *e;
846   MDeformVert *origdvert, *outdvert;
847   int totorigvert, totorigedge;
848   int totsubd, *degree, *edge_subd;
849   int i, j, k, u, v;
850   float radrat;
851
852   orignode = CustomData_get_layer(&orig->vdata, CD_MVERT_SKIN);
853   origvert = orig->mvert;
854   origedge = orig->medge;
855   origdvert = orig->dvert;
856   totorigvert = orig->totvert;
857   totorigedge = orig->totedge;
858
859   /* Get degree of all vertices */
860   degree = MEM_calloc_arrayN(totorigvert, sizeof(int), "degree");
861   for (i = 0; i < totorigedge; i++) {
862     degree[origedge[i].v1]++;
863     degree[origedge[i].v2]++;
864   }
865
866   /* Per edge, store how many subdivisions are needed */
867   edge_subd = MEM_calloc_arrayN((uint)totorigedge, sizeof(int), "edge_subd");
868   for (i = 0, totsubd = 0; i < totorigedge; i++) {
869     edge_subd[i] += calc_edge_subdivisions(origvert, orignode, &origedge[i], degree);
870     BLI_assert(edge_subd[i] >= 0);
871     totsubd += edge_subd[i];
872   }
873
874   MEM_freeN(degree);
875
876   /* Allocate output mesh */
877   result = BKE_mesh_new_nomain_from_template(
878       orig, totorigvert + totsubd, totorigedge + totsubd, 0, 0, 0);
879
880   outvert = result->mvert;
881   outedge = result->medge;
882   outnode = CustomData_get_layer(&result->vdata, CD_MVERT_SKIN);
883   outdvert = result->dvert;
884
885   /* Copy original vertex data */
886   CustomData_copy_data(&orig->vdata, &result->vdata, 0, 0, totorigvert);
887
888   /* Subdivide edges */
889   for (i = 0, v = totorigvert; i < totorigedge; i++) {
890     struct {
891       /* Vertex group number */
892       int def_nr;
893       float w1, w2;
894     } *vgroups = NULL, *vg;
895     int totvgroup = 0;
896
897     e = &origedge[i];
898
899     if (origdvert) {
900       const MDeformVert *dv1 = &origdvert[e->v1];
901       const MDeformVert *dv2 = &origdvert[e->v2];
902       vgroups = MEM_calloc_arrayN(dv1->totweight, sizeof(*vgroups), "vgroup");
903
904       /* Only want vertex groups used by both vertices */
905       for (j = 0; j < dv1->totweight; j++) {
906         vg = NULL;
907         for (k = 0; k < dv2->totweight; k++) {
908           if (dv1->dw[j].def_nr == dv2->dw[k].def_nr) {
909             vg = &vgroups[totvgroup];
910             totvgroup++;
911             break;
912           }
913         }
914
915         if (vg) {
916           vg->def_nr = dv1->dw[j].def_nr;
917           vg->w1 = dv1->dw[j].weight;
918           vg->w2 = dv2->dw[k].weight;
919         }
920       }
921     }
922
923     u = e->v1;
924     radrat = (half_v2(outnode[e->v2].radius) / half_v2(outnode[e->v1].radius));
925     radrat = (radrat + 1) / 2;
926
927     /* Add vertices and edge segments */
928     for (j = 0; j < edge_subd[i]; j++, v++, outedge++) {
929       float r = (j + 1) / (float)(edge_subd[i] + 1);
930       float t = powf(r, radrat);
931
932       /* Interpolate vertex coord */
933       interp_v3_v3v3(outvert[v].co, outvert[e->v1].co, outvert[e->v2].co, t);
934
935       /* Interpolate skin radii */
936       interp_v3_v3v3(outnode[v].radius, orignode[e->v1].radius, orignode[e->v2].radius, t);
937
938       /* Interpolate vertex group weights */
939       for (k = 0; k < totvgroup; k++) {
940         float weight;
941
942         vg = &vgroups[k];
943         weight = interpf(vg->w2, vg->w1, t);
944
945         if (weight > 0) {
946           defvert_add_index_notest(&outdvert[v], vg->def_nr, weight);
947         }
948       }
949
950       outedge->v1 = u;
951       outedge->v2 = v;
952       u = v;
953     }
954
955     if (vgroups) {
956       MEM_freeN(vgroups);
957     }
958
959     /* Link up to final vertex */
960     outedge->v1 = u;
961     outedge->v2 = e->v2;
962     outedge++;
963   }
964
965   MEM_freeN(edge_subd);
966
967   return result;
968 }
969
970 /******************************* Output *******************************/
971
972 /* Can be either quad or triangle */
973 static void add_poly(SkinOutput *so, BMVert *v1, BMVert *v2, BMVert *v3, BMVert *v4)
974 {
975   BMVert *verts[4] = {v1, v2, v3, v4};
976   BMFace *f;
977
978   BLI_assert(v1 != v2 && v1 != v3 && v1 != v4);
979   BLI_assert(v2 != v3 && v2 != v4);
980   BLI_assert(v3 != v4);
981   BLI_assert(v1 && v2 && v3);
982
983   f = BM_face_create_verts(so->bm, verts, v4 ? 4 : 3, NULL, BM_CREATE_NO_DOUBLE, true);
984   BM_face_normal_update(f);
985   if (so->smd->flag & MOD_SKIN_SMOOTH_SHADING) {
986     BM_elem_flag_enable(f, BM_ELEM_SMOOTH);
987   }
988   f->mat_nr = so->mat_nr;
989 }
990
991 static void connect_frames(SkinOutput *so, BMVert *frame1[4], BMVert *frame2[4])
992 {
993   BMVert *q[4][4] = {
994       {frame2[0], frame2[1], frame1[1], frame1[0]},
995       {frame2[1], frame2[2], frame1[2], frame1[1]},
996       {frame2[2], frame2[3], frame1[3], frame1[2]},
997       {frame2[3], frame2[0], frame1[0], frame1[3]},
998   };
999   int i;
1000   bool swap;
1001
1002   /* Check if frame normals need swap */
1003 #if 0
1004   {
1005     /* simple method, works mostly */
1006     float p[3], no[3];
1007     sub_v3_v3v3(p, q[3][0]->co, q[0][0]->co);
1008     normal_quad_v3(no, q[0][0]->co, q[0][1]->co, q[0][2]->co, q[0][3]->co);
1009     swap = dot_v3v3(no, p) > 0;
1010   }
1011 #else
1012   {
1013     /* comprehensive method, accumulate flipping of all faces */
1014     float cent_sides[4][3];
1015     float cent[3];
1016     float dot = 0.0f;
1017
1018     for (i = 0; i < 4; i++) {
1019       mid_v3_v3v3v3v3(cent_sides[i], UNPACK4_EX(, q[i], ->co));
1020     }
1021     mid_v3_v3v3v3v3(cent, UNPACK4(cent_sides));
1022
1023     for (i = 0; i < 4; i++) {
1024       float p[3], no[3];
1025       normal_quad_v3(no, UNPACK4_EX(, q[i], ->co));
1026       sub_v3_v3v3(p, cent, cent_sides[i]);
1027       dot += dot_v3v3(no, p);
1028     }
1029
1030     swap = dot > 0;
1031   }
1032 #endif
1033
1034   for (i = 0; i < 4; i++) {
1035     if (swap) {
1036       add_poly(so, q[i][3], q[i][2], q[i][1], q[i][0]);
1037     }
1038     else {
1039       add_poly(so, q[i][0], q[i][1], q[i][2], q[i][3]);
1040     }
1041   }
1042 }
1043
1044 static void output_frames(BMesh *bm, SkinNode *sn, const MDeformVert *input_dvert)
1045 {
1046   Frame *f;
1047   int i, j;
1048
1049   /* Output all frame verts */
1050   for (i = 0; i < sn->totframe; i++) {
1051     f = &sn->frames[i];
1052     for (j = 0; j < 4; j++) {
1053       if (!f->merge[j].frame) {
1054         BMVert *v = f->verts[j] = BM_vert_create(bm, f->co[j], NULL, BM_CREATE_NOP);
1055
1056         if (input_dvert) {
1057           MDeformVert *dv;
1058           dv = CustomData_bmesh_get(&bm->vdata, v->head.data, CD_MDEFORMVERT);
1059
1060           BLI_assert(dv->totweight == 0);
1061           defvert_copy(dv, input_dvert);
1062         }
1063       }
1064     }
1065   }
1066 }
1067
1068 #define PRINT_HOLE_INFO 0
1069
1070 static void calc_frame_center(float center[3], const Frame *frame)
1071 {
1072   add_v3_v3v3(center, frame->verts[0]->co, frame->verts[1]->co);
1073   add_v3_v3(center, frame->verts[2]->co);
1074   add_v3_v3(center, frame->verts[3]->co);
1075   mul_v3_fl(center, 0.25f);
1076 }
1077
1078 /* Does crappy fan triangulation of poly, may not be so accurate for
1079  * concave faces */
1080 static int isect_ray_poly(const float ray_start[3],
1081                           const float ray_dir[3],
1082                           BMFace *f,
1083                           float *r_lambda)
1084 {
1085   BMVert *v, *v_first = NULL, *v_prev = NULL;
1086   BMIter iter;
1087   float best_dist = FLT_MAX;
1088   bool hit = false;
1089
1090   BM_ITER_ELEM (v, &iter, f, BM_VERTS_OF_FACE) {
1091     if (!v_first) {
1092       v_first = v;
1093     }
1094     else if (v_prev != v_first) {
1095       float dist;
1096       bool curhit;
1097
1098       curhit = isect_ray_tri_v3(ray_start, ray_dir, v_first->co, v_prev->co, v->co, &dist, NULL);
1099       if (curhit && dist < best_dist) {
1100         hit = true;
1101         best_dist = dist;
1102       }
1103     }
1104
1105     v_prev = v;
1106   }
1107
1108   *r_lambda = best_dist;
1109   return hit;
1110 }
1111
1112 /* Reduce the face down to 'n' corners by collapsing the edges;
1113  * returns the new face.
1114  *
1115  * The orig_verts should contain the vertices of 'f'
1116  */
1117 static BMFace *collapse_face_corners(BMesh *bm, BMFace *f, int n, BMVert **orig_verts)
1118 {
1119   int orig_len = f->len;
1120
1121   BLI_assert(n >= 3);
1122   BLI_assert(f->len > n);
1123   if (f->len <= n) {
1124     return f;
1125   }
1126
1127   /* Collapse shortest edge for now */
1128   while (f->len > n) {
1129     BMFace *vf;
1130     BMEdge *shortest_edge;
1131     BMVert *v_safe, *v_merge;
1132     BMOperator op;
1133     BMIter iter;
1134     int i;
1135     BMOpSlot *slot_targetmap;
1136
1137     shortest_edge = BM_face_find_shortest_loop(f)->e;
1138     BMO_op_initf(bm, &op, (BMO_FLAG_DEFAULTS & ~BMO_FLAG_RESPECT_HIDE), "weld_verts");
1139
1140     slot_targetmap = BMO_slot_get(op.slots_in, "targetmap");
1141
1142     /* Note: could probably calculate merges in one go to be
1143      * faster */
1144
1145     v_safe = shortest_edge->v1;
1146     v_merge = shortest_edge->v2;
1147     mid_v3_v3v3(v_safe->co, v_safe->co, v_merge->co);
1148     BMO_slot_map_elem_insert(&op, slot_targetmap, v_merge, v_safe);
1149     BMO_op_exec(bm, &op);
1150     BMO_op_finish(bm, &op);
1151
1152     /* Find the new face */
1153     f = NULL;
1154     BM_ITER_ELEM (vf, &iter, v_safe, BM_FACES_OF_VERT) {
1155       bool wrong_face = false;
1156
1157       for (i = 0; i < orig_len; i++) {
1158         if (orig_verts[i] == v_merge) {
1159           orig_verts[i] = NULL;
1160         }
1161         else if (orig_verts[i] && !BM_vert_in_face(orig_verts[i], vf)) {
1162           wrong_face = true;
1163           break;
1164         }
1165       }
1166
1167       if (!wrong_face) {
1168         f = vf;
1169         break;
1170       }
1171     }
1172
1173     BLI_assert(f);
1174   }
1175
1176   return f;
1177 }
1178
1179 /* Choose a good face to merge the frame with, used in case the frame
1180  * is completely inside the hull. */
1181 static BMFace *skin_hole_target_face(BMesh *bm, Frame *frame)
1182 {
1183   BMFace *f, *isect_target_face, *center_target_face;
1184   BMIter iter;
1185   float frame_center[3];
1186   float frame_normal[3];
1187   float best_isect_dist = FLT_MAX;
1188   float best_center_dist = FLT_MAX;
1189
1190   calc_frame_center(frame_center, frame);
1191   normal_quad_v3(frame_normal,
1192                  frame->verts[3]->co,
1193                  frame->verts[2]->co,
1194                  frame->verts[1]->co,
1195                  frame->verts[0]->co);
1196
1197   /* Use a line intersection test and nearest center test against
1198    * all faces */
1199   isect_target_face = center_target_face = NULL;
1200   BM_ITER_MESH (f, &iter, bm, BM_FACES_OF_MESH) {
1201     float dist, poly_center[3];
1202     int hit;
1203
1204     /* Intersection test */
1205     hit = isect_ray_poly(frame_center, frame_normal, f, &dist);
1206     if (hit && dist < best_isect_dist) {
1207       isect_target_face = f;
1208       best_isect_dist = dist;
1209     }
1210
1211     /* Nearest test */
1212     BM_face_calc_center_median(f, poly_center);
1213     dist = len_v3v3(frame_center, poly_center);
1214     if (dist < best_center_dist) {
1215       center_target_face = f;
1216       best_center_dist = dist;
1217     }
1218   }
1219
1220   f = isect_target_face;
1221   if (!f || best_center_dist < best_isect_dist / 2) {
1222     f = center_target_face;
1223   }
1224
1225   /* This case is unlikely now, but could still happen. Should look
1226    * into splitting edges to make new faces. */
1227 #if PRINT_HOLE_INFO
1228   if (!f) {
1229     printf("no good face found\n");
1230   }
1231 #endif
1232
1233   return f;
1234 }
1235
1236 /* Use edge-length heuristic to choose from eight possible polygon bridges */
1237 static void skin_choose_quad_bridge_order(BMVert *a[4], BMVert *b[4], int best_order[4])
1238 {
1239   int orders[8][4];
1240   float shortest_len;
1241   int i, j;
1242
1243   /* Enumerate all valid orderings */
1244   for (i = 0; i < 4; i++) {
1245     for (j = 0; j < 4; j++) {
1246       orders[i][j] = (j + i) % 4;
1247       orders[i + 4][j] = 3 - ((j + i) % 4);
1248     }
1249   }
1250
1251   shortest_len = FLT_MAX;
1252   for (i = 0; i < 8; i++) {
1253     float len = 0;
1254
1255     /* Get total edge length for this configuration */
1256     for (j = 0; j < 4; j++) {
1257       len += len_squared_v3v3(a[j]->co, b[orders[i][j]]->co);
1258     }
1259
1260     if (len < shortest_len) {
1261       shortest_len = len;
1262       memcpy(best_order, orders[i], sizeof(int) * 4);
1263     }
1264   }
1265 }
1266
1267 static void skin_fix_hole_no_good_verts(BMesh *bm, Frame *frame, BMFace *split_face)
1268 {
1269   BMFace *f;
1270   BMVert *verts[4];
1271   BMVert **vert_buf = NULL;
1272   BLI_array_declare(vert_buf);
1273   BMOIter oiter;
1274   BMOperator op;
1275   int i, best_order[4];
1276   BMOpSlot *slot_targetmap;
1277
1278   BLI_assert(split_face->len >= 3);
1279
1280   /* Extrude the split face */
1281   BM_mesh_elem_hflag_disable_all(bm, BM_FACE, BM_ELEM_TAG, false);
1282   BM_elem_flag_enable(split_face, BM_ELEM_TAG);
1283   BMO_op_initf(bm,
1284                &op,
1285                (BMO_FLAG_DEFAULTS & ~BMO_FLAG_RESPECT_HIDE),
1286                "extrude_discrete_faces faces=%hf",
1287                BM_ELEM_TAG);
1288   BMO_op_exec(bm, &op);
1289
1290   /* Update split face (should only be one new face created
1291    * during extrusion) */
1292   split_face = NULL;
1293   BMO_ITER (f, &oiter, op.slots_out, "faces.out", BM_FACE) {
1294     BLI_assert(!split_face);
1295     split_face = f;
1296   }
1297
1298   BMO_op_finish(bm, &op);
1299
1300   if (split_face->len == 3) {
1301     BMEdge *longest_edge;
1302
1303     /* Need at least four ring edges, so subdivide longest edge if
1304      * face is a triangle */
1305     longest_edge = BM_face_find_longest_loop(split_face)->e;
1306
1307     BM_mesh_elem_hflag_disable_all(bm, BM_EDGE, BM_ELEM_TAG, false);
1308     BM_elem_flag_enable(longest_edge, BM_ELEM_TAG);
1309
1310     BMO_op_callf(bm,
1311                  BMO_FLAG_DEFAULTS,
1312                  "subdivide_edges edges=%he cuts=%i quad_corner_type=%i",
1313                  BM_ELEM_TAG,
1314                  1,
1315                  SUBD_CORNER_STRAIGHT_CUT);
1316   }
1317   else if (split_face->len > 4) {
1318     /* Maintain a dynamic vert array containing the split_face's
1319      * vertices, avoids frequent allocs in collapse_face_corners() */
1320     if (BLI_array_len(vert_buf) < split_face->len) {
1321       BLI_array_grow_items(vert_buf, (split_face->len - BLI_array_len(vert_buf)));
1322     }
1323
1324     /* Get split face's verts */
1325     BM_iter_as_array(bm, BM_VERTS_OF_FACE, split_face, (void **)vert_buf, split_face->len);
1326
1327     /* Earlier edge split operations may have turned some quads
1328      * into higher-degree faces */
1329     split_face = collapse_face_corners(bm, split_face, 4, vert_buf);
1330   }
1331
1332   /* Done with dynamic array, split_face must now be a quad */
1333   BLI_array_free(vert_buf);
1334   BLI_assert(split_face->len == 4);
1335   if (split_face->len != 4) {
1336     return;
1337   }
1338
1339   /* Get split face's verts */
1340   // BM_iter_as_array(bm, BM_VERTS_OF_FACE, split_face, (void **)verts, 4);
1341   BM_face_as_array_vert_quad(split_face, verts);
1342   skin_choose_quad_bridge_order(verts, frame->verts, best_order);
1343
1344   /* Delete split face and merge */
1345   BM_face_kill(bm, split_face);
1346   BMO_op_init(bm, &op, (BMO_FLAG_DEFAULTS & ~BMO_FLAG_RESPECT_HIDE), "weld_verts");
1347   slot_targetmap = BMO_slot_get(op.slots_in, "targetmap");
1348   for (i = 0; i < 4; i++) {
1349     BMO_slot_map_elem_insert(&op, slot_targetmap, verts[i], frame->verts[best_order[i]]);
1350   }
1351   BMO_op_exec(bm, &op);
1352   BMO_op_finish(bm, &op);
1353 }
1354
1355 /* If the frame has some vertices that are inside the hull (detached)
1356  * and some attached, duplicate the attached vertices and take the
1357  * whole frame off the hull. */
1358 static void skin_hole_detach_partially_attached_frame(BMesh *bm, Frame *frame)
1359 {
1360   int i, attached[4], totattached = 0;
1361
1362   /* Get/count attached frame corners */
1363   for (i = 0; i < 4; i++) {
1364     if (!frame->inside_hull[i]) {
1365       attached[totattached++] = i;
1366     }
1367   }
1368
1369   /* Detach everything */
1370   for (i = 0; i < totattached; i++) {
1371     BMVert **av = &frame->verts[attached[i]];
1372     (*av) = BM_vert_create(bm, (*av)->co, *av, BM_CREATE_NOP);
1373   }
1374 }
1375
1376 static void quad_from_tris(BMEdge *e, BMFace *adj[2], BMVert *ndx[4])
1377 {
1378   BMVert *tri[2][3];
1379   BMVert *opp = NULL;
1380   int i, j;
1381
1382   BLI_assert(adj[0]->len == 3 && adj[1]->len == 3);
1383
1384 #if 0
1385   BM_iter_as_array(bm, BM_VERTS_OF_FACE, adj[0], (void **)tri[0], 3);
1386   BM_iter_as_array(bm, BM_VERTS_OF_FACE, adj[1], (void **)tri[1], 3);
1387 #else
1388   BM_face_as_array_vert_tri(adj[0], tri[0]);
1389   BM_face_as_array_vert_tri(adj[1], tri[1]);
1390 #endif
1391
1392   /* Find what the second tri has that the first doesn't */
1393   for (i = 0; i < 3; i++) {
1394     if (tri[1][i] != tri[0][0] && tri[1][i] != tri[0][1] && tri[1][i] != tri[0][2]) {
1395       opp = tri[1][i];
1396       break;
1397     }
1398   }
1399   BLI_assert(opp);
1400
1401   for (i = 0, j = 0; i < 3; i++, j++) {
1402     ndx[j] = tri[0][i];
1403     /* When the triangle edge cuts across our quad-to-be,
1404      * throw in the second triangle's vertex */
1405     if ((tri[0][i] == e->v1 || tri[0][i] == e->v2) &&
1406         (tri[0][(i + 1) % 3] == e->v1 || tri[0][(i + 1) % 3] == e->v2)) {
1407       j++;
1408       ndx[j] = opp;
1409     }
1410   }
1411 }
1412
1413 static void add_quad_from_tris(SkinOutput *so, BMEdge *e, BMFace *adj[2])
1414 {
1415   BMVert *quad[4];
1416
1417   quad_from_tris(e, adj, quad);
1418
1419   add_poly(so, quad[0], quad[1], quad[2], quad[3]);
1420 }
1421
1422 static void hull_merge_triangles(SkinOutput *so, const SkinModifierData *smd)
1423 {
1424   BMIter iter;
1425   BMEdge *e;
1426   HeapSimple *heap;
1427   float score;
1428
1429   heap = BLI_heapsimple_new();
1430
1431   BM_mesh_elem_hflag_disable_all(so->bm, BM_FACE, BM_ELEM_TAG, false);
1432
1433   /* Build heap */
1434   BM_ITER_MESH (e, &iter, so->bm, BM_EDGES_OF_MESH) {
1435     BMFace *adj[2];
1436
1437     /* Only care if the edge is used by exactly two triangles */
1438     if (BM_edge_face_pair(e, &adj[0], &adj[1])) {
1439       if (adj[0]->len == 3 && adj[1]->len == 3) {
1440         BMVert *quad[4];
1441
1442         BLI_assert(BM_face_is_normal_valid(adj[0]));
1443         BLI_assert(BM_face_is_normal_valid(adj[1]));
1444
1445         /* Construct quad using the two triangles adjacent to
1446          * the edge */
1447         quad_from_tris(e, adj, quad);
1448
1449         /* Calculate a score for the quad, higher score for
1450          * triangles being closer to coplanar */
1451         score = ((BM_face_calc_area(adj[0]) + BM_face_calc_area(adj[1])) *
1452                  dot_v3v3(adj[0]->no, adj[1]->no));
1453
1454         /* Check if quad crosses the axis of symmetry */
1455         if (quad_crosses_symmetry_plane(quad, smd)) {
1456           /* Increase score if the triangles form a
1457            * symmetric quad, otherwise don't use it */
1458           if (is_quad_symmetric(quad, smd)) {
1459             score *= 10;
1460           }
1461           else {
1462             continue;
1463           }
1464         }
1465
1466         /* Don't use the quad if it's concave */
1467         if (!is_quad_convex_v3(quad[0]->co, quad[1]->co, quad[2]->co, quad[3]->co)) {
1468           continue;
1469         }
1470
1471         BLI_heapsimple_insert(heap, -score, e);
1472       }
1473     }
1474   }
1475
1476   while (!BLI_heapsimple_is_empty(heap)) {
1477     BMFace *adj[2];
1478
1479     e = BLI_heapsimple_pop_min(heap);
1480
1481     if (BM_edge_face_pair(e, &adj[0], &adj[1])) {
1482       /* If both triangles still free, and if they don't already
1483        * share a border with another face, output as a quad */
1484       if (!BM_elem_flag_test(adj[0], BM_ELEM_TAG) && !BM_elem_flag_test(adj[1], BM_ELEM_TAG) &&
1485           !BM_face_share_face_check(adj[0], adj[1])) {
1486         add_quad_from_tris(so, e, adj);
1487         BM_elem_flag_enable(adj[0], BM_ELEM_TAG);
1488         BM_elem_flag_enable(adj[1], BM_ELEM_TAG);
1489         BM_elem_flag_enable(e, BM_ELEM_TAG);
1490       }
1491     }
1492   }
1493
1494   BLI_heapsimple_free(heap, NULL);
1495
1496   BM_mesh_delete_hflag_tagged(so->bm, BM_ELEM_TAG, BM_EDGE | BM_FACE);
1497 }
1498
1499 static void skin_merge_close_frame_verts(SkinNode *skin_nodes,
1500                                          int totvert,
1501                                          const MeshElemMap *emap,
1502                                          const MEdge *medge)
1503 {
1504   Frame **hull_frames;
1505   int v, tothullframe;
1506
1507   for (v = 0; v < totvert; v++) {
1508     /* Only check branch nodes */
1509     if (!skin_nodes[v].totframe) {
1510       hull_frames = collect_hull_frames(v, skin_nodes, emap, medge, &tothullframe);
1511       merge_frame_corners(hull_frames, tothullframe);
1512       MEM_freeN(hull_frames);
1513     }
1514   }
1515 }
1516
1517 static void skin_update_merged_vertices(SkinNode *skin_nodes, int totvert)
1518 {
1519   int v;
1520
1521   for (v = 0; v < totvert; v++) {
1522     SkinNode *sn = &skin_nodes[v];
1523     int i, j;
1524
1525     for (i = 0; i < sn->totframe; i++) {
1526       Frame *f = &sn->frames[i];
1527
1528       for (j = 0; j < 4; j++) {
1529         if (f->merge[j].frame) {
1530           /* Merge chaining not allowed */
1531           BLI_assert(!f->merge[j].frame->merge[f->merge[j].corner].frame);
1532
1533           f->verts[j] = f->merge[j].frame->verts[f->merge[j].corner];
1534         }
1535       }
1536     }
1537   }
1538 }
1539
1540 static void skin_fix_hull_topology(BMesh *bm, SkinNode *skin_nodes, int totvert)
1541 {
1542   int v;
1543
1544   for (v = 0; v < totvert; v++) {
1545     SkinNode *sn = &skin_nodes[v];
1546     int j;
1547
1548     for (j = 0; j < sn->totframe; j++) {
1549       Frame *f = &sn->frames[j];
1550
1551       if (f->detached) {
1552         BMFace *target_face;
1553
1554         skin_hole_detach_partially_attached_frame(bm, f);
1555
1556         target_face = skin_hole_target_face(bm, f);
1557         if (target_face) {
1558           skin_fix_hole_no_good_verts(bm, f, target_face);
1559         }
1560       }
1561     }
1562   }
1563 }
1564
1565 static void skin_output_end_nodes(SkinOutput *so, SkinNode *skin_nodes, int totvert)
1566 {
1567   int v;
1568
1569   for (v = 0; v < totvert; v++) {
1570     SkinNode *sn = &skin_nodes[v];
1571     /* Assuming here just two frames */
1572     if (sn->flag & SEAM_FRAME) {
1573       BMVert *v_order[4];
1574       int i, order[4];
1575
1576       skin_choose_quad_bridge_order(sn->frames[0].verts, sn->frames[1].verts, order);
1577       for (i = 0; i < 4; i++) {
1578         v_order[i] = sn->frames[1].verts[order[i]];
1579       }
1580       connect_frames(so, sn->frames[0].verts, v_order);
1581     }
1582     else if (sn->totframe == 2) {
1583       connect_frames(so, sn->frames[0].verts, sn->frames[1].verts);
1584     }
1585
1586     if (sn->flag & CAP_START) {
1587       if (sn->flag & FLIP_NORMAL) {
1588         add_poly(so,
1589                  sn->frames[0].verts[0],
1590                  sn->frames[0].verts[1],
1591                  sn->frames[0].verts[2],
1592                  sn->frames[0].verts[3]);
1593       }
1594       else {
1595         add_poly(so,
1596                  sn->frames[0].verts[3],
1597                  sn->frames[0].verts[2],
1598                  sn->frames[0].verts[1],
1599                  sn->frames[0].verts[0]);
1600       }
1601     }
1602     if (sn->flag & CAP_END) {
1603       add_poly(so,
1604                sn->frames[1].verts[0],
1605                sn->frames[1].verts[1],
1606                sn->frames[1].verts[2],
1607                sn->frames[1].verts[3]);
1608     }
1609   }
1610 }
1611
1612 static void skin_output_connections(SkinOutput *so,
1613                                     SkinNode *skin_nodes,
1614                                     const MEdge *medge,
1615                                     int totedge)
1616 {
1617   int e;
1618
1619   for (e = 0; e < totedge; e++) {
1620     SkinNode *a, *b;
1621     a = &skin_nodes[medge[e].v1];
1622     b = &skin_nodes[medge[e].v2];
1623
1624     if (a->totframe && b->totframe) {
1625       if ((a->flag & SEAM_FRAME) || (b->flag & SEAM_FRAME)) {
1626         Frame *fr[2] = {&a->frames[0], &b->frames[0]};
1627         BMVert *v_order[4];
1628         int i, order[4];
1629
1630         if ((a->flag & SEAM_FRAME) && (e != a->seam_edges[0])) {
1631           fr[0]++;
1632         }
1633         if ((b->flag & SEAM_FRAME) && (e != b->seam_edges[0])) {
1634           fr[1]++;
1635         }
1636
1637         skin_choose_quad_bridge_order(fr[0]->verts, fr[1]->verts, order);
1638         for (i = 0; i < 4; i++) {
1639           v_order[i] = fr[1]->verts[order[i]];
1640         }
1641         connect_frames(so, fr[0]->verts, v_order);
1642       }
1643       else {
1644         connect_frames(so, a->frames[0].verts, b->frames[0].verts);
1645       }
1646     }
1647   }
1648 }
1649
1650 static void skin_smooth_hulls(BMesh *bm,
1651                               SkinNode *skin_nodes,
1652                               int totvert,
1653                               const SkinModifierData *smd)
1654 {
1655   BMIter iter, eiter;
1656   BMVert *v;
1657   int i, j, k, skey;
1658
1659   if (smd->branch_smoothing == 0) {
1660     return;
1661   }
1662
1663   /* Mark all frame vertices */
1664   BM_mesh_elem_hflag_disable_all(bm, BM_VERT, BM_ELEM_TAG, false);
1665   for (i = 0; i < totvert; i++) {
1666     for (j = 0; j < skin_nodes[i].totframe; j++) {
1667       Frame *frame = &skin_nodes[i].frames[j];
1668
1669       for (k = 0; k < 4; k++) {
1670         BM_elem_flag_enable(frame->verts[k], BM_ELEM_TAG);
1671       }
1672     }
1673   }
1674
1675   /* Add temporary shapekey layer to store original coordinates */
1676   BM_data_layer_add(bm, &bm->vdata, CD_SHAPEKEY);
1677   skey = CustomData_number_of_layers(&bm->vdata, CD_SHAPEKEY) - 1;
1678   BM_ITER_MESH (v, &iter, bm, BM_VERTS_OF_MESH) {
1679     copy_v3_v3(CustomData_bmesh_get_n(&bm->vdata, v->head.data, CD_SHAPEKEY, skey), v->co);
1680   }
1681
1682   /* Smooth vertices, weight unmarked vertices more strongly (helps
1683    * to smooth frame vertices, but don't want to alter them too
1684    * much) */
1685   BM_ITER_MESH (v, &iter, bm, BM_VERTS_OF_MESH) {
1686     BMEdge *e;
1687     float avg[3];
1688     float weight = smd->branch_smoothing;
1689     int totv = 1;
1690
1691     if (BM_elem_flag_test(v, BM_ELEM_TAG)) {
1692       weight *= 0.5f;
1693     }
1694
1695     copy_v3_v3(avg, v->co);
1696     BM_ITER_ELEM (e, &eiter, v, BM_EDGES_OF_VERT) {
1697       BMVert *other = BM_edge_other_vert(e, v);
1698
1699       add_v3_v3(avg, CustomData_bmesh_get_n(&bm->vdata, other->head.data, CD_SHAPEKEY, skey));
1700       totv++;
1701     }
1702
1703     if (totv > 1) {
1704       mul_v3_fl(avg, 1.0f / (float)totv);
1705       interp_v3_v3v3(v->co, v->co, avg, weight);
1706     }
1707   }
1708
1709   /* Done with original coordinates */
1710   BM_data_layer_free_n(bm, &bm->vdata, CD_SHAPEKEY, skey);
1711 }
1712
1713 /* Returns true if all hulls are successfully built, false otherwise */
1714 static bool skin_output_branch_hulls(
1715     SkinOutput *so, SkinNode *skin_nodes, int totvert, const MeshElemMap *emap, const MEdge *medge)
1716 {
1717   bool result = true;
1718   int v;
1719
1720   for (v = 0; v < totvert; v++) {
1721     SkinNode *sn = &skin_nodes[v];
1722
1723     /* Branch node hulls */
1724     if (!sn->totframe) {
1725       Frame **hull_frames;
1726       int tothullframe;
1727
1728       hull_frames = collect_hull_frames(v, skin_nodes, emap, medge, &tothullframe);
1729       if (!build_hull(so, hull_frames, tothullframe)) {
1730         result = false;
1731       }
1732
1733       MEM_freeN(hull_frames);
1734     }
1735   }
1736
1737   return result;
1738 }
1739
1740 static BMesh *build_skin(SkinNode *skin_nodes,
1741                          int totvert,
1742                          const MeshElemMap *emap,
1743                          const MEdge *medge,
1744                          int totedge,
1745                          const MDeformVert *input_dvert,
1746                          SkinModifierData *smd)
1747 {
1748   SkinOutput so;
1749   int v;
1750
1751   so.smd = smd;
1752   so.bm = BM_mesh_create(&bm_mesh_allocsize_default,
1753                          &((struct BMeshCreateParams){
1754                              .use_toolflags = true,
1755                          }));
1756   so.mat_nr = 0;
1757
1758   /* BMESH_TODO: bumping up the stack level (see MOD_array.c) */
1759   BM_mesh_elem_toolflags_ensure(so.bm);
1760   BMO_push(so.bm, NULL);
1761   bmesh_edit_begin(so.bm, 0);
1762
1763   if (input_dvert) {
1764     BM_data_layer_add(so.bm, &so.bm->vdata, CD_MDEFORMVERT);
1765   }
1766
1767   /* Check for mergeable frame corners around hulls before
1768    * outputting vertices */
1769   skin_merge_close_frame_verts(skin_nodes, totvert, emap, medge);
1770
1771   /* Write out all frame vertices to the mesh */
1772   for (v = 0; v < totvert; v++) {
1773     if (skin_nodes[v].totframe) {
1774       output_frames(so.bm, &skin_nodes[v], input_dvert ? &input_dvert[v] : NULL);
1775     }
1776   }
1777
1778   /* Update vertex pointers for merged frame corners */
1779   skin_update_merged_vertices(skin_nodes, totvert);
1780
1781   if (!skin_output_branch_hulls(&so, skin_nodes, totvert, emap, medge)) {
1782     modifier_setError(&smd->modifier, "Hull error");
1783   }
1784
1785   /* Merge triangles here in the hope of providing better target
1786    * faces for skin_fix_hull_topology() to connect to */
1787   hull_merge_triangles(&so, smd);
1788
1789   /* Using convex hulls may not generate a nice manifold mesh. Two
1790    * problems can occur: an input frame's edges may be inside the
1791    * hull, and/or an input frame's vertices may be inside the hull.
1792    *
1793    * General fix to produce manifold mesh: for any frame that is
1794    * partially detached, first detach it fully, then find a suitable
1795    * existing face to merge with. (Note that we do this after
1796    * creating all hull faces, but before creating any other
1797    * faces.
1798    */
1799   skin_fix_hull_topology(so.bm, skin_nodes, totvert);
1800
1801   skin_smooth_hulls(so.bm, skin_nodes, totvert, smd);
1802
1803   skin_output_end_nodes(&so, skin_nodes, totvert);
1804   skin_output_connections(&so, skin_nodes, medge, totedge);
1805   hull_merge_triangles(&so, smd);
1806
1807   bmesh_edit_end(so.bm, 0);
1808   BMO_pop(so.bm);
1809
1810   return so.bm;
1811 }
1812
1813 static void skin_set_orig_indices(Mesh *mesh)
1814 {
1815   int *orig, totpoly;
1816
1817   totpoly = mesh->totpoly;
1818   orig = CustomData_add_layer(&mesh->pdata, CD_ORIGINDEX, CD_CALLOC, NULL, totpoly);
1819   copy_vn_i(orig, totpoly, ORIGINDEX_NONE);
1820 }
1821
1822 /*
1823  * 0) Subdivide edges (in caller)
1824  * 1) Generate good edge matrices (uses root nodes)
1825  * 2) Generate node frames
1826  * 3) Output vertices and polygons from frames, connections, and hulls
1827  */
1828 static Mesh *base_skin(Mesh *origmesh, SkinModifierData *smd)
1829 {
1830   Mesh *result;
1831   MVertSkin *nodes;
1832   BMesh *bm;
1833   EMat *emat;
1834   SkinNode *skin_nodes;
1835   MeshElemMap *emap;
1836   int *emapmem;
1837   MVert *mvert;
1838   MEdge *medge;
1839   MDeformVert *dvert;
1840   int totvert, totedge;
1841   bool has_valid_root = false;
1842
1843   nodes = CustomData_get_layer(&origmesh->vdata, CD_MVERT_SKIN);
1844
1845   mvert = origmesh->mvert;
1846   dvert = origmesh->dvert;
1847   medge = origmesh->medge;
1848   totvert = origmesh->totvert;
1849   totedge = origmesh->totedge;
1850
1851   BKE_mesh_vert_edge_map_create(&emap, &emapmem, medge, totvert, totedge);
1852
1853   emat = build_edge_mats(nodes, mvert, totvert, medge, emap, totedge, &has_valid_root);
1854   skin_nodes = build_frames(mvert, totvert, nodes, emap, emat);
1855   MEM_freeN(emat);
1856   emat = NULL;
1857
1858   bm = build_skin(skin_nodes, totvert, emap, medge, totedge, dvert, smd);
1859
1860   MEM_freeN(skin_nodes);
1861   MEM_freeN(emap);
1862   MEM_freeN(emapmem);
1863
1864   if (!has_valid_root) {
1865     modifier_setError(
1866         &smd->modifier,
1867         "No valid root vertex found (you need one per mesh island you want to skin)");
1868   }
1869
1870   if (!bm) {
1871     return NULL;
1872   }
1873
1874   result = BKE_mesh_from_bmesh_for_eval_nomain(bm, NULL, origmesh);
1875   BM_mesh_free(bm);
1876
1877   result->runtime.cd_dirty_vert |= CD_MASK_NORMAL;
1878
1879   skin_set_orig_indices(result);
1880
1881   return result;
1882 }
1883
1884 static Mesh *final_skin(SkinModifierData *smd, Mesh *mesh)
1885 {
1886   Mesh *result;
1887
1888   /* Skin node layer is required */
1889   if (!CustomData_get_layer(&mesh->vdata, CD_MVERT_SKIN)) {
1890     return mesh;
1891   }
1892
1893   mesh = subdivide_base(mesh);
1894   result = base_skin(mesh, smd);
1895
1896   BKE_id_free(NULL, mesh);
1897   return result;
1898 }
1899
1900 /**************************** Skin Modifier ***************************/
1901
1902 static void initData(ModifierData *md)
1903 {
1904   SkinModifierData *smd = (SkinModifierData *)md;
1905
1906   /* Enable in editmode by default */
1907   md->mode |= eModifierMode_Editmode;
1908
1909   smd->branch_smoothing = 0;
1910   smd->flag = 0;
1911   smd->symmetry_axes = MOD_SKIN_SYMM_X;
1912 }
1913
1914 static Mesh *applyModifier(ModifierData *md, const ModifierEvalContext *UNUSED(ctx), Mesh *mesh)
1915 {
1916   Mesh *result;
1917
1918   if (!(result = final_skin((SkinModifierData *)md, mesh))) {
1919     return mesh;
1920   }
1921   return result;
1922 }
1923
1924 static void requiredDataMask(Object *UNUSED(ob),
1925                              ModifierData *UNUSED(md),
1926                              CustomData_MeshMasks *r_cddata_masks)
1927 {
1928   r_cddata_masks->vmask |= CD_MASK_MVERT_SKIN | CD_MASK_MDEFORMVERT;
1929 }
1930
1931 ModifierTypeInfo modifierType_Skin = {
1932     /* name */ "Skin",
1933     /* structName */ "SkinModifierData",
1934     /* structSize */ sizeof(SkinModifierData),
1935     /* type */ eModifierTypeType_Constructive,
1936     /* flags */ eModifierTypeFlag_AcceptsMesh | eModifierTypeFlag_SupportsEditmode,
1937
1938     /* copyData */ modifier_copyData_generic,
1939
1940     /* deformVerts */ NULL,
1941     /* deformMatrices */ NULL,
1942     /* deformVertsEM */ NULL,
1943     /* deformMatricesEM */ NULL,
1944     /* applyModifier */ applyModifier,
1945
1946     /* initData */ initData,
1947     /* requiredDataMask */ requiredDataMask,
1948     /* freeData */ NULL,
1949     /* isDisabled */ NULL,
1950     /* updateDepsgraph */ NULL,
1951     /* dependsOnTime */ NULL,
1952     /* dependsOnNormals */ NULL,
1953     /* foreachObjectLink */ NULL,
1954     /* foreachIDLink */ NULL,
1955     /* freeRuntimeData */ NULL,
1956 };