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