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