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