Bevel: change 'adj' pattern to use subdivision.
[blender.git] / source / blender / bmesh / tools / bmesh_bevel.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):
19  *         Joseph Eagar,
20  *         Aleksandr Mokhov,
21  *         Howard Trickey,
22  *         Campbell Barton
23  *
24  * ***** END GPL LICENSE BLOCK *****
25  */
26
27 /** \file blender/bmesh/tools/bmesh_bevel.c
28  *  \ingroup bmesh
29  */
30
31 #include "MEM_guardedalloc.h"
32
33 #include "DNA_object_types.h"
34 #include "DNA_meshdata_types.h"
35
36 #include "BLI_array.h"
37 #include "BLI_alloca.h"
38 #include "BLI_gsqueue.h"
39 #include "BLI_math.h"
40 #include "BLI_memarena.h"
41
42 #include "BKE_customdata.h"
43 #include "BKE_deform.h"
44
45 #include "bmesh.h"
46 #include "bmesh_bevel.h"  /* own include */
47
48 #include "./intern/bmesh_private.h"
49
50 /* use new way of doing ADJ pattern */
51 #define USE_ADJ_SUBDIV
52
53 #define BEVEL_EPSILON_D  1e-6
54 #define BEVEL_EPSILON    1e-6f
55 #define BEVEL_EPSILON_SQ 1e-12f
56
57 /* happens far too often, uncomment for development */
58 // #define BEVEL_ASSERT_PROJECT
59
60 /* for testing */
61 // #pragma GCC diagnostic error "-Wpadded"
62
63 /* Constructed vertex, sometimes later instantiated as BMVert */
64 typedef struct NewVert {
65         BMVert *v;
66         float co[3];
67 //      int _pad;
68 } NewVert;
69
70 struct BoundVert;
71
72 /* Data for one end of an edge involved in a bevel */
73 typedef struct EdgeHalf {
74         struct EdgeHalf *next, *prev;   /* in CCW order */
75         BMEdge *e;                  /* original mesh edge */
76         BMFace *fprev;              /* face between this edge and previous, if any */
77         BMFace *fnext;              /* face between this edge and next, if any */
78         struct BoundVert *leftv;    /* left boundary vert (looking along edge to end) */
79         struct BoundVert *rightv;   /* right boundary vert, if beveled */
80         int   seg;                  /* how many segments for the bevel */
81         float offset_l;             /* offset for this edge, on left side */
82         float offset_r;             /* offset for this edge, on right side */
83         float offset_l_spec;        /* user specification for offset_l */
84         float offset_r_spec;        /* user specification for offset_r */
85         bool is_bev;                /* is this edge beveled? */
86         bool is_rev;                /* is e->v2 the vertex at this end? */
87         bool is_seam;               /* is e a seam for custom loopdata (e.g., UVs)? */
88 //      int _pad;
89 } EdgeHalf;
90
91 /* Profile specification.
92  * For now, only have round profiles so only need midpoint.
93  * The start and end points of the profile are stored separately.
94  * TODO: generalize to superellipse profiles.
95  */
96 typedef struct Profile {
97         float midco[3];      /* mid control point for profile */
98 } Profile;
99
100 /* An element in a cyclic boundary of a Vertex Mesh (VMesh) */
101 typedef struct BoundVert {
102         struct BoundVert *next, *prev;  /* in CCW order */
103         NewVert nv;
104         EdgeHalf *efirst;   /* first of edges attached here: in CCW order */
105         EdgeHalf *elast;
106         EdgeHalf *ebev;     /* beveled edge whose left side is attached here, if any */
107         int index;          /* used for vmesh indexing */
108         Profile profile;    /* edge profile between this and next BoundVert */
109         bool any_seam;      /* are any of the edges attached here seams? */
110 //      int _pad;
111 } BoundVert;
112
113 /* Mesh structure replacing a vertex */
114 typedef struct VMesh {
115         NewVert *mesh;           /* allocated array - size and structure depends on kind */
116         BoundVert *boundstart;   /* start of boundary double-linked list */
117         int count;               /* number of vertices in the boundary */
118         int seg;                 /* common # of segments for segmented edges */
119         enum {
120                 M_NONE,         /* no polygon mesh needed */
121                 M_POLY,         /* a simple polygon */
122                 M_ADJ,          /* "adjacent edges" mesh pattern */
123                 M_ADJ_SUBDIV,   /* like M_ADJ, but using subdivision */
124                 M_TRI_FAN,      /* a simple polygon - fan filled */
125                 M_QUAD_STRIP,   /* a simple polygon - cut into parallel strips */
126         } mesh_kind;
127 //      int _pad;
128 } VMesh;
129
130 /* Data for a vertex involved in a bevel */
131 typedef struct BevVert {
132         BMVert *v;          /* original mesh vertex */
133         int edgecount;          /* total number of edges around the vertex */
134         int selcount;           /* number of selected edges around the vertex */
135         float offset;           /* offset for this vertex, if vertex_only bevel */
136         bool any_seam;                  /* any seams on attached edges? */
137         bool visited;           /* used in graph traversal */
138         EdgeHalf *edges;        /* array of size edgecount; CCW order from vertex normal side */
139         VMesh *vmesh;           /* mesh structure for replacing vertex */
140 } BevVert;
141
142 /* Bevel parameters and state */
143 typedef struct BevelParams {
144         /* hash of BevVert for each vertex involved in bevel
145          * GHash: (key=(BMVert *), value=(BevVert *)) */
146         GHash    *vert_hash;
147         MemArena *mem_arena;    /* use for all allocs while bevel runs, if we need to free we can switch to mempool */
148
149         float offset;           /* blender units to offset each side of a beveled edge */
150         int offset_type;        /* how offset is measured; enum defined in bmesh_operators.h */
151         int seg;                /* number of segments in beveled edge profile */
152         bool vertex_only;       /* bevel vertices only */
153         bool use_weights;       /* bevel amount affected by weights on edges or verts */
154         bool preserve_widths;   /* should bevel prefer widths over angles, if forced to choose? */
155         bool limit_offset;      /* should offsets be limited by collisions? */
156         const struct MDeformVert *dvert; /* vertex group array, maybe set if vertex_only */
157         int vertex_group;       /* vertex group index, maybe set if vertex_only */
158 } BevelParams;
159
160 // #pragma GCC diagnostic ignored "-Wpadded"
161
162 // #include "bevdebug.c"
163
164 /* Make a new BoundVert of the given kind, insert it at the end of the circular linked
165  * list with entry point bv->boundstart, and return it. */
166 static BoundVert *add_new_bound_vert(MemArena *mem_arena, VMesh *vm, const float co[3])
167 {
168         BoundVert *ans = (BoundVert *)BLI_memarena_alloc(mem_arena, sizeof(BoundVert));
169
170         copy_v3_v3(ans->nv.co, co);
171         if (!vm->boundstart) {
172                 ans->index = 0;
173                 vm->boundstart = ans;
174                 ans->next = ans->prev = ans;
175         }
176         else {
177                 BoundVert *tail = vm->boundstart->prev;
178                 ans->index = tail->index + 1;
179                 ans->prev = tail;
180                 ans->next = vm->boundstart;
181                 tail->next = ans;
182                 vm->boundstart->prev = ans;
183         }
184         vm->count++;
185         return ans;
186 }
187
188 BLI_INLINE void adjust_bound_vert(BoundVert *bv, const float co[3])
189 {
190         copy_v3_v3(bv->nv.co, co);
191 }
192
193 /* Mesh verts are indexed (i, j, k) where
194  * i = boundvert index (0 <= i < nv)
195  * j = ring index (0 <= j <= ns2)
196  * k = segment index (0 <= k <= ns)
197  * Not all of these are used, and some will share BMVerts */
198 static NewVert *mesh_vert(VMesh *vm, int i, int j, int k)
199 {
200         int nj = (vm->seg / 2) + 1;
201         int nk = vm->seg + 1;
202
203         return &vm->mesh[i * nk * nj  + j * nk + k];
204 }
205
206 static void create_mesh_bmvert(BMesh *bm, VMesh *vm, int i, int j, int k, BMVert *eg)
207 {
208         NewVert *nv = mesh_vert(vm, i, j, k);
209         nv->v = BM_vert_create(bm, nv->co, eg, BM_CREATE_NOP);
210         BM_elem_flag_disable(nv->v, BM_ELEM_TAG);
211 }
212
213 static void copy_mesh_vert(VMesh *vm, int ito, int jto, int kto,
214                            int ifrom, int jfrom, int kfrom)
215 {
216         NewVert *nvto, *nvfrom;
217
218         nvto = mesh_vert(vm, ito, jto, kto);
219         nvfrom = mesh_vert(vm, ifrom, jfrom, kfrom);
220         nvto->v = nvfrom->v;
221         copy_v3_v3(nvto->co, nvfrom->co);
222 }
223
224 /* find the EdgeHalf in bv's array that has edge bme */
225 static EdgeHalf *find_edge_half(BevVert *bv, BMEdge *bme)
226 {
227         int i;
228
229         for (i = 0; i < bv->edgecount; i++) {
230                 if (bv->edges[i].e == bme)
231                         return &bv->edges[i];
232         }
233         return NULL;
234 }
235
236 /* find the BevVert corresponding to BMVert bmv */
237 static BevVert *find_bevvert(BevelParams *bp, BMVert *bmv)
238 {
239         return BLI_ghash_lookup(bp->vert_hash, bmv);
240 }
241
242 /* Find the EdgeHalf representing the other end of e->e.
243  * Return other end's BevVert in *bvother, if r_bvother is provided.
244  * That may not have been constructed yet, in which case return NULL. */
245 static EdgeHalf *find_other_end_edge_half(BevelParams *bp, EdgeHalf *e, BevVert **r_bvother)
246 {
247         BevVert *bvo;
248         EdgeHalf *eother;
249
250         bvo = find_bevvert(bp, e->is_rev ? e->e->v1 : e->e->v2);
251         if (bvo) {
252                 if (r_bvother)
253                         *r_bvother = bvo;
254                 eother = find_edge_half(bvo, e->e);
255                 BLI_assert(eother != NULL);
256                 return eother;
257         }
258         else if (r_bvother) {
259                 *r_bvother = NULL;
260         }
261         return NULL;
262 }
263
264 static bool other_edge_half_visited(BevelParams *bp, EdgeHalf *e)
265 {
266         BevVert *bvo;
267
268         bvo = find_bevvert(bp, e->is_rev ? e->e->v1 : e->e->v2);
269         if (bvo)
270                 return bvo->visited;
271         else
272                 return false;
273 }
274
275 static bool edge_half_offset_changed(EdgeHalf *e)
276 {
277         return e->offset_l != e->offset_l_spec ||
278                e->offset_r != e->offset_r_spec;
279 }
280
281 static bool any_edge_half_offset_changed(BevVert *bv)
282 {
283         int i;
284
285         for (i = 0; i < bv->edgecount; i++) {
286                 if (edge_half_offset_changed(&bv->edges[i]))
287                         return true;
288         }
289         return false;
290 }
291
292 /* Return the next EdgeHalf after from_e that is beveled.
293  * If from_e is NULL, find the first beveled edge. */
294 static EdgeHalf *next_bev(BevVert *bv, EdgeHalf *from_e)
295 {
296         EdgeHalf *e;
297
298         if (from_e == NULL)
299                 from_e = &bv->edges[bv->edgecount - 1];
300         e = from_e;
301         do {
302                 if (e->is_bev) {
303                         return e;
304                 }
305         } while ((e = e->next) != from_e);
306         return NULL;
307 }
308
309 /* Return a good representative face (for materials, etc.) for faces
310  * created around/near BoundVert v */
311 static BMFace *boundvert_rep_face(BoundVert *v)
312 {
313         BLI_assert(v->efirst != NULL && v->elast != NULL);
314         if (v->efirst->fnext == v->elast->fprev)
315                 return v->efirst->fnext;
316         else if (v->efirst->fnext)
317                 return v->efirst->fnext;
318         else
319                 return v->elast->fprev;
320 }
321
322 /**
323  * Make ngon from verts alone.
324  * Make sure to properly copy face attributes and do custom data interpolation from
325  * corresponding elements of face_arr, if that is non-NULL, else from facerep.
326  *
327  * \note ALL face creation goes through this function, this is important to keep!
328  */
329 static BMFace *bev_create_ngon(BMesh *bm, BMVert **vert_arr, const int totv,
330                                BMFace **face_arr, BMFace *facerep, bool do_interp)
331 {
332         BMIter iter;
333         BMLoop *l;
334         BMFace *f, *interp_f;
335         int i;
336
337         f = BM_face_create_verts(bm, vert_arr, totv, facerep, BM_CREATE_NOP, true);
338
339         if ((facerep || (face_arr && face_arr[0])) && f) {
340                 BM_elem_attrs_copy(bm, bm, facerep ? facerep : face_arr[0], f);
341                 if (do_interp) {
342                         i = 0;
343                         BM_ITER_ELEM (l, &iter, f, BM_LOOPS_OF_FACE) {
344                                 if (face_arr) {
345                                         /* assume loops of created face are in same order as verts */
346                                         BLI_assert(l->v == vert_arr[i]);
347                                         interp_f = face_arr[i];
348                                 }
349                                 else {
350                                         interp_f = facerep;
351                                 }
352                                 if (interp_f)
353                                         BM_loop_interp_from_face(bm, l, interp_f, TRUE, TRUE);
354                                 i++;
355                         }
356                 }
357         }
358
359         /* not essential for bevels own internal logic,
360          * this is done so the operator can select newly created faces */
361         if (f) {
362                 BM_elem_flag_enable(f, BM_ELEM_TAG);
363         }
364
365         return f;
366 }
367
368 static BMFace *bev_create_quad_tri(BMesh *bm, BMVert *v1, BMVert *v2, BMVert *v3, BMVert *v4,
369                                    BMFace *facerep, bool do_interp)
370 {
371         BMVert *varr[4] = {v1, v2, v3, v4};
372         return bev_create_ngon(bm, varr, v4 ? 4 : 3, NULL, facerep, do_interp);
373 }
374
375 static BMFace *bev_create_quad_tri_ex(BMesh *bm, BMVert *v1, BMVert *v2, BMVert *v3, BMVert *v4,
376                                       BMFace *f1, BMFace *f2, BMFace *f3, BMFace *f4)
377 {
378         BMVert *varr[4] = {v1, v2, v3, v4};
379         BMFace *farr[4] = {f1, f2, f3, f4};
380         return bev_create_ngon(bm, varr, v4 ? 4 : 3, farr, f1, true);
381 }
382
383
384 /* Is Loop layer layer_index contiguous across shared vertex of l1 and l2? */
385 static bool contig_ldata_across_loops(BMesh *bm, BMLoop *l1, BMLoop *l2,
386                                       int layer_index)
387 {
388         const int offset = bm->ldata.layers[layer_index].offset;
389         const int type = bm->ldata.layers[layer_index].type;
390
391         return CustomData_data_equals(type,
392                                       (char *)l1->head.data + offset,
393                                       (char *)l2->head.data + offset);
394 }
395
396 /* Are all loop layers with have math (e.g., UVs) contiguous from face f1 to face f2 across edge e? */
397 static bool contig_ldata_across_edge(BMesh *bm, BMEdge *e, BMFace *f1, BMFace *f2)
398 {
399         BMLoop *lef1, *lef2;
400         BMLoop *lv1f1, *lv1f2, *lv2f1, *lv2f2;
401         BMVert *v1, *v2;
402         int i;
403
404         if (bm->ldata.totlayer == 0)
405                 return true;
406
407         v1 = e->v1;
408         v2 = e->v2;
409         if (!BM_edge_loop_pair(e, &lef1, &lef2))
410                 return false;
411         if (lef1->f == f2) {
412                 SWAP(BMLoop *, lef1, lef2);
413         }
414
415         if (lef1->v == v1) {
416                 lv1f1 = lef1;
417                 lv2f1 = BM_face_other_edge_loop(f1, e, v2);
418         }
419         else {
420                 lv2f1 = lef1;
421                 lv1f1 = BM_face_other_edge_loop(f1, e, v1);
422         }
423
424         if (lef2->v == v1) {
425                 lv1f2 = lef2;
426                 lv2f2 = BM_face_other_edge_loop(f2, e, v2);
427         }
428         else {
429                 lv2f2 = lef2;
430                 lv1f2 = BM_face_other_edge_loop(f2, e, v1);
431         }
432
433         for (i = 0; i < bm->ldata.totlayer; i++) {
434                 if (CustomData_layer_has_math(&bm->ldata, i) &&
435                     (!contig_ldata_across_loops(bm, lv1f1, lv1f2, i) ||
436                      !contig_ldata_across_loops(bm, lv2f1, lv2f2, i)))
437                 {
438                         return false;
439                 }
440         }
441         return true;
442 }
443
444 /* Like bev_create_quad_tri, but when verts straddle an old edge.
445  *        e
446  *        |
447  *  v1+---|---+v4
448  *    |   |   |
449  *    |   |   |
450  *  v2+---|---+v3
451  *        |
452  *    f1  |  f2
453  *
454  * Most CustomData for loops can be interpolated in their respective
455  * faces' loops, but for UVs and other 'has_math_cd' layers, only
456  * do this if the UVs are continuous across the edge e, otherwise pick
457  * one side (f1, arbitrarily), and interpolate them all on that side.
458  * For face data, use f1 (arbitrarily) as face representative. */
459 static BMFace *bev_create_quad_straddle(BMesh *bm, BMVert *v1, BMVert *v2, BMVert *v3, BMVert *v4,
460         BMFace *f1, BMFace *f2, bool is_seam)
461 {
462         BMFace *f, *facerep;
463         BMLoop *l;
464         BMIter iter;
465
466         f = bev_create_quad_tri(bm, v1, v2, v3, v4, f1, false);
467
468         if (!f)
469                 return NULL;
470
471         BM_ITER_ELEM (l, &iter, f, BM_LOOPS_OF_FACE) {
472                 if (is_seam || l->v == v1 || l->v == v2)
473                         facerep = f1;
474                 else
475                         facerep = f2;
476                 if (facerep)
477                         BM_loop_interp_from_face(bm, l, facerep, TRUE, TRUE);
478         }
479         return f;
480 }
481
482 /* Merge (using average) all the UV values for loops of v's faces.
483  * Caller should ensure that no seams are violated by doing this. */
484 static void bev_merge_uvs(BMesh *bm, BMVert *v)
485 {
486         BMIter iter;
487         MLoopUV *luv;
488         BMLoop *l;
489         float uv[2];
490         int n;
491         int cd_loop_uv_offset = CustomData_get_offset(&bm->ldata, CD_MLOOPUV);
492
493         if (cd_loop_uv_offset == -1)
494                 return;
495
496         n = 0;
497         zero_v2(uv);
498         BM_ITER_ELEM (l, &iter, v, BM_LOOPS_OF_VERT) {
499                 luv = BM_ELEM_CD_GET_VOID_P(l, cd_loop_uv_offset);
500                 add_v2_v2(uv, luv->uv);
501                 n++;
502         }
503         if (n > 1) {
504                 mul_v2_fl(uv, 1.0f / (float)n);
505                 BM_ITER_ELEM (l, &iter, v, BM_LOOPS_OF_VERT) {
506                         luv = BM_ELEM_CD_GET_VOID_P(l, cd_loop_uv_offset);
507                         copy_v2_v2(luv->uv, uv);
508                 }
509         }
510 }
511
512 /* Calculate coordinates of a point a distance d from v on e->e and return it in slideco */
513 static void slide_dist(EdgeHalf *e, BMVert *v, float d, float slideco[3])
514 {
515         float dir[3], len;
516
517         sub_v3_v3v3(dir, v->co, BM_edge_other_vert(e->e, v)->co);
518         len = normalize_v3(dir);
519         if (d > len)
520                 d = len - (float)(50.0 * BEVEL_EPSILON_D);
521         copy_v3_v3(slideco, v->co);
522         madd_v3_v3fl(slideco, dir, -d);
523 }
524
525 /*
526  * Calculate the meeting point between the offset edges for e1 and e2, putting answer in meetco.
527  * e1 and e2 share vertex v and face f (may be NULL) and viewed from the normal side of
528  * the bevel vertex,  e1 precedes e2 in CCW order.
529  * Offset edge is on right of both edges, where e1 enters v and e2 leave it.
530  * When offsets are equal, the new point is on the edge bisector, with length offset/sin(angle/2),
531  * but if the offsets are not equal (allowing for this, as bevel modifier has edge weights that may
532  * lead to different offsets) then meeting point can be found be intersecting offset lines.
533  * If making the meeting point significantly changes the left or right offset from the user spec,
534  * record the change in offset_l (or offset_r); later we can tell that a change has happened because
535  * the offset will differ from its original value in offset_l_spec (or offset_r_spec).
536  */
537 static void offset_meet(EdgeHalf *e1, EdgeHalf *e2, BMVert *v, BMFace *f, float meetco[3])
538 {
539         float dir1[3], dir2[3], norm_v[3], norm_perp1[3], norm_perp2[3],
540               off1a[3], off1b[3], off2a[3], off2b[3], isect2[3], ang, d;
541
542         /* get direction vectors for two offset lines */
543         sub_v3_v3v3(dir1, v->co, BM_edge_other_vert(e1->e, v)->co);
544         sub_v3_v3v3(dir2, BM_edge_other_vert(e2->e, v)->co, v->co);
545
546         ang = angle_v3v3(dir1, dir2);
547         if (ang < 100.0f * BEVEL_EPSILON) {
548                 /* special case: e1 and e2 are parallel; put offset point perp to both, from v.
549                  * need to find a suitable plane.
550                  * if offsets are different, we're out of luck: just use e1->offset_r */
551                 if (f)
552                         copy_v3_v3(norm_v, f->no);
553                 else
554                         copy_v3_v3(norm_v, v->no);
555                 cross_v3_v3v3(norm_perp1, dir1, norm_v);
556                 normalize_v3(norm_perp1);
557                 copy_v3_v3(off1a, v->co);
558                 madd_v3_v3fl(off1a, norm_perp1, e1->offset_r);
559                 if (e2->offset_l != e1->offset_r)
560                         e2->offset_l = e1->offset_r;
561                 copy_v3_v3(meetco, off1a);
562         }
563         else if (fabsf(ang - (float)M_PI) < 100.0f * BEVEL_EPSILON) {
564                 /* special case e1 and e2 are antiparallel, so bevel is into
565                  * a zero-area face.  Just make the offset point on the
566                  * common line, at offset distance from v. */
567                 slide_dist(e2, v, e1->offset_r, meetco);
568                 if (e2->offset_l != e1->offset_r)
569                         e2->offset_l = e1->offset_r;
570         }
571         else {
572                 /* Get normal to plane where meet point should be,
573                  * using cross product instead of f->no in case f is non-planar.
574                  * If e1-v-e2 is a reflex angle (viewed from vertex normal side), need to flip*/
575                 cross_v3_v3v3(norm_v, dir2, dir1);
576                 normalize_v3(norm_v);
577                 if (dot_v3v3(norm_v, v->no) < 0.0f)
578                         negate_v3(norm_v);
579
580                 /* get vectors perp to each edge, perp to norm_v, and pointing into face */
581                 cross_v3_v3v3(norm_perp1, dir1, norm_v);
582                 cross_v3_v3v3(norm_perp2, dir2, norm_v);
583                 normalize_v3(norm_perp1);
584                 normalize_v3(norm_perp2);
585
586                 /* get points that are offset distances from each line, then another point on each line */
587                 copy_v3_v3(off1a, v->co);
588                 madd_v3_v3fl(off1a, norm_perp1, e1->offset_r);
589                 add_v3_v3v3(off1b, off1a, dir1);
590                 copy_v3_v3(off2a, v->co);
591                 madd_v3_v3fl(off2a, norm_perp2, e2->offset_l);
592                 add_v3_v3v3(off2b, off2a, dir2);
593
594                 /* intersect the lines; by construction they should be on the same plane and not parallel */
595                 if (!isect_line_line_v3(off1a, off1b, off2a, off2b, meetco, isect2)) {
596 #ifdef BEVEL_ASSERT_PROJECT
597                         BLI_assert(!"offset_meet failure");
598 #endif
599                         copy_v3_v3(meetco, off1a);  /* just to do something */
600                         d = dist_to_line_v3(meetco, v->co, BM_edge_other_vert(e2->e, v)->co);
601                         if (fabsf(d - e2->offset_l) > BEVEL_EPSILON)
602                                 e2->offset_l = d;
603                 }
604         }
605 }
606
607 /* Calculate the meeting point between e1 and e2 (one of which should have zero offsets),
608  * where e1 precedes e2 in CCW order around their common vertex v (viewed from normal side).
609  * If r_angle is provided, return the angle between e and emeet in *r_angle.
610  * If the angle is 0, or it is 180 degrees or larger, there will be no meeting point;
611  * return false in that case, else true */
612 static bool offset_meet_edge(EdgeHalf *e1, EdgeHalf *e2, BMVert *v,  float meetco[3], float *r_angle)
613 {
614         float dir1[3], dir2[3], fno[3], ang, sinang;
615
616         sub_v3_v3v3(dir1, BM_edge_other_vert(e1->e, v)->co, v->co);
617         sub_v3_v3v3(dir2, BM_edge_other_vert(e2->e, v)->co, v->co);
618         normalize_v3(dir1);
619         normalize_v3(dir2);
620
621         /* find angle from dir1 to dir2 as viewed from vertex normal side */
622         ang = angle_normalized_v3v3(dir1, dir2);
623         if (ang < BEVEL_EPSILON) {
624                 if (r_angle)
625                         *r_angle = 0.0f;
626                 return false;
627         }
628         cross_v3_v3v3(fno, dir1, dir2);
629         if (dot_v3v3(fno, v->no) < 0.0f)
630                 ang = 2.0f * (float)M_PI - ang;  /* angle is reflex */
631         if (r_angle)
632                 *r_angle = ang;
633
634         if (ang - (float)M_PI > BEVEL_EPSILON)
635                 return false;
636
637         sinang = sinf(ang);
638         copy_v3_v3(meetco, v->co);
639         if (e1->offset_r == 0.0f)
640                 madd_v3_v3fl(meetco, dir1, e2->offset_l / sinang);
641         else
642                 madd_v3_v3fl(meetco, dir2, e1->offset_r / sinang);
643         return true;
644 }
645
646 /* Calculate the best place for a meeting point for the offsets from edges e1 and e2
647  * on the in-between edge emid.  Viewed from the vertex normal side, the CCW
648  * order of these edges is e1, emid, e2.
649  * The offsets probably do not meet at a common point on emid, so need to pick
650  * one that causes the least problems. If the other end of one of e1 or e2 has been visited
651  * already, prefer to keep the offset the same on this end.
652  * Otherwise, pick a point between the two intersection points on emid that minimizes
653  * the sum of squares of errors from desired offset. */
654 static void offset_on_edge_between(BevelParams *bp, EdgeHalf *e1, EdgeHalf *e2, EdgeHalf *emid,
655                                    BMVert *v, float meetco[3])
656 {
657         float d, ang1, ang2, sina1, sina2, lambda;
658         float meet1[3], meet2[3];
659         bool visited1, visited2, ok1, ok2;
660
661         BLI_assert(e1->is_bev && e2->is_bev && !emid->is_bev);
662
663         visited1 = other_edge_half_visited(bp, e1);
664         visited2 = other_edge_half_visited(bp, e2);
665
666         ok1 = offset_meet_edge(e1, emid, v, meet1, &ang1);
667         ok2 = offset_meet_edge(emid, e2, v, meet2, &ang2);
668         if (ok1 && ok2) {
669                 if (visited1 && !visited2) {
670                         copy_v3_v3(meetco, meet1);
671                 }
672                 else if (!visited1 && visited2) {
673                         copy_v3_v3(meetco, meet2);
674                 }
675                 else {
676                         /* find best compromise meet point */
677                         sina1 = sinf(ang1);
678                         sina2 = sinf(ang2);
679                         lambda = sina2 * sina2 / (sina1 * sina1 + sina2 * sina2);
680                         interp_v3_v3v3(meetco, meet1, meet2, lambda);
681                 }
682         }
683         else if (ok1 && !ok2) {
684                 copy_v3_v3(meetco, meet1);
685         }
686         else if (!ok1 && ok2) {
687                 copy_v3_v3(meetco, meet2);
688         }
689         else {
690                 /* Neither offset line met emid.
691                  * This should only happen if all three lines are on top of each other */
692                 slide_dist(emid, v, e1->offset_r, meetco);
693         }
694
695         /* offsets may have changed now */
696         d = dist_to_line_v3(meetco, v->co, BM_edge_other_vert(e1->e, v)->co);
697         if (fabsf(d - e1->offset_r) > BEVEL_EPSILON)
698                 e1->offset_r = d;
699         d = dist_to_line_v3(meetco, v->co, BM_edge_other_vert(e2->e, v)->co);
700         if (fabsf(d - e2->offset_l) > BEVEL_EPSILON)
701                 e2->offset_l = d;
702 }
703
704 /* Calculate the best place for a meeting point for the offsets from edges e1 and e2
705  * when there is an in-between edge emid, and we prefer to have a point that may not
706  * be on emid if that does a better job of keeping offsets at the user spec.
707  * Viewed from the vertex normal side, the CCW order of the edges is e1, emid, e2.
708  * The offset lines may not meet exactly: the lines may be angled so that they can't meet.
709  * In that case, pick  the the offset_on_edge_between. */
710 static void offset_in_two_planes(BevelParams *bp, EdgeHalf *e1, EdgeHalf *e2, EdgeHalf *emid,
711                                  BMVert *v,  float meetco[3])
712 {
713         float dir1[3], dir2[3], dirmid[3], norm_perp1[3], norm_perp2[3],
714               off1a[3], off1b[3], off2a[3], off2b[3], isect2[3],
715               f1no[3], f2no[3], ang, d;
716         int iret;
717
718         /* get direction vectors for two offset lines */
719         sub_v3_v3v3(dir1, v->co, BM_edge_other_vert(e1->e, v)->co);
720         sub_v3_v3v3(dir2, BM_edge_other_vert(e2->e, v)->co, v->co);
721         sub_v3_v3v3(dirmid, BM_edge_other_vert(emid->e, v)->co, v->co);
722
723         /* get directions into offset planes */
724         /* calculate face normals at corner in case faces are nonplanar */
725         cross_v3_v3v3(f1no, dirmid, dir1);
726         cross_v3_v3v3(f2no, dirmid, dir2);
727
728         /* if e1-v-emid or emid-v-e2 are reflex angles, need to flip corner normals */
729         if (dot_v3v3(f1no, v->no) < 0.0f)
730                 negate_v3(f1no);
731         if (dot_v3v3(f2no, v->no) < 0.0f)
732                 negate_v3(f2no);
733
734         /* get vectors perpendicular to e1 and e2, pointing into the proper faces */
735         cross_v3_v3v3(norm_perp1, dir1, f1no);
736         normalize_v3(norm_perp1);
737         cross_v3_v3v3(norm_perp2, dir2, f2no);
738         normalize_v3(norm_perp2);
739
740         /* get points that are offset distances from each line, then another point on each line */
741         copy_v3_v3(off1a, v->co);
742         madd_v3_v3fl(off1a, norm_perp1, e1->offset_r);
743         sub_v3_v3v3(off1b, off1a, dir1);
744         copy_v3_v3(off2a, v->co);
745         madd_v3_v3fl(off2a, norm_perp2, e2->offset_l);
746         add_v3_v3v3(off2b, off2a, dir2);
747
748         ang = angle_v3v3(dir1, dir2);
749         if (ang < 100.0f * BEVEL_EPSILON) {
750                 /* lines are parallel; put intersection on emid */
751                 offset_on_edge_between(bp, e1, e2, emid, v, meetco);
752         }
753         else if (fabsf(ang - (float)M_PI) < 100.0f * BEVEL_EPSILON) {
754                 slide_dist(e2, v, e2->offset_l, meetco);
755                 d = dist_to_line_v3(meetco, v->co, BM_edge_other_vert(e1->e, v)->co);
756                 if (fabsf(d - e1->offset_r) > BEVEL_EPSILON)
757                         e1->offset_r = d;
758         }
759         else {
760                 iret = isect_line_line_v3(off1a, off1b, off2a, off2b, meetco, isect2);
761                 if (iret == 0) {
762                         /* lines colinear: another test says they are parallel. so shouldn't happen */
763                         copy_v3_v3(meetco, off1a);
764                         d = dist_to_line_v3(meetco, v->co, BM_edge_other_vert(e2->e, v)->co);
765                         if (fabsf(d - e2->offset_l) > BEVEL_EPSILON)
766                                 e2->offset_l = d;
767                 }
768                 else if (iret == 2) {
769                         /* lines are not coplanar and don't meet; meetco and isect2 are nearest to first and second lines */
770                         if (len_squared_v3v3(meetco, isect2) > 100.0f * BEVEL_EPSILON_SQ) {
771                                 /* offset lines don't meet so can't preserve widths */
772                                 offset_on_edge_between(bp, e1, e2, emid, v, meetco);
773                         }
774                 }
775                 /* else iret == 1 and the lines are coplanar so meetco has the intersection */
776         }
777 }
778
779 /* Offset by e->offset in plane with normal plane_no, on left if left==TRUE,
780  * else on right.  If no is NULL, choose an arbitrary plane different
781  * from eh's direction. */
782 static void offset_in_plane(EdgeHalf *e, const float plane_no[3], int left, float r[3])
783 {
784         float dir[3], no[3], fdir[3];
785         BMVert *v;
786
787         v = e->is_rev ? e->e->v2 : e->e->v1;
788
789         sub_v3_v3v3(dir, BM_edge_other_vert(e->e, v)->co, v->co);
790         normalize_v3(dir);
791         if (plane_no) {
792                 copy_v3_v3(no, plane_no);
793         }
794         else {
795                 zero_v3(no);
796                 if (fabsf(dir[0]) < fabsf(dir[1]))
797                         no[0] = 1.0f;
798                 else
799                         no[1] = 1.0f;
800         }
801         if (left)
802                 cross_v3_v3v3(fdir, dir, no);
803         else
804                 cross_v3_v3v3(fdir, no, dir);
805         normalize_v3(fdir);
806         copy_v3_v3(r, v->co);
807         madd_v3_v3fl(r, fdir, left? e->offset_l : e->offset_r);
808 }
809
810 /* Calculate the point on e where line (co_a, co_b) comes closest to and return it in projco */
811 static void project_to_edge(BMEdge *e, const float co_a[3], const float co_b[3], float projco[3])
812 {
813         float otherco[3];
814
815         if (!isect_line_line_v3(e->v1->co, e->v2->co, co_a, co_b, projco, otherco)) {
816 #ifdef BEVEL_ASSERT_PROJECT
817                 BLI_assert(!"project meet failure");
818 #endif
819                 copy_v3_v3(projco, e->v1->co);
820         }
821 }
822
823 /* If there is a bndv->ebev edge, find the mid control point if necessary.
824  * It is the closest point on the beveled edge to the line segment between
825  * bndv and bndv->next.  */
826 static void set_profile_params(BoundVert *bndv)
827 {
828         EdgeHalf *e;
829
830         e = bndv->ebev;
831         if (e) {
832                 project_to_edge(e->e, bndv->nv.co, bndv->next->nv.co,
833                                 bndv->profile.midco);
834         }
835 }
836
837 /* project co along the direction of edir to the plane containing
838  * e1 and e2 (which share a vert) */
839 static void project_to_edges_plane(float co[3], EdgeHalf *edir, EdgeHalf *e1, EdgeHalf *e2)
840 {
841         float dir[3], co2[3], d1[3], d2[3], no[3], snap[3];
842
843         sub_v3_v3v3(dir, edir->e->v1->co, edir->e->v2->co);
844         sub_v3_v3v3(d1, e1->e->v1->co, e1->e->v2->co);
845         sub_v3_v3v3(d2, e2->e->v1->co, e2->e->v2->co);
846         cross_v3_v3v3(no, d1, d2);
847         if (len_squared_v3(no) < BEVEL_EPSILON_SQ)
848                 return; /* e1 and e2 parallel so plane undefined -- don't snap */
849         add_v3_v3v3(co2, co, dir);
850         if (isect_line_plane_v3(snap, co, co2, e1->e->v1->co, no))
851                 copy_v3_v3(co, snap);
852 }
853
854 /* return 1 if a and b are in CCW order on the normal side of f,
855  * and -1 if they are reversed, and 0 if there is no shared face f */
856 static int bev_ccw_test(BMEdge *a, BMEdge *b, BMFace *f)
857 {
858         BMLoop *la, *lb;
859
860         if (!f)
861                 return 0;
862         la = BM_face_edge_share_loop(f, a);
863         lb = BM_face_edge_share_loop(f, b);
864         if (!la || !lb)
865                 return 0;
866         return lb->next == la ? 1 : -1;
867 }
868
869 /* Fill matrix r_mat so that a point in the sheared parallelogram with corners
870  * va, vmid, vb (and the 4th that is implied by it being a parallelogram)
871  * is the result of transforming the unit square by multiplication with r_mat.
872  * If it can't be done because the parallelogram is degenerate, return FALSE
873  * else return TRUE.
874  * Method:
875  * Find vo, the origin of the parallelogram with other three points va, vmid, vb.
876  * Also find vd, which is in direction normal to parallelogram and 1 unit away
877  * from the origin.
878  * The quarter circle in first quadrant of unit square will be mapped to the
879  * quadrant of a sheared ellipse in the parallelogram, using a matrix.
880  * The matrix mat is calculated to map:
881  *    (0,1,0) -> va
882  *    (1,1,0) -> vmid
883  *    (1,0,0) -> vb
884  *    (0,1,1) -> vd
885  * We want M to make M*A=B where A has the left side above, as columns
886  * and B has the right side as columns - both extended into homogeneous coords.
887  * So M = B*(Ainverse).  Doing Ainverse by hand gives the code below.
888  */
889 static int make_unit_square_map(const float va[3], const float vmid[3], const float vb[3],
890                                 float r_mat[4][4])
891 {
892         float vo[3], vd[3], vb_vmid[3], va_vmid[3], vddir[3];
893
894         sub_v3_v3v3(va_vmid, vmid, va);
895         sub_v3_v3v3(vb_vmid, vmid, vb);
896         if (fabsf(angle_v3v3(va_vmid, vb_vmid) - (float)M_PI) > 100.0f * BEVEL_EPSILON) {
897                 sub_v3_v3v3(vo, va, vb_vmid);
898                 cross_v3_v3v3(vddir, vb_vmid, va_vmid);
899                 normalize_v3(vddir);
900                 add_v3_v3v3(vd, vo, vddir);
901
902                 /* The cols of m are: {vmid - va, vmid - vb, vmid + vd - va -vb, va + vb - vmid;
903                  * blender transform matrices are stored such that m[i][*] is ith column;
904                  * the last elements of each col remain as they are in unity matrix */
905                 sub_v3_v3v3(&r_mat[0][0], vmid, va);
906                 r_mat[0][3] = 0.0f;
907                 sub_v3_v3v3(&r_mat[1][0], vmid, vb);
908                 r_mat[1][3] = 0.0f;
909                 add_v3_v3v3(&r_mat[2][0], vmid, vd);
910                 sub_v3_v3(&r_mat[2][0], va);
911                 sub_v3_v3(&r_mat[2][0], vb);
912                 r_mat[2][3] = 0.0f;
913                 add_v3_v3v3(&r_mat[3][0], va, vb);
914                 sub_v3_v3(&r_mat[3][0], vmid);
915                 r_mat[3][3] = 1.0f;
916
917                 return TRUE;
918         }
919         else
920                 return FALSE;
921 }
922
923 /* Like make_unit_square_map, but this one makes a matrix that transforms the
924  * (1,1,1) corner of a unit cube into an arbitrary corner with corner vert d
925  * and verts around it a, b, c (in ccw order, viewed from d normal dir).
926  * The matrix mat is calculated to map:
927  *    (1,0,0) -> va
928  *    (0,1,0) -> vb
929  *    (0,0,1) -> vc
930  *    (1,1,1) -> vd
931  * We want M to make M*A=B where A has the left side above, as columns
932  * and B has the right side as columns - both extended into homogeneous coords.
933  * So M = B*(Ainverse).  Doing Ainverse by hand gives the code below.
934  * The cols of M are 1/2{va-vb+vc-vd}, 1/2{-va+vb-vc+vd}, 1/2{-va-vb+vc+vd},
935  * and 1/2{va+vb+vc-vd}
936  * and Blender matrices have cols at m[i][*].
937  */
938 static void make_unit_cube_map(const float va[3], const float vb[3], const float vc[3],
939                               const float vd[3], float r_mat[4][4])
940 {
941         copy_v3_v3(r_mat[0], va);
942         sub_v3_v3(r_mat[0], vb);
943         sub_v3_v3(r_mat[0], vc);
944         add_v3_v3(r_mat[0], vd);
945         mul_v3_fl(r_mat[0], 0.5f);
946         r_mat[0][3] = 0.0f;
947         copy_v3_v3(r_mat[1], vb);
948         sub_v3_v3(r_mat[1], va);
949         sub_v3_v3(r_mat[1], vc);
950         add_v3_v3(r_mat[1], vd);
951         mul_v3_fl(r_mat[1], 0.5f);
952         r_mat[1][3] = 0.0f;
953         copy_v3_v3(r_mat[2], vc);
954         sub_v3_v3(r_mat[2], va);
955         sub_v3_v3(r_mat[2], vb);
956         add_v3_v3(r_mat[2], vd);
957         mul_v3_fl(r_mat[2], 0.5f);
958         r_mat[2][3] = 0.0f;
959         copy_v3_v3(r_mat[3], va);
960         add_v3_v3(r_mat[3], vb);
961         add_v3_v3(r_mat[3], vc);
962         sub_v3_v3(r_mat[3], vd);
963         mul_v3_fl(r_mat[3], 0.5f);
964         r_mat[3][3] = 1.0f;
965 }
966
967 #ifndef USE_ADJ_SUBDIV
968 /*
969  * Find the point (/n) of the way around the round profile for e,
970  * where start point is va, midarc point is vmid, and end point is vb.
971  * Return the answer in r_co.
972  * If va -- vmid -- vb is approximately a straight line, just
973  * interpolate along the line.
974  */
975 static void get_point_on_round_edge(EdgeHalf *e, int k,
976                                     const float va[3], const float vmid[3], const float vb[3],
977                                     float r_co[3])
978 {
979         float p[3], angle;
980         float m[4][4];
981         int n = e->seg;
982
983         if (make_unit_square_map(va, vmid, vb, m)) {
984                 /* Find point k/(e->seg) along quarter circle from (0,1,0) to (1,0,0) */
985                 angle = (float)M_PI * (float)k / (2.0f * (float)n);  /* angle from y axis */
986                 p[0] = sinf(angle);
987                 p[1] = cosf(angle);
988                 p[2] = 0.0f;
989                 mul_v3_m4v3(r_co, m, p);
990         }
991         else {
992                 /* degenerate case: profile is a line */
993                 interp_v3_v3v3(r_co, va, vb, (float)k / (float)n);
994         }
995 }
996 #endif
997
998 /* Find the point on given profile at parameter u which goes from 0 to 2 as
999  * the profile is moved from va to vb. */
1000 static void get_profile_point(const Profile *pro, const float va[3], const float vb[3], float u, float r_co[3])
1001  {
1002         float p[3], angle;
1003         float m[4][4];
1004
1005         if (u <= 0.0f)
1006                 copy_v3_v3(r_co, va);
1007         else if (u >= 2.0f)
1008                 copy_v3_v3(r_co, vb);
1009         else if (!make_unit_square_map(va, pro->midco, vb, m)) {
1010                 interp_v3_v3v3(r_co, va, vb, u / 2.0f);
1011         }
1012         else {
1013                 angle = u * (float)M_PI / 4.0f;  /* angle from y axis */
1014                 p[0] = sinf(angle);
1015                 p[1] = cosf(angle);
1016                 p[2] = 0.0f;
1017                 mul_v3_m4v3(r_co, m, p);
1018         }
1019 }
1020
1021 BLI_INLINE void get_bndv_profile_point(const BoundVert *bndv, const float u, float r_co[3])
1022 {
1023         get_profile_point(&bndv->profile, bndv->nv.co, bndv->next->nv.co, u, r_co);
1024 }
1025
1026 #ifndef USE_ADJ_SUBDIV
1027 /* Calculate a snapped point to the transformed profile of edge e, extended as
1028  * in a cylinder-like surface in the direction of e.
1029  * co is the point to snap and is modified in place.
1030  * va and vb are the limits of the profile (with peak on e). */
1031 static void snap_to_edge_profile(EdgeHalf *e, const float va[3], const float vb[3],
1032                                  float co[3])
1033 {
1034         float m[4][4], minv[4][4];
1035         float edir[3], va0[3], vb0[3], vmid0[3], p[3], snap[3], plane[4];
1036
1037         sub_v3_v3v3(edir, e->e->v1->co, e->e->v2->co);
1038
1039         /* project va and vb onto plane P, with normal edir and containing co */
1040         plane_from_point_normal_v3(plane, co, edir);
1041         closest_to_plane_v3(va0, plane, va);
1042         closest_to_plane_v3(vb0, plane, vb);
1043         project_to_edge(e->e, va0, vb0, vmid0);
1044         if (make_unit_square_map(va0, vmid0, vb0, m)) {
1045                 /* Transform co and project it onto the unit circle.
1046                  * Projecting is in fact just normalizing the transformed co */
1047                 if (!invert_m4_m4(minv, m)) {
1048                         /* shouldn't happen, by angle test and construction of vd */
1049                         BLI_assert(!"failed inverse during profile snap");
1050                         return;
1051                 }
1052                 mul_v3_m4v3(p, minv, co);
1053                 normalize_v3(p);
1054                 mul_v3_m4v3(snap, m, p);
1055                 copy_v3_v3(co, snap);
1056         }
1057         else {
1058                 /* planar case: just snap to line va--vb */
1059                 closest_to_line_segment_v3(p, co, va, vb);
1060                 copy_v3_v3(co, p);
1061         }
1062 }
1063 #endif
1064
1065 /* Snapping a direction co to a unit radius sphere is just normalizing co.
1066  * TODO: generalize to superellipsoid */
1067 static void snap_to_sphere(float co[3])
1068 {
1069         normalize_v3(co);
1070 }
1071
1072 static void snap_to_profile(BoundVert *bndv, EdgeHalf *e, float co[3])
1073 {
1074         float va[3], vb[3], edir[3], va0[3], vb0[3], vmid0[3];
1075         float plane[4], m[4][4], minv[4][4], p[3], snap[3];
1076
1077         copy_v3_v3(va, bndv->nv.co);
1078         copy_v3_v3(vb, bndv->next->nv.co);
1079
1080         sub_v3_v3v3(edir, e->e->v1->co, e->e->v2->co);
1081
1082         plane_from_point_normal_v3(plane, co, edir);
1083         closest_to_plane_v3(va0, plane, va);
1084         closest_to_plane_v3(vb0, plane, vb);
1085         closest_to_plane_v3(vmid0, plane, bndv->profile.midco);
1086         if (make_unit_square_map(va0, vmid0, vb0, m)) {
1087                 /* Transform co and project it onto sphere */
1088                 if (!invert_m4_m4(minv, m)) {
1089                         /* shouldn't happen */
1090                         BLI_assert(!"failed inverse during profile snap");
1091                         return;
1092                 }
1093                 mul_v3_m4v3(p, minv, co);
1094                 snap_to_sphere(p);
1095                 mul_v3_m4v3(snap, m, p);
1096                 copy_v3_v3(co, snap);
1097         }
1098         else {
1099                 /* planar case: just snap to line va--vb */
1100                 closest_to_line_segment_v3(p, co, va, vb);
1101                 copy_v3_v3(co, p);
1102         }
1103 }
1104
1105 /* Set the any_seam property for a BevVert and all its BoundVerts */
1106 static void set_bound_vert_seams(BevVert *bv)
1107 {
1108         BoundVert *v;
1109         EdgeHalf *e;
1110
1111         bv->any_seam = false;
1112         v = bv->vmesh->boundstart;
1113         do {
1114                 v->any_seam = false;
1115                 for (e = v->efirst; e; e = e->next) {
1116                         v->any_seam |= e->is_seam;
1117                         if (e == v->elast)
1118                                 break;
1119                 }
1120                 bv->any_seam |= v->any_seam;
1121         } while ((v = v->next) != bv->vmesh->boundstart);
1122 }
1123
1124 /* Make a circular list of BoundVerts for bv, each of which has the coordinates
1125  * of a vertex on the the boundary of the beveled vertex bv->v.
1126  * This may adjust some EdgeHalf widths, and there might have to be
1127  * a subsequent pass to make the widths as consistent as possible.
1128  * The first time through, construct will be true and we are making the BoundVerts
1129  * and setting up the BoundVert and EdgeHalf pointers appropriately.
1130  * For a width consistency path, we just recalculate the coordinates of the
1131  * BoundVerts. If the other ends have been (re)built already, then we
1132  * copy the offsets from there to match, else we use the ideal (user-specified)
1133  * widths.
1134  * Also, if construct, decide on the mesh pattern that will be used inside the boundary.
1135  * Doesn't make the actual BMVerts */
1136 static void build_boundary(BevelParams *bp, BevVert *bv, bool construct)
1137 {
1138         MemArena *mem_arena = bp->mem_arena;
1139         EdgeHalf *efirst, *e, *eother;
1140         BoundVert *v;
1141         BevVert *bvother;
1142         VMesh *vm;
1143         float co[3];
1144         const float  *no;
1145         float lastd;
1146
1147         vm = bv->vmesh;
1148
1149         if (bp->vertex_only) {
1150                 e = efirst = &bv->edges[0];
1151         }
1152         else {
1153                 e = efirst = next_bev(bv, NULL);
1154                 do {
1155                         eother = find_other_end_edge_half(bp, e, &bvother);
1156                         if (eother && bvother->visited && bp->offset_type != BEVEL_AMT_PERCENT) {
1157                                 /* try to keep bevel even by matching other end offsets */
1158                                 e->offset_l = eother->offset_r;
1159                                 e->offset_r = eother->offset_l;
1160                         }
1161                         else {
1162                                 /* reset to user spec */
1163                                 e->offset_l = e->offset_l_spec;
1164                                 e->offset_r = e->offset_r_spec;
1165                         }
1166                 } while ((e = e->next) != efirst);
1167                 e = efirst;
1168         }
1169
1170         BLI_assert(bv->edgecount >= 2);  /* since bevel edges incident to 2 faces */
1171
1172         if (bv->edgecount == 2 && bv->selcount == 1) {
1173                 /* special case: beveled edge meets non-beveled one at valence 2 vert */
1174                 no = e->fprev ? e->fprev->no : (e->fnext ? e->fnext->no : NULL);
1175                 offset_in_plane(e, no, TRUE, co);
1176                 if (construct) {
1177                         v = add_new_bound_vert(mem_arena, vm, co);
1178                         v->efirst = v->elast = v->ebev = e;
1179                         e->leftv = v;
1180                 }
1181                 else {
1182                         adjust_bound_vert(e->leftv, co);
1183                 }
1184                 no = e->fnext ? e->fnext->no : (e->fprev ? e->fprev->no : NULL);
1185                 offset_in_plane(e, no, FALSE, co);
1186                 if (construct) {
1187                         v = add_new_bound_vert(mem_arena, vm, co);
1188                         v->efirst = v->elast = e;
1189                         e->rightv = v;
1190                 }
1191                 else {
1192                         adjust_bound_vert(e->rightv, co);
1193                 }
1194                 /* make artifical extra point along unbeveled edge, and form triangle */
1195                 slide_dist(e->next, bv->v, e->offset_l, co);
1196                 if (construct) {
1197                         v = add_new_bound_vert(mem_arena, vm, co);
1198                         v->efirst = v->elast = e->next;
1199                         e->next->leftv = e->next->rightv = v;
1200                         /* could use M_POLY too, but tri-fan looks nicer)*/
1201                         vm->mesh_kind = M_TRI_FAN;
1202                         set_bound_vert_seams(bv);
1203                 }
1204                 else {
1205                         adjust_bound_vert(e->next->leftv, co);
1206                 }
1207                 set_profile_params(vm->boundstart);
1208                 return;
1209         }
1210
1211         lastd = bp->vertex_only ? bv->offset : e->offset_l;
1212         do {
1213                 if (e->is_bev) {
1214                         /* handle only left side of beveled edge e here: next iteration should do right side */
1215                         if (e->prev->is_bev) {
1216                                 BLI_assert(e->prev != e);  /* see: wire edge special case */
1217                                 offset_meet(e->prev, e, bv->v, e->fprev, co);
1218                                 if (construct) {
1219                                         v = add_new_bound_vert(mem_arena, vm, co);
1220                                         v->efirst = e->prev;
1221                                         v->elast = v->ebev = e;
1222                                         e->leftv = v;
1223                                         e->prev->rightv = v;
1224                                 }
1225                                 else {
1226                                         v = e->leftv;
1227                                         adjust_bound_vert(v, co);
1228                                 }
1229                         }
1230                         else {
1231                                 /* e->prev is not beveled */
1232                                 if (e->prev->prev->is_bev) {
1233                                         BLI_assert(e->prev->prev != e); /* see: edgecount 2, selcount 1 case */
1234                                         /* find meet point between e->prev->prev and e and attach e->prev there */
1235                                         if (bp->preserve_widths)
1236                                                 offset_in_two_planes(bp, e->prev->prev, e, e->prev, bv->v, co);
1237                                         else
1238                                                 offset_on_edge_between(bp, e->prev->prev, e, e->prev, bv->v, co);
1239                                         if (construct) {
1240                                                 v = add_new_bound_vert(mem_arena, vm, co);
1241                                                 v->efirst = e->prev->prev;
1242                                                 v->elast = v->ebev = e;
1243                                                 e->leftv = v;
1244                                                 e->prev->leftv = v;
1245                                                 e->prev->prev->rightv = v;
1246                                         }
1247                                         else {
1248                                                 v = e->leftv;
1249                                                 adjust_bound_vert(v, co);
1250                                         }
1251                                 }
1252                                 else {
1253                                         /* neither e->prev nor e->prev->prev are beveled: make on-edge on e->prev */
1254                                         offset_meet(e->prev, e, bv->v, e->fprev, co);
1255                                         if (construct) {
1256                                                 v = add_new_bound_vert(mem_arena, vm, co);
1257                                                 v->efirst = e->prev;
1258                                                 v->elast = v->ebev = e;
1259                                                 e->leftv = v;
1260                                                 e->prev->leftv = v;
1261                                         }
1262                                         else {
1263                                                 v = e->leftv;
1264                                                 adjust_bound_vert(v, co);
1265                                         }
1266                                 }
1267                         }
1268                         lastd = len_v3v3(bv->v->co, v->nv.co);
1269                 }
1270                 else {
1271                         /* e is not beveled */
1272                         if (e->next->is_bev) {
1273                                 /* next iteration will place e between beveled previous and next edges */
1274                                 /* do nothing... */
1275                         }
1276                         else if (e->prev->is_bev) {
1277                                 /* on-edge meet between e->prev and e */
1278                                 offset_meet(e->prev, e, bv->v, e->fprev, co);
1279                                 if (construct) {
1280                                         v = add_new_bound_vert(mem_arena, vm, co);
1281                                         v->efirst = e->prev;
1282                                         v->elast = e;
1283                                         e->leftv = v;
1284                                         e->prev->rightv = v;
1285                                 }
1286                                 else {
1287                                         adjust_bound_vert(e->leftv, co);
1288                                 }
1289                         }
1290                         else {
1291                                 /* None of e->prev, e, e->next are beveled.
1292                                  * could either leave alone or add slide points to make
1293                                  * one polygon around bv->v.  For now, we choose latter.
1294                                  * Could slide to make an even bevel plane but for now will
1295                                  * just use last distance a meet point moved from bv->v. */
1296                                 slide_dist(e, bv->v, lastd, co);
1297                                 if (construct) {
1298                                         v = add_new_bound_vert(mem_arena, vm, co);
1299                                         v->efirst = v->elast = e;
1300                                         e->leftv = v;
1301                                 }
1302                                 else {
1303                                         adjust_bound_vert(e->leftv, co);
1304                                 }
1305                         }
1306                 }
1307         } while ((e = e->next) != efirst);
1308
1309         v = vm->boundstart;
1310         do {
1311                 set_profile_params(v);
1312         } while ((v = v->next) != vm->boundstart);
1313
1314         if (bv->selcount == 1 && bv->edgecount == 3) {
1315                 /* special case: snap profile to third face */
1316                 v = vm->boundstart;
1317                 BLI_assert(v->ebev != NULL);
1318                 project_to_edges_plane(v->profile.midco, v->ebev, v->efirst, v->next->elast);
1319         }
1320
1321         if (construct) {
1322                 set_bound_vert_seams(bv);
1323
1324                 BLI_assert(vm->count >= 2);
1325                 if (bp->vertex_only) {
1326                         if (vm->count == 2)
1327                                 vm->mesh_kind = M_NONE;
1328                         else if (bp->seg > 1)
1329                                 vm->mesh_kind = M_ADJ_SUBDIV;
1330                         else
1331                                 vm->mesh_kind = M_POLY;
1332                 }
1333                 else if (vm->count == 2 && bv->edgecount == 3) {
1334                         vm->mesh_kind = M_NONE;
1335                 }
1336                 else if (bv->selcount == 2) {
1337                         vm->mesh_kind = M_QUAD_STRIP;
1338                 }
1339                 else if (efirst->seg == 1 || bv->selcount == 1) {
1340                         if (vm->count == 3 && bv->selcount == 1) {
1341                                 vm->mesh_kind = M_TRI_FAN;
1342                         }
1343                         else {
1344                                 vm->mesh_kind = M_POLY;
1345                         }
1346                 }
1347                 else {
1348 #ifdef USE_ADJ_SUBDIV
1349                         vm->mesh_kind = M_ADJ_SUBDIV;
1350 #else
1351                         vm->mesh_kind = M_ADJ;
1352 #endif
1353                 }
1354         }
1355 }
1356
1357 /* Do a global pass to try to make offsets as even as possible.
1358  * Consider this graph:
1359  *   nodes = BevVerts
1360  *   edges = { (u,v) } where u and v are nodes such that u and v
1361  *        are connected by a mesh edge that has at least one end
1362  *        whose offset does not match the user spec.
1363  *
1364  * Do a breadth-first search on this graph, starting from nodes
1365  * that have any_adjust=true, and changing all
1366  * not-already-changed offsets on EdgeHalfs to match the
1367  * corresponding ones that changed on the other end.
1368  * The graph is dynamic in the sense that having an offset that
1369  * doesn't meet the user spec can be added as the search proceeds.
1370  * We want this search to be deterministic (not dependendent
1371  * on order of processing through hash table), so as to avoid
1372  * flicker to to different decisions made if search is different
1373  * while dragging the offset number in the UI.  So look for the
1374  * lower vertex number when there is a choice of where to start.
1375  *
1376  * Note that this might not process all BevVerts, only the ones
1377  * that need adjustment.
1378  */
1379 static void adjust_offsets(BevelParams *bp)
1380 {
1381         BevVert *bv, *searchbv, *bvother;
1382         int i, searchi;
1383         GHashIterator giter;
1384         EdgeHalf *e, *efirst, *eother;
1385         GSQueue *q;
1386
1387         BLI_assert(!bp->vertex_only);
1388         GHASH_ITER(giter, bp->vert_hash) {
1389                 bv = BLI_ghashIterator_getValue(&giter);
1390                 bv->visited = false;
1391         }
1392
1393         q = BLI_gsqueue_new((int)sizeof(BevVert *));
1394         /* the following loop terminates because at least one node is visited each time */
1395         for (;;) {
1396                 /* look for root of a connected component in search graph */
1397                 searchbv = NULL;
1398                 searchi = -1;
1399                 GHASH_ITER(giter, bp->vert_hash) {
1400                         bv = BLI_ghashIterator_getValue(&giter);
1401                         if (!bv->visited && any_edge_half_offset_changed(bv)) {
1402                                 i = BM_elem_index_get(bv->v);
1403                                 if (!searchbv || i < searchi) {
1404                                         searchbv = bv;
1405                                         searchi = i;
1406                                 }
1407                         }
1408                 }
1409                 if (searchbv == NULL)
1410                         break;
1411
1412                 BLI_gsqueue_push(q, &searchbv);
1413                 while (!BLI_gsqueue_is_empty(q)) {
1414                         BLI_gsqueue_pop(q, &bv);
1415                         /* If do this check, don't have to check for already-on-queue before push, below */
1416                         if (bv->visited)
1417                                 continue;
1418                         bv->visited = true;
1419                         build_boundary(bp, bv, false);
1420
1421                         e = efirst = &bv->edges[0];
1422                         do {
1423                                 eother = find_other_end_edge_half(bp, e, &bvother);
1424                                 if (eother && !bvother->visited && edge_half_offset_changed(e)) {
1425                                         BLI_gsqueue_push(q, &bvother);
1426                                 }
1427                         } while ((e = e->next) != efirst);
1428                 }
1429         }
1430         BLI_gsqueue_free(q);
1431 }
1432
1433 /* Do the edges at bv form a "pipe"?
1434  * Current definition: at least three beveled edges,
1435  * two in line, and sharing a face. */
1436 static EdgeHalf *pipe_test(BevVert *bv)
1437 {
1438         EdgeHalf *e1, *e2, *epipe;
1439
1440         epipe = NULL;
1441         if (bv->selcount > 2) {
1442                 for (e1 = &bv->edges[0]; epipe == NULL && e1 != &bv->edges[bv->edgecount]; e1++) {
1443                         if (e1->is_bev) {
1444                                 for (e2 = &bv->edges[0]; e2 != &bv->edges[bv->edgecount]; e2++) {
1445                                         if (e1 != e2 && e2->is_bev) {
1446                                                 if ((e1->fnext == e2->fprev) || (e1->fprev == e2->fnext)) {
1447                                                         float dir1[3], dir2[3];
1448                                                         sub_v3_v3v3(dir1, bv->v->co, BM_edge_other_vert(e1->e, bv->v)->co);
1449                                                         sub_v3_v3v3(dir2, BM_edge_other_vert(e2->e, bv->v)->co, bv->v->co);
1450                                                         if (angle_v3v3(dir1, dir2) < 100.0f * BEVEL_EPSILON) {
1451                                                                 epipe = e1;
1452                                                                 break;
1453                                                         }
1454                                                 }
1455                                         }
1456                                 }
1457                         }
1458                 }
1459         }
1460         return epipe;
1461 }
1462
1463 #ifndef USE_ADJ_SUBDIV
1464 /*
1465  * Given that the boundary is built and the boundary BMVerts have been made,
1466  * calculate the positions of the interior mesh points for the M_ADJ pattern,
1467  * then make the BMVerts and the new faces. */
1468 static void bevel_build_rings(BMesh *bm, BevVert *bv)
1469 {
1470         int k, ring, i, n, ns, ns2, nn, odd;
1471         VMesh *vm = bv->vmesh;
1472         BoundVert *v, *vprev, *vnext;
1473         NewVert *nv, *nvprev, *nvnext;
1474         EdgeHalf *epipe;
1475         BMVert *bmv, *bmv1, *bmv2, *bmv3, *bmv4;
1476         BMFace *f, *f2, *f23;
1477         float co[3], coa[3], cob[3], midco[3];
1478         float va_pipe[3], vb_pipe[3];
1479
1480         n = vm->count;
1481         ns = vm->seg;
1482         ns2 = ns / 2;
1483         odd = (ns % 2) != 0;
1484         BLI_assert(n > 2 && ns > 1);
1485
1486         epipe = pipe_test(bv);
1487
1488         /* Make initial rings, going between points on neighbors.
1489          * After this loop, will have coords for all (i, r, k) where
1490          * BoundVert for i has a bevel, 0 <= r <= ns2, 0 <= k <= ns */
1491         for (ring = 1; ring <= ns2; ring++) {
1492                 v = vm->boundstart;
1493
1494                 do {
1495                         i = v->index;
1496                         if (v->ebev) {
1497                                 /* get points coords of points a and b, on outer rings
1498                                  * of prev and next edges, k away from this edge */
1499                                 vprev = v->prev;
1500                                 vnext = v->next;
1501
1502                                 if (vprev->ebev)
1503                                         nvprev = mesh_vert(vm, vprev->index, 0, ns - ring);
1504                                 else
1505                                         nvprev = mesh_vert(vm, vprev->index, 0, ns);
1506                                 copy_v3_v3(coa, nvprev->co);
1507                                 nv = mesh_vert(vm, i, ring, 0);
1508                                 copy_v3_v3(nv->co, coa);
1509                                 nv->v = nvprev->v;
1510
1511                                 if (vnext->ebev)
1512                                         nvnext = mesh_vert(vm, vnext->index, 0, ring);
1513                                 else
1514                                         nvnext = mesh_vert(vm, vnext->index, 0, 0);
1515                                 copy_v3_v3(cob, nvnext->co);
1516                                 nv = mesh_vert(vm, i, ring, ns);
1517                                 copy_v3_v3(nv->co, cob);
1518                                 nv->v = nvnext->v;
1519
1520                                 /* TODO: better calculation of new midarc point? */
1521                                 project_to_edge(v->ebev->e, coa, cob, midco);
1522
1523                                 for (k = 1; k < ns; k++) {
1524                                         get_point_on_round_edge(v->ebev, k, coa, midco, cob, co);
1525                                         copy_v3_v3(mesh_vert(vm, i, ring, k)->co, co);
1526                                 }
1527
1528                                 if (v->ebev == epipe) {
1529                                         /* save profile extremes for later snapping */
1530                                         copy_v3_v3(va_pipe, mesh_vert(vm, i, 0, 0)->co);
1531                                         copy_v3_v3(vb_pipe, mesh_vert(vm, i, 0, ns)->co);
1532                                 }
1533                         }
1534                 } while ((v = v->next) != vm->boundstart);
1535         }
1536
1537         /* Now make sure cross points of rings share coordinates and vertices.
1538          * After this loop, will have BMVerts for all (i, r, k) where
1539          * i is for a BoundVert that is beveled and has either a predecessor or
1540          * successor BoundVert beveled too, and
1541          * for odd ns: 0 <= r <= ns2, 0 <= k <= ns
1542          * for even ns: 0 <= r < ns2, 0 <= k <= ns except k=ns2 */
1543         v = vm->boundstart;
1544         do {
1545                 i = v->index;
1546                 if (v->ebev) {
1547                         vprev = v->prev;
1548                         vnext = v->next;
1549                         if (vprev->ebev) {
1550                                 for (ring = 1; ring <= ns2; ring++) {
1551                                         for (k = 1; k <= ns2; k++) {
1552                                                 if (!odd && (k == ns2 || ring == ns2))
1553                                                         continue;  /* center line is special case: do after the rest are done */
1554                                                 nv = mesh_vert(vm, i, ring, k);
1555                                                 nvprev = mesh_vert(vm, vprev->index, k, ns - ring);
1556                                                 mid_v3_v3v3(co, nv->co, nvprev->co);
1557                                                 if (epipe)
1558                                                         snap_to_edge_profile(epipe, va_pipe, vb_pipe, co);
1559
1560                                                 copy_v3_v3(nv->co, co);
1561                                                 BLI_assert(nv->v == NULL && nvprev->v == NULL);
1562                                                 create_mesh_bmvert(bm, vm, i, ring, k, bv->v);
1563                                                 copy_mesh_vert(vm, vprev->index, k, ns - ring, i, ring, k);
1564                                         }
1565                                 }
1566                                 if (!vprev->prev->ebev) {
1567                                         for (ring = 1; ring <= ns2; ring++) {
1568                                                 for (k = 1; k <= ns2; k++) {
1569                                                         if (!odd && (k == ns2 || ring == ns2))
1570                                                                 continue;
1571                                                         create_mesh_bmvert(bm, vm, vprev->index, ring, k, bv->v);
1572                                                 }
1573                                         }
1574                                 }
1575                                 if (!vnext->ebev) {
1576                                         for (ring = 1; ring <= ns2; ring++) {
1577                                                 for (k = ns - ns2; k < ns; k++) {
1578                                                         if (!odd && (k == ns2 || ring == ns2))
1579                                                                 continue;
1580                                                         create_mesh_bmvert(bm, vm, i, ring, k, bv->v);
1581                                                 }
1582                                         }
1583                                 }
1584                         }
1585                 }
1586         } while ((v = v->next) != vm->boundstart);
1587
1588         if (!odd) {
1589                 /* Do special case center lines.
1590                  * This loop makes verts for (i, ns2, k) for 1 <= k <= ns-1, k!=ns2
1591                  * and for (i, r, ns2) for 1 <= r <= ns2-1,
1592                  * whenever i is in a sequence of at least two beveled verts */
1593                 v = vm->boundstart;
1594                 do {
1595                         i = v->index;
1596                         if (v->ebev) {
1597                                 vprev = v->prev;
1598                                 vnext = v->next;
1599                                 for (k = 1; k < ns2; k++) {
1600                                         nv = mesh_vert(vm, i, k, ns2);
1601                                         if (vprev->ebev)
1602                                                 nvprev = mesh_vert(vm, vprev->index, ns2, ns - k);
1603                                         if (vnext->ebev)
1604                                                 nvnext = mesh_vert(vm, vnext->index, ns2, k);
1605                                         if (vprev->ebev && vnext->ebev) {
1606                                                 mid_v3_v3v3v3(co, nvprev->co, nv->co, nvnext->co);
1607                                                 if (epipe)
1608                                                         snap_to_edge_profile(epipe, va_pipe, vb_pipe, co);
1609                                                 copy_v3_v3(nv->co, co);
1610                                                 create_mesh_bmvert(bm, vm, i, k, ns2, bv->v);
1611                                                 copy_mesh_vert(vm, vprev->index, ns2, ns - k, i, k, ns2);
1612                                                 copy_mesh_vert(vm, vnext->index, ns2, k, i, k, ns2);
1613
1614                                         }
1615                                         else if (vprev->ebev) {
1616                                                 mid_v3_v3v3(co, nvprev->co, nv->co);
1617                                                 if (epipe)
1618                                                         snap_to_edge_profile(epipe, va_pipe, vb_pipe, co);
1619                                                 copy_v3_v3(nv->co, co);
1620                                                 create_mesh_bmvert(bm, vm, i, k, ns2, bv->v);
1621                                                 copy_mesh_vert(vm, vprev->index, ns2, ns - k, i, k, ns2);
1622
1623                                                 create_mesh_bmvert(bm, vm, i, ns2, ns - k, bv->v);
1624                                         }
1625                                         else if (vnext->ebev) {
1626                                                 mid_v3_v3v3(co, nv->co, nvnext->co);
1627                                                 if (epipe)
1628                                                         snap_to_edge_profile(epipe, va_pipe, vb_pipe, co);
1629                                                 copy_v3_v3(nv->co, co);
1630                                                 create_mesh_bmvert(bm, vm, i, k, ns2, bv->v);
1631                                                 copy_mesh_vert(vm, vnext->index, ns2, k, i, k, ns2);
1632
1633                                                 create_mesh_bmvert(bm, vm, i, ns2, k, bv->v);
1634                                         }
1635                                 }
1636                         }
1637                 } while ((v = v->next) != vm->boundstart);
1638
1639                 /* center point need to be average of all centers of rings */
1640                 /* TODO: this is wrong if not all verts have ebev: could have
1641                  * several disconnected sections of mesh. */
1642                 zero_v3(midco);
1643                 nn = 0;
1644                 v = vm->boundstart;
1645                 do {
1646                         i = v->index;
1647                         if (v->ebev) {
1648                                 nv = mesh_vert(vm, i, ns2, ns2);
1649                                 add_v3_v3(midco, nv->co);
1650                                 nn++;
1651                         }
1652                 } while ((v = v->next) != vm->boundstart);
1653                 mul_v3_fl(midco, 1.0f / nn);
1654                 if (epipe)
1655                         snap_to_edge_profile(epipe, va_pipe, vb_pipe, midco);
1656                 bmv = BM_vert_create(bm, midco, NULL, BM_CREATE_NOP);
1657                 v = vm->boundstart;
1658                 do {
1659                         i = v->index;
1660                         if (v->ebev) {
1661                                 nv = mesh_vert(vm, i, ns2, ns2);
1662                                 copy_v3_v3(nv->co, midco);
1663                                 nv->v = bmv;
1664                         }
1665                 } while ((v = v->next) != vm->boundstart);
1666         }
1667
1668         /* Make the ring quads */
1669         for (ring = 0; ring < ns2; ring++) {
1670                 v = vm->boundstart;
1671                 do {
1672                         i = v->index;
1673                         f = boundvert_rep_face(v);
1674                         f2 = boundvert_rep_face(v->next);
1675                         if (v->ebev && (v->prev->ebev || v->next->ebev)) {
1676                                 for (k = 0; k < ns2 + odd; k++) {
1677                                         bmv1 = mesh_vert(vm, i, ring, k)->v;
1678                                         bmv2 = mesh_vert(vm, i, ring, k + 1)->v;
1679                                         bmv3 = mesh_vert(vm, i, ring + 1, k + 1)->v;
1680                                         bmv4 = mesh_vert(vm, i, ring + 1, k)->v;
1681                                         BLI_assert(bmv1 && bmv2 && bmv3 && bmv4);
1682                                         if (bmv3 == bmv4 || bmv1 == bmv4)
1683                                                 bmv4 = NULL;
1684                                         /* f23 is interp face for bmv2 and bmv3 */
1685                                         f23 = f;
1686                                         if (odd && k == ns2 && f2 && !v->any_seam)
1687                                                 f23 = f2;
1688                                         bev_create_quad_tri_ex(bm, bmv1, bmv2, bmv3, bmv4,
1689                                                                f, f23, f23, f);
1690                                 }
1691                         }
1692                         else if (v->prev->ebev && v->prev->prev->ebev) {
1693                                 /* finish off a sequence of beveled edges */
1694                                 i = v->prev->index;
1695                                 f = boundvert_rep_face(v->prev);
1696                                 f2 = boundvert_rep_face(v);
1697                                 for (k = ns2 + (ns % 2); k < ns; k++) {
1698                                         bmv1 = mesh_vert(vm, i, ring, k)->v;
1699                                         bmv2 = mesh_vert(vm, i, ring, k + 1)->v;
1700                                         bmv3 = mesh_vert(vm, i, ring + 1, k + 1)->v;
1701                                         bmv4 = mesh_vert(vm, i, ring + 1, k)->v;
1702                                         BLI_assert(bmv1 && bmv2 && bmv3 && bmv4);
1703                                         if (bmv2 == bmv3) {
1704                                                 bmv3 = bmv4;
1705                                                 bmv4 = NULL;
1706                                         }
1707                                         f23 = f;
1708                                         if (odd && k == ns2 && f2 && !v->any_seam)
1709                                                 f23 = f2;
1710                                         bev_create_quad_tri_ex(bm, bmv1, bmv2, bmv3, bmv4,
1711                                                                f, f23, f23, f);
1712                                 }
1713                         }
1714                 } while ((v = v->next) != vm->boundstart);
1715         }
1716
1717         /* Fix UVs along center lines if even number of segments */
1718         if (!odd) {
1719                 v = vm->boundstart;
1720                 do {
1721                         i = v->index;
1722                         f = boundvert_rep_face(v);
1723                         f2 = boundvert_rep_face(v->next);
1724                         if (!v->any_seam) {
1725                                 for (ring = 1; ring < ns2; ring++) {
1726                                         BMVert *v_uv = mesh_vert(vm, i, ring, ns2)->v;
1727                                         if (v_uv) {
1728                                                 bev_merge_uvs(bm, v_uv);
1729                                         }
1730                                 }
1731                         }
1732                 } while ((v = v->next) != vm->boundstart);
1733                 if (!bv->any_seam)
1734                         bev_merge_uvs(bm, mesh_vert(vm, 0, ns2, ns2)->v);
1735         }
1736
1737         /* Make center ngon if odd number of segments and fully beveled */
1738         if (odd && vm->count == bv->selcount) {
1739                 BMVert **vv = NULL;
1740                 BMFace **vf = NULL;
1741                 BLI_array_staticdeclare(vv, BM_DEFAULT_NGON_STACK_SIZE);
1742                 BLI_array_staticdeclare(vf, BM_DEFAULT_NGON_STACK_SIZE);
1743
1744                 v = vm->boundstart;
1745                 do {
1746                         i = v->index;
1747                         BLI_assert(v->ebev);
1748                         BLI_array_append(vv, mesh_vert(vm, i, ns2, ns2)->v);
1749                         BLI_array_append(vf, bv->any_seam ? f: boundvert_rep_face(v));
1750                 } while ((v = v->next) != vm->boundstart);
1751                 f = boundvert_rep_face(vm->boundstart);
1752                 bev_create_ngon(bm, vv, BLI_array_count(vv), vf, f, true);
1753
1754                 BLI_array_free(vv);
1755         }
1756
1757         /* Make 'rest-of-vmesh' polygon if not fully beveled */
1758         /* TODO: use interpolation face array here too */
1759         if (vm->count > bv->selcount) {
1760                 int j;
1761                 BMVert **vv = NULL;
1762                 BLI_array_staticdeclare(vv, BM_DEFAULT_NGON_STACK_SIZE);
1763
1764                 v = vm->boundstart;
1765                 f = boundvert_rep_face(v);
1766                 j = 0;
1767                 do {
1768                         i = v->index;
1769                         if (v->ebev) {
1770                                 if (!v->prev->ebev) {
1771                                         for (k = 0; k < ns2; k++) {
1772                                                 bmv1 = mesh_vert(vm, i, ns2, k)->v;
1773                                                 if (!bmv1)
1774                                                         bmv1 = mesh_vert(vm, i, 0, k)->v;
1775                                                 if (!(j > 0 && bmv1 == vv[j - 1])) {
1776                                                         BLI_assert(bmv1 != NULL);
1777                                                         BLI_array_append(vv, bmv1);
1778                                                         j++;
1779                                                 }
1780                                         }
1781                                 }
1782                                 bmv1 = mesh_vert(vm, i, ns2, ns2)->v;
1783                                 if (!bmv1)
1784                                         bmv1 = mesh_vert(vm, i, 0, ns2)->v;
1785                                 if (!(j > 0 && bmv1 == vv[j - 1])) {
1786                                         BLI_assert(bmv1 != NULL);
1787                                         BLI_array_append(vv, bmv1);
1788                                         j++;
1789                                 }
1790                                 if (!v->next->ebev) {
1791                                         for (k = ns - ns2; k < ns; k++) {
1792                                                 bmv1 = mesh_vert(vm, i, ns2, k)->v;
1793                                                 if (!bmv1)
1794                                                         bmv1 = mesh_vert(vm, i, 0, k)->v;
1795                                                 if (!(j > 0 && bmv1 == vv[j - 1])) {
1796                                                         BLI_assert(bmv1 != NULL);
1797                                                         BLI_array_append(vv, bmv1);
1798                                                         j++;
1799                                                 }
1800                                         }
1801                                 }
1802                         }
1803                         else {
1804                                 BLI_assert(mesh_vert(vm, i, 0, 0)->v != NULL);
1805                                 BLI_array_append(vv, mesh_vert(vm, i, 0, 0)->v);
1806                                 j++;
1807                         }
1808                 } while ((v = v->next) != vm->boundstart);
1809                 if (vv[0] == vv[j - 1])
1810                         j--;
1811                 bev_create_ngon(bm, vv, j, NULL, f, true);
1812
1813                 BLI_array_free(vv);
1814         }
1815 }
1816 #endif
1817
1818 static VMesh *new_adj_subdiv_vmesh(MemArena *mem_arena, int count, int seg, BoundVert *bounds)
1819 {
1820         VMesh *vm;
1821
1822         vm = (VMesh *)BLI_memarena_alloc(mem_arena, sizeof(VMesh));
1823         vm->count = count;
1824         vm->seg = seg;
1825         vm->boundstart = bounds;
1826         vm->mesh = (NewVert *)BLI_memarena_alloc(mem_arena, count * (1 + seg / 2) * (1 + seg) * sizeof(NewVert));
1827         vm->mesh_kind = M_ADJ_SUBDIV;
1828         return vm;
1829 }
1830
1831 /* VMesh verts for vertex i have data for (i, 0 <= j <= ns2, 0 <= k <= ns), where ns2 = floor(nseg / 2).
1832  * But these overlap data from previous and next i: there are some forced equivalences.
1833  * Let's call these indices the canonical ones: we will just calculate data for these
1834  *    0 <= j <= ns2, 0 <= k < ns2  (for odd ns2)
1835  *    0 <= j < ns2, 0 <= k <= ns2  (for even ns2)
1836  *        also (j=ns2, k=ns2) at i=0 (for even ns2)
1837  * This function returns the canonical one for any i, j, k in [0,n],[0,ns],[0,ns] */
1838 static NewVert *mesh_vert_canon(VMesh *vm, int i, int j, int k)
1839 {
1840         int n, ns, ns2, odd;
1841         NewVert *ans;
1842
1843         n = vm->count;
1844         ns = vm->seg;
1845         ns2 = ns / 2;
1846         odd = ns % 2;
1847         BLI_assert(0 <= i && i <= n && 0 <= j && j <= ns && 0 <= k && k <= ns);
1848
1849         if (!odd && j == ns2 && k == ns2)
1850                 ans = mesh_vert(vm, 0, j, k);
1851         else if (j <= ns2 - 1 + odd && k <= ns2)
1852                 ans = mesh_vert(vm, i, j, k);
1853         else if (k <= ns2)
1854                 ans = mesh_vert(vm, (i + n - 1) % n, k, ns - j);
1855         else
1856                 ans = mesh_vert(vm, (i + 1) % n, ns - k, j);
1857         return ans;
1858 }
1859
1860 static int is_canon(VMesh *vm, int i, int j, int k)
1861 {
1862         int ns2 = vm->seg / 2;
1863         if (vm->seg % 2 == 1)
1864                 return (j <= ns2 && k <= ns2);
1865         else
1866                 return ((j < ns2 && k <= ns2) || (j == ns2 && k == ns2 && i == 0));
1867 }
1868
1869 /* Copy the vertex data to all of vm verts from canonical ones */
1870 static void vmesh_copy_equiv_verts(VMesh *vm)
1871 {
1872         int n, ns, ns2, i, j, k;
1873         NewVert *v0, *v1;
1874
1875         n = vm->count;
1876         ns = vm->seg;
1877         ns2 = ns / 2;
1878         for (i = 0; i < n; i++) {
1879                 for (j = 0; j <= ns2; j++) {
1880                         for (k = 0; k <= ns; k++) {
1881                                 if (is_canon(vm, i, j, k))
1882                                         continue;
1883                                 v1 = mesh_vert(vm, i, j, k);
1884                                 v0 = mesh_vert_canon(vm, i, j, k);
1885                                 copy_v3_v3(v1->co, v0->co);
1886                                 v1->v = v0->v;
1887                         }
1888                 }
1889         }
1890 }
1891
1892 /* Calculate and return in r_cent the centroid of the center poly */
1893 static void vmesh_center(VMesh *vm, float r_cent[3])
1894 {
1895         int n, ns2, i;
1896
1897         n = vm->count;
1898         ns2 = vm->seg / 2;
1899         if (vm->seg % 2) {
1900                 zero_v3(r_cent);
1901                 for (i = 0; i < n; i++) {
1902                         add_v3_v3(r_cent, mesh_vert(vm, i, ns2, ns2)->co);
1903                 }
1904                 mul_v3_fl(r_cent, 1.0f / (float) n);
1905         }
1906         else {
1907                 copy_v3_v3(r_cent, mesh_vert(vm, 0, ns2, ns2)->co);
1908         }
1909 }
1910
1911 static void avg4(float co[3], const NewVert *v0, const NewVert *v1,
1912                               const NewVert *v2, const NewVert *v3)
1913 {
1914         add_v3_v3v3(co, v0->co, v1->co);
1915         add_v3_v3(co, v2->co);
1916         add_v3_v3(co, v3->co);
1917         mul_v3_fl(co, 0.25f);
1918 }
1919
1920 /* gamma needed for smooth Catmull-Clark, Sabin modification */
1921 static float sabin_gamma(int n)
1922 {
1923         double ans, k, k2, k4, k6, x, y;
1924
1925         /* precalculated for common cases of n */
1926         if (n < 3)
1927                 return 0.0f;
1928         else if (n == 3)
1929                 ans = 0.065247584f;
1930         else if (n == 4)
1931                 ans = 0.25f;
1932         else if (n == 5)
1933                 ans = 0.401983447f;
1934         else if (n == 6)
1935                 ans = 0.523423277f;
1936         else {
1937                 k = cos(M_PI / (double)n);
1938                 /* need x, real root of x^3 + (4k^2 - 3)x - 2k = 0.
1939                  * answer calculated via Wolfram Alpha */
1940                 k2 = k * k;
1941                 k4 = k2 * k2;
1942                 k6 = k4 * k2;
1943                 y = pow(1.73205080756888 * sqrt(64.0 * k6 - 144.0 * k4 + 135.0 * k2 - 27.0) + 9.0 * k,
1944                  1.0/3.0);
1945                 x = 0.480749856769136 * y - (0.231120424783545 * (12.0 * k2 - 9.0)) / y;
1946                 ans = (k * x + 2.0 * k2 - 1.0) / (x * x * (k * x + 1.0));
1947         }
1948         return (float)ans;
1949 }
1950
1951 /* Fill frac with fractions of way along ring 0 for vertex i, for use with interp_range function */
1952 static void fill_vmesh_fracs(VMesh *vm, float *frac, int i)
1953 {
1954         int k, ns;
1955         float total = 0.0f;
1956
1957         ns = vm->seg;
1958         frac[0] = 0.0f;
1959         for (k = 0; k < ns; k++) {
1960                 total += len_v3v3(mesh_vert(vm, i, 0, k)->co, mesh_vert(vm, i, 0, k + 1)->co);
1961                 frac[k + 1] = total;
1962         }
1963         if (total > BEVEL_EPSILON) {
1964                 for (k = 1; k <= ns; k++)
1965                         frac[k] /= total;
1966         }
1967 }
1968
1969 /* Like fill_vmesh_fracs but want fractions for profile points of bndv, with ns segments */
1970 static void fill_profile_fracs(BoundVert *bndv, float *frac, int ns)
1971 {
1972         int k;
1973         float co[3], nextco[3];
1974         float total = 0.0f;
1975
1976         frac[0] = 0.0f;
1977         copy_v3_v3(co, bndv->nv.co);
1978         for (k = 0; k < ns; k++) {
1979                 get_bndv_profile_point(bndv, 2.0f * (float) (k + 1) / (float) ns, nextco);
1980                 total += len_v3v3(co, nextco);
1981                 frac[k + 1] = total;
1982                 copy_v3_v3(co, nextco);
1983         }
1984         if (total > BEVEL_EPSILON) {
1985                 for (k = 1; k <= ns; k++)
1986                 frac[k] /= total;
1987         }
1988 }
1989
1990 /* Return i such that frac[i] <= f <= frac[i + 1], where frac[n] == 1.0
1991  * and put fraction of rest of way between frac[i] and frac[i + 1] into r_rest */
1992 static int interp_range(const float *frac, int n, const float f, float *r_rest)
1993 {
1994         int i;
1995         float rest;
1996
1997         /* could binary search in frac, but expect n to be reasonably small */
1998         for (i = 0; i < n; i++) {
1999                 if (f <= frac[i + 1]) {
2000                         rest = f - frac[i];
2001                         if (rest == 0)
2002                                 *r_rest = 0.0f;
2003                         else
2004                                 *r_rest = rest / (frac[i + 1] - frac[i]);
2005                         return i;
2006                 }
2007         }
2008         *r_rest = 0.0f;
2009         return n;
2010 }
2011
2012 /* Interpolate given vmesh to make one with target nseg border vertices on the profiles */
2013 static VMesh *interp_vmesh(MemArena *mem_arena, VMesh *vm0, int nseg)
2014 {
2015         int n, ns0, nseg2, odd, i, j, k, j0, k0;
2016         float *prev_frac, *frac, *new_frac, f, restj, restk;
2017         float quad[4][3], co[3], center[3];
2018         VMesh *vm1;
2019         BoundVert *bndv;
2020
2021         n = vm0->count;
2022         ns0 = vm0->seg;
2023         nseg2 = nseg / 2;
2024         odd = nseg % 2;
2025         vm1 = new_adj_subdiv_vmesh(mem_arena, n, nseg, vm0->boundstart);
2026         prev_frac = (float *)BLI_memarena_alloc(mem_arena, (ns0 + 1) * sizeof(float));
2027         frac = (float *)BLI_memarena_alloc(mem_arena, (ns0 + 1) * sizeof(float));
2028         new_frac = (float *)BLI_memarena_alloc(mem_arena, (nseg + 1) * sizeof(float));
2029
2030         fill_vmesh_fracs(vm0, prev_frac, n - 1);
2031         bndv = vm0->boundstart;
2032         for (i = 0; i < n; i++) {
2033                 fill_vmesh_fracs(vm0, frac, i);
2034                 fill_profile_fracs(bndv, new_frac, nseg);
2035                 for (j = 0; j <= nseg2 -1 + odd; j++) {
2036                         for (k = 0; k <= nseg2; k++) {
2037                                 f = new_frac[k];
2038                                 k0 = interp_range(frac, ns0, f, &restk);
2039                                 f = 1.0f - new_frac[j];
2040                                 j0 = interp_range(prev_frac, ns0, f, &restj);
2041                                 if (restj < BEVEL_EPSILON) {
2042                                         j0 = ns0 - j0;
2043                                         restj = 0.0f;
2044                                 }
2045                                 else {
2046                                         j0 = ns0 - j0 - 1;
2047                                         restj = 1.0f - restj;
2048                                 }
2049                                 /* Use bilinear interpolation within the source quad; could be smarter here */
2050                                 if (restj < BEVEL_EPSILON && restk < BEVEL_EPSILON) {
2051                                         copy_v3_v3(co, mesh_vert_canon(vm0, i, j0, k0)->co);
2052                                 }
2053                                 else {
2054                                         copy_v3_v3(quad[0], mesh_vert_canon(vm0, i, j0, k0)->co);
2055                                         copy_v3_v3(quad[1], mesh_vert_canon(vm0, i, j0, k0 + 1)->co);
2056                                         copy_v3_v3(quad[2], mesh_vert_canon(vm0, i, j0 + 1, k0 + 1)->co);
2057                                         copy_v3_v3(quad[3], mesh_vert_canon(vm0, i, j0 + 1, k0)->co);
2058                                         interp_bilinear_quad_v3(quad, restk, restj, co);
2059                                 }
2060                                 copy_v3_v3(mesh_vert(vm1, i, j, k)->co, co);
2061                         }
2062                 }
2063                 bndv = bndv->next;
2064                 memcpy(prev_frac, frac, (ns0 + 1) * sizeof(float));
2065         }
2066         if (!odd) {
2067                 vmesh_center(vm0, center);
2068                 copy_v3_v3(mesh_vert(vm1, 0, nseg2, nseg2)->co, center);
2069         }
2070         vmesh_copy_equiv_verts(vm1);
2071         return vm1;
2072 }
2073
2074 /* Do one step of cubic subdivision (Catmull-Clark), with special rules at boundaries.
2075  * For now, this is written assuming vm0->nseg is even and > 0.
2076  * We are allowed to modify vm0, as it will not be used after this call.
2077  * See Levin 1999 paper: "Filling an N-sided hole using combined subdivision schemes". */
2078 static VMesh *cubic_subdiv(MemArena *mem_arena, VMesh *vm0)
2079 {
2080         int n, ns0, ns20, ns1;
2081         int i, j, k, inext;
2082         float co[3], co1[3], co2[3], acc[3];
2083         float beta, gamma, u;
2084         VMesh *vm1;
2085         BoundVert *bndv;
2086         
2087         n = vm0->count;
2088         ns0 = vm0->seg;
2089         ns20 = ns0 / 2;
2090         BLI_assert(ns0 % 2 == 0);
2091         ns1 = 2 * ns0;
2092         vm1 = new_adj_subdiv_vmesh(mem_arena, n, ns1, vm0->boundstart);
2093
2094         /* First we adjust the boundary vertices of the input mesh, storing in output mesh */
2095         for (i = 0; i < n; i++) {
2096                 copy_v3_v3(mesh_vert(vm1, i, 0, 0)->co, mesh_vert(vm0, i, 0, 0)->co);
2097                 for (k = 1; k < ns0; k++) {
2098                         /* smooth boundary rule */
2099                         copy_v3_v3(co, mesh_vert(vm0, i, 0, k)->co);
2100                         copy_v3_v3(co1, mesh_vert(vm0, i, 0, k - 1)->co);
2101                         copy_v3_v3(co2, mesh_vert(vm0, i, 0, k + 1)->co);
2102
2103                         add_v3_v3v3(acc, co1, co2);
2104                         madd_v3_v3fl(acc, co, -2.0f);
2105                         madd_v3_v3fl(co, acc, -1.0f/6.0f);
2106                         
2107                         copy_v3_v3(mesh_vert_canon(vm1, i, 0, 2 * k)->co, co);
2108                 }
2109         }
2110         /* now do odd ones in output mesh, based on even ones */
2111         bndv = vm1->boundstart;
2112         for (i = 0; i < n; i++) {
2113                 for (k = 1; k < ns1; k += 2) {
2114                         get_bndv_profile_point(bndv, 2.0f * (float) k / (float) ns1, co);
2115                         copy_v3_v3(co1, mesh_vert_canon(vm1, i, 0, k - 1)->co);
2116                         copy_v3_v3(co2, mesh_vert_canon(vm1, i, 0, k + 1)->co);
2117
2118                         add_v3_v3v3(acc, co1, co2);
2119                         madd_v3_v3fl(acc, co, -2.0f);
2120                         madd_v3_v3fl(co, acc, -1.0f/6.0f);
2121                         
2122                         copy_v3_v3(mesh_vert_canon(vm1, i, 0, k)->co, co);
2123                 }
2124                 bndv = bndv->next;
2125         }
2126         vmesh_copy_equiv_verts(vm1);
2127
2128         /* Copy adjusted verts back into vm0 */
2129         for (i = 0; i < n; i++) {
2130                 for (k = 0; k < ns0; k++) {
2131                         copy_v3_v3(mesh_vert(vm0, i, 0, k)->co,
2132                                    mesh_vert(vm1, i, 0, 2 * k)->co);
2133                 }
2134         }
2135
2136         vmesh_copy_equiv_verts(vm0);
2137
2138         /* Now we do the internal vertices, using standard Catmull-Clark
2139          * and assuming all boundary vertices have valence 4 */
2140         
2141         /* The new face vertices */
2142         for (i = 0; i < n; i++) {
2143                 for (j = 0; j < ns20; j++) {
2144                         for (k = 0; k < ns20; k++) {
2145                                 /* face up and right from (j, k) */
2146                                 avg4(co, mesh_vert(vm0, i, j, k),
2147                                          mesh_vert(vm0, i, j, k + 1),
2148                                          mesh_vert(vm0, i, j + 1, k),
2149                                          mesh_vert(vm0, i, j + 1, k + 1));
2150                                 copy_v3_v3(mesh_vert(vm1, i, 2 * j + 1, 2 * k + 1)->co, co);
2151                         }
2152                 }
2153         }
2154
2155         /* The new vertical edge vertices  */
2156         for (i = 0; i < n; i++) {
2157                 for (j = 0; j < ns20; j++) {
2158                         for (k = 1; k <= ns20; k++) {
2159                                 /* vertical edge between (j, k) and (j+1, k) */
2160                                 avg4(co, mesh_vert(vm0, i, j, k),
2161                                          mesh_vert(vm0, i, j + 1, k),
2162                                          mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k - 1),
2163                                          mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k + 1));
2164                                 copy_v3_v3(mesh_vert(vm1, i, 2 * j + 1, 2 * k)->co, co);
2165                         }
2166                 }
2167         }
2168
2169         /* The new horizontal edge vertices  */
2170         for (i = 0; i < n; i++) {
2171                 for (j = 1; j < ns20; j++) {
2172                         for (k = 0; k < ns20; k++) {
2173                                 /* horizontal edge between (j, k) and (j, k+1) */
2174                                 avg4(co, mesh_vert(vm0, i, j, k),
2175                                          mesh_vert(vm0, i, j, k + 1),
2176                                          mesh_vert_canon(vm1, i, 2 * j - 1, 2 * k + 1),
2177                                          mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k + 1));
2178                                 copy_v3_v3(mesh_vert(vm1, i, 2 * j, 2 * k + 1)->co, co);
2179                         }
2180                 }
2181         }
2182
2183         /* The new vertices, not on border */
2184         gamma = 0.25f;
2185         beta = -gamma;
2186         for (i = 0; i < n; i++) {
2187                 for (j = 1; j < ns20; j++) {
2188                         for (k = 1; k <= ns20; k++) {
2189                                 /* co1 = centroid of adjacent new edge verts */
2190                                 avg4(co1, mesh_vert_canon(vm1, i, 2 * j, 2 * k - 1),
2191                                           mesh_vert_canon(vm1, i, 2 * j, 2 * k + 1),
2192                                           mesh_vert_canon(vm1, i, 2 * j - 1, 2 * k),
2193                                           mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k));
2194                                 /* co2 = centroid of adjacent new face verts */
2195                                 avg4(co2, mesh_vert_canon(vm1, i, 2 * j - 1, 2 * k - 1),
2196                                           mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k - 1),
2197                                           mesh_vert_canon(vm1, i, 2 * j - 1, 2 * k + 1),
2198                                           mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k + 1));
2199                                 /* combine with original vert with alpha, beta, gamma factors */
2200                                 copy_v3_v3(co, co1);  /* alpha = 1.0 */
2201                                 madd_v3_v3fl(co, co2, beta);
2202                                 madd_v3_v3fl(co, mesh_vert(vm0, i, j, k)->co, gamma);
2203                                 copy_v3_v3(mesh_vert(vm1, i, 2 * j, 2 * k)->co, co);
2204                         }
2205                 }
2206         }
2207
2208         vmesh_copy_equiv_verts(vm1);
2209
2210         /* The center vertex is special */
2211         gamma = sabin_gamma(n);
2212         beta = -gamma;
2213         /* accumulate edge verts in co1, face verts in co2 */
2214         zero_v3(co1);
2215         zero_v3(co2);
2216         for (i = 0; i < n; i++) {
2217                 add_v3_v3(co1, mesh_vert(vm1, i, ns0, ns0 - 1)->co);
2218                 add_v3_v3(co2, mesh_vert(vm1, i, ns0 - 1, ns0 - 1)->co);
2219                 add_v3_v3(co2, mesh_vert(vm1, i, ns0 - 1, ns0 + 1)->co);
2220         }
2221         copy_v3_v3(co, co1);
2222         mul_v3_fl(co, 1.0f / (float)n);
2223         madd_v3_v3fl(co, co2, beta / (2.0f *(float)n));
2224         madd_v3_v3fl(co, mesh_vert(vm0, 0, ns20, ns20)->co, gamma);
2225         for (i = 0; i < n; i++)
2226                 copy_v3_v3(mesh_vert(vm1, i, ns0, ns0)->co, co);
2227
2228         /* Final step: sample the boundary vertices at even parameter spacing */
2229         bndv = vm1->boundstart;
2230         for (i = 0; i < n; i++) {
2231                 inext = (i + 1) % n;
2232                 for (k = 0; k <= ns1; k++) {
2233                         u = 2.0f * (float)k / (float)ns1;
2234                         get_bndv_profile_point(bndv, u, co);
2235                         copy_v3_v3(mesh_vert(vm1, i, 0, k)->co, co);
2236                         if (k >= ns0 && k < ns1) {
2237                                 copy_v3_v3(mesh_vert(vm1, inext, ns1 - k, 0)->co, co);
2238                         }
2239                 }
2240                 bndv = bndv->next;
2241         }
2242
2243         return vm1;
2244 }
2245
2246 /* Make a VMesh with nseg segments that covers the unit radius sphere octant
2247  * with center at (0,0,0).
2248  * This has BoundVerts at (1,0,0), (0,1,0) and (0,0,1), with quarter circle arcs
2249  * on the faces for the orthogonal planes through the origin.
2250  */
2251 static VMesh *make_cube_corner_adj_vmesh(MemArena *mem_arena, int nseg)
2252 {
2253         VMesh *vm0, *vm1;
2254         BoundVert *bndv;
2255         int i, j, k, ns2;
2256         float co[3], coa[3], cob[3], coc[3];
2257         float w;
2258
2259         /* initial mesh has 3 sides, 2 segments */
2260         vm0 = new_adj_subdiv_vmesh(mem_arena, 3, 2, NULL);
2261         vm0->count = 0;  // reset, so following loop will end up with correct count
2262         for (i = 0; i < 3; i++) {
2263                 zero_v3(co);
2264                 co[i] = 1.0f;
2265                 add_new_bound_vert(mem_arena, vm0, co);
2266         }
2267         bndv = vm0->boundstart;
2268         for (i = 0; i < 3; i++) {
2269                 /* Get point, 1/2 of the way around profile, on arc between this and next */
2270                 copy_v3_v3(coa, bndv->nv.co);
2271                 copy_v3_v3(cob, bndv->next->nv.co);
2272                 coc[i] = 1.0f;
2273                 coc[(i + 1) % 3] = 1.0f;
2274                 coc[(i + 2) % 3] = 0.0f;
2275                 copy_v3_v3(bndv->profile.midco, coc);
2276                 copy_v3_v3(mesh_vert(vm0, i, 0, 0)->co, coa);
2277                 get_profile_point(&bndv->profile, coa, cob, 1.0f, co);
2278                 copy_v3_v3(mesh_vert(vm0, i, 0, 1)->co, co);
2279                 
2280                 bndv = bndv->next;
2281         }
2282         /* center vertex */
2283         w = 0.57735027f;  /* 1/sqrt(3) */
2284         co[0] = w;
2285         co[1] = w;
2286         co[2] = w;
2287         if (nseg > 2) {
2288                 mul_v3_fl(co, 1.4f);
2289         }
2290         copy_v3_v3(mesh_vert(vm0, 0, 1, 1)->co, co);
2291
2292         vmesh_copy_equiv_verts(vm0);
2293
2294         vm1 = vm0;
2295         while (vm1->seg < nseg) {
2296                 vm1 = cubic_subdiv(mem_arena, vm1);
2297         }
2298         if (vm1->seg != nseg)
2299                 vm1 = interp_vmesh(mem_arena, vm1, nseg);
2300
2301         /* Now snap each vertex to the superellipsoid */
2302         ns2 = nseg / 2;
2303         for (i = 0; i < 3; i++) {
2304                 for (j = 0; j <= ns2; j++) {
2305                         for (k = 0; k <= nseg; k++) {
2306                                 snap_to_sphere(mesh_vert(vm1, i, j, k)->co);
2307                         }
2308                 }
2309         }
2310         return vm1;
2311 }
2312
2313 /* Is this a good candidate for using tri_corner_adj_vmesh? */
2314 static bool tri_corner_test(BevVert *bv)
2315 {
2316         float ang, totang, angdiff;
2317         EdgeHalf *e;
2318         int i;
2319
2320         if (bv->edgecount != 3 || bv->selcount != 3)
2321                 return false;
2322         totang = 0.0f;
2323         for (i = 0; i < 3; i++) {
2324                 e = &bv->edges[i];
2325                 ang = BM_edge_calc_face_angle_signed_ex(e->e, 0.0f);
2326                 if (ang <= (float) M_PI_4 || ang >= 3.0f * (float) M_PI_4)
2327                         return false;
2328                 totang += ang;
2329         }
2330         angdiff = fabsf(totang - 3.0f * (float)M_PI_2);
2331         if (angdiff > (float)M_PI_4)
2332                 return false;
2333         return true;
2334 }
2335
2336 static VMesh *tri_corner_adj_vmesh(BevelParams *bp, BevVert *bv)
2337 {
2338         int i, j, k, ns, ns2;
2339         float co0[3], co1[3], co2[3];
2340         float mat[4][4], v[4];
2341         VMesh *vm;
2342         BoundVert *bndv;
2343
2344         BLI_assert(bv->edgecount == 3 && bv->selcount == 3);
2345         bndv = bv->vmesh->boundstart;
2346         copy_v3_v3(co0, bndv->nv.co);
2347         bndv = bndv->next;
2348         copy_v3_v3(co1, bndv->nv.co);
2349         bndv = bndv->next;
2350         copy_v3_v3(co2, bndv->nv.co);
2351         make_unit_cube_map(co0, co1, co2, bv->v->co, mat);
2352         ns = bp->seg;
2353         ns2 = ns / 2;
2354         vm = make_cube_corner_adj_vmesh(bp->mem_arena, bp->seg);
2355         for (i = 0; i < 3; i++) {
2356                 for (j = 0; j <= ns2; j++) {
2357                         for (k = 0; k <= ns; k++) {
2358                                 copy_v3_v3(v, mesh_vert(vm, i, j, k)->co);
2359                                 v[3] = 1.0f;
2360                                 mul_m4_v4(mat, v);
2361                                 copy_v3_v3(mesh_vert(vm, i, j, k)->co, v);
2362                         }
2363                 }
2364         }
2365
2366         return vm;
2367 }
2368
2369 static VMesh *adj_vmesh(BevelParams *bp, BevVert *bv)
2370 {
2371         int n, ns, i;
2372         VMesh *vm0, *vm1;
2373         float co[3], coa[3], cob[3], dir[3];
2374         BoundVert *bndv;
2375         MemArena *mem_arena = bp->mem_arena;
2376         float fac, fullness;
2377
2378         /* First construct an initial control mesh, with nseg==2 */
2379         n = bv->vmesh->count;
2380         ns = bv->vmesh->seg;
2381         vm0 = new_adj_subdiv_vmesh(mem_arena, n, 2, bv->vmesh->boundstart);
2382
2383         bndv = vm0->boundstart;
2384         zero_v3(co);
2385         for (i = 0; i < n; i++) {
2386                 /* Boundaries just divide input polygon edges into 2 even segments */
2387                 copy_v3_v3(mesh_vert(vm0, i, 0, 0)->co, bndv->nv.co);
2388                 get_bndv_profile_point(bndv, 1.0f, mesh_vert(vm0, i, 0, 1)->co);
2389                 add_v3_v3(co, bndv->nv.co);
2390                 bndv = bndv->next;
2391         }
2392         /* To place center vertex:
2393          * coa is original vertex
2394          * co is centroid of boundary corners
2395          * cob is reflection of coa in across co.
2396          * Calculate 'fullness' = fraction of way
2397          * from co to coa (if positive) or to cob (if negative).
2398          */
2399         copy_v3_v3(coa, bv->v->co);
2400         mul_v3_fl(co, 1.0f / (float)n);
2401         sub_v3_v3v3(cob, co, coa);
2402         add_v3_v3(cob, co);
2403         if (bp->vertex_only)
2404                 fac = 0.25f;
2405         else
2406                 fac = 0.5f;
2407         fullness = 1.0f - fac / 2.0f;
2408         sub_v3_v3v3(dir, coa, co);
2409         if (len_squared_v3(dir) > BEVEL_EPSILON_SQ)
2410                 madd_v3_v3fl(co, dir, fullness);
2411         copy_v3_v3(mesh_vert(vm0, 0, 1, 1)->co, co);
2412         vmesh_copy_equiv_verts(vm0);
2413
2414         vm1 = vm0;
2415         do {
2416                 vm1 = cubic_subdiv(mem_arena, vm1);
2417         } while (vm1->seg < ns);
2418         if (vm1->seg != ns)
2419                 vm1 = interp_vmesh(mem_arena, vm1, ns);
2420         return vm1;
2421 }
2422
2423 static VMesh *pipe_adj_vmesh(BevelParams *bp, BevVert *bv, EdgeHalf *epipe)
2424 {
2425         int i, j, k, n, ns, ns2;
2426         VMesh *vm;
2427         BoundVert *bndv;
2428
2429         vm = adj_vmesh(bp, bv);
2430
2431         /* Now snap all interior coordinates to be on the epipe profile */
2432         n = bv->vmesh->count;
2433         ns = bv->vmesh->seg;
2434         ns2 = ns / 2;
2435         bndv = vm->boundstart;
2436         for (i = 0; i < n; i++) {
2437                 if (bndv->ebev == epipe)
2438                         break;
2439                 bndv = bndv->next;
2440         }
2441         for (i = 0; i < n; i++) {
2442                 for (j = 1; j <= ns2; j++) {
2443                         for (k = 0; k <= ns2; k++) {
2444                                 if (!is_canon(vm, i, j, k))
2445                                         continue;
2446                                 snap_to_profile(bndv, epipe, mesh_vert(vm, i, j, k)->co);
2447                         }
2448                 }
2449         }
2450
2451         return vm;
2452 }
2453
2454 /*
2455  * Given that the boundary is built and the boundary BMVerts have been made,
2456  * calculate the positions of the interior mesh points for the M_ADJ_SUBDIV pattern,
2457  * using cubic subdivision, then make the BMVerts and the new faces. */
2458 static void bevel_build_rings_subdiv(BevelParams *bp, BMesh *bm, BevVert *bv)
2459 {
2460         int n, ns, ns2, odd, i, j, k, ring;
2461         VMesh *vm1, *vm;
2462         BoundVert *v;
2463         BMVert *bmv1, *bmv2, *bmv3, *bmv4;
2464         BMFace *f, *f2, *f23;
2465         EdgeHalf *epipe;
2466
2467         n = bv->vmesh->count;
2468         ns = bv->vmesh->seg;
2469         ns2 = ns / 2;
2470         odd = ns % 2;
2471         BLI_assert(n >= 3 && ns > 1);
2472
2473         epipe = pipe_test(bv);
2474
2475         if (epipe)
2476                 vm1 = pipe_adj_vmesh(bp, bv, epipe);
2477         else if (tri_corner_test(bv))
2478                 vm1 = tri_corner_adj_vmesh(bp, bv);
2479         else
2480                 vm1 = adj_vmesh(bp, bv);
2481
2482         /* copy final vmesh into bv->vmesh, make BMVerts and BMFaces */
2483         vm = bv->vmesh;
2484         for (i = 0; i < n; i ++) {
2485                 for (j = 0; j <= ns2; j++) {
2486                         for (k = 0; k <= ns; k++) {
2487                                 if (j == 0 && (k == 0 || k == ns))
2488                                         continue;  /* boundary corners already made */
2489                                 if (!is_canon(vm, i, j, k))
2490                                         continue;
2491                                 copy_v3_v3(mesh_vert(vm, i, j, k)->co, mesh_vert(vm1, i, j, k)->co);
2492                                 create_mesh_bmvert(bm, vm, i, j, k, bv->v);
2493                         }
2494                 }
2495         }
2496         vmesh_copy_equiv_verts(vm);
2497         /* make the polygons */
2498         v = vm->boundstart;
2499         do {
2500                 i = v->index;
2501                 f = boundvert_rep_face(v);
2502                 f2 = boundvert_rep_face(v->next);
2503                 /* For odd ns, make polys with lower left corner at (i,j,k) for
2504                  *    j in [0, ns2-1], k in [0, ns2].  And then the center ngon.
2505                  * For even ns,
2506                  *    j in [0, ns2-1], k in [0, ns2-1] */
2507                 for (j = 0; j < ns2; j++) {
2508                         for (k = 0; k < ns2 + odd; k++) {
2509                                 bmv1 = mesh_vert(vm, i, j, k)->v;
2510                                 bmv2 = mesh_vert(vm, i, j, k + 1)->v;
2511                                 bmv3 = mesh_vert(vm, i, j + 1, k + 1)->v;
2512                                 bmv4 = mesh_vert(vm, i, j + 1, k)->v;
2513                                 BLI_assert(bmv1 && bmv2 && bmv3 && bmv4);
2514                                 f23 = f;
2515                                 if (odd && k == ns2 && f2 && !v->any_seam)
2516                                         f23 = f2;
2517                                 bev_create_quad_tri_ex(bm, bmv1, bmv2, bmv3, bmv4,
2518                                                        f, f23, f23, f);
2519                         }
2520                 }
2521         } while ((v = v->next) != vm->boundstart);
2522
2523         /* Fix UVs along center lines if even number of segments */
2524         if (!odd) {
2525                 v = vm->boundstart;
2526                 do {
2527                         i = v->index;
2528                         f = boundvert_rep_face(v);
2529                         f2 = boundvert_rep_face(v->next);
2530                         if (!v->any_seam) {
2531                                 for (ring = 1; ring < ns2; ring++) {
2532                                         BMVert *v_uv = mesh_vert(vm, i, ring, ns2)->v;
2533                                         if (v_uv) {
2534                                                 bev_merge_uvs(bm, v_uv);
2535                                         }
2536                                 }
2537                         }
2538                 } while ((v = v->next) != vm->boundstart);
2539                 if (!bv->any_seam)
2540                         bev_merge_uvs(bm, mesh_vert(vm, 0, ns2, ns2)->v);
2541         }
2542
2543         /* center ngon */
2544         if (odd) {
2545                 BMVert **vv = NULL;
2546                 BMFace **vf = NULL;
2547                 BLI_array_staticdeclare(vv, BM_DEFAULT_NGON_STACK_SIZE);
2548                 BLI_array_staticdeclare(vf, BM_DEFAULT_NGON_STACK_SIZE);
2549
2550                 v = vm->boundstart;
2551                 do {
2552                         i = v->index;
2553                         BLI_array_append(vv, mesh_vert(vm, i, ns2, ns2)->v);
2554                         BLI_array_append(vf, v->any_seam ? f : boundvert_rep_face(v));
2555                 } while ((v = v->next) != vm->boundstart);
2556                 f = boundvert_rep_face(vm->boundstart);
2557                 bev_create_ngon(bm, vv, BLI_array_count(vv), vf, f, true);
2558
2559                 BLI_array_free(vv);
2560         }
2561 }
2562
2563 static BMFace *bevel_build_poly(BMesh *bm, BevVert *bv)
2564 {
2565         BMFace *f;
2566         int n, k;
2567         VMesh *vm = bv->vmesh;
2568         BoundVert *v;
2569         BMFace *frep;
2570         BMVert **vv = NULL;
2571         BMFace **vf = NULL;
2572         BLI_array_staticdeclare(vv, BM_DEFAULT_NGON_STACK_SIZE);
2573         BLI_array_staticdeclare(vf, BM_DEFAULT_NGON_STACK_SIZE);
2574
2575         frep = boundvert_rep_face(vm->boundstart);
2576         v = vm->boundstart;
2577         n = 0;
2578         do {
2579                 /* accumulate vertices for vertex ngon */
2580                 /* also accumulate faces in which uv interpolation is to happen for each */
2581                 BLI_array_append(vv, v->nv.v);
2582                 BLI_array_append(vf, bv->any_seam ? frep : boundvert_rep_face(v));
2583                 n++;
2584                 if (v->ebev && v->ebev->seg > 1) {
2585                         for (k = 1; k < v->ebev->seg; k++) {
2586                                 BLI_array_append(vv, mesh_vert(vm, v->index, 0, k)->v);
2587                                 BLI_array_append(vf, bv->any_seam ? frep : boundvert_rep_face(v));
2588                                 n++;
2589                         }
2590                 }
2591         } while ((v = v->next) != vm->boundstart);
2592         if (n > 2) {
2593                 f = bev_create_ngon(bm, vv, n, vf, boundvert_rep_face(v), true);
2594         }
2595         else {
2596                 f = NULL;
2597         }
2598         BLI_array_free(vv);
2599         return f;
2600 }
2601
2602 static void bevel_build_trifan(BMesh *bm, BevVert *bv)
2603 {
2604         BMFace *f;
2605         BLI_assert(next_bev(bv, NULL)->seg == 1 || bv->selcount == 1);
2606
2607         f = bevel_build_poly(bm, bv);
2608
2609         if (f) {
2610                 /* we have a polygon which we know starts at the previous vertex, make it into a fan */
2611                 BMLoop *l_fan = BM_FACE_FIRST_LOOP(f)->prev;
2612                 BMVert *v_fan = l_fan->v;
2613
2614                 while (f->len > 3) {
2615                         BMLoop *l_new;
2616                         BMFace *f_new;
2617                         BLI_assert(v_fan == l_fan->v);
2618                         f_new = BM_face_split(bm, f, l_fan, l_fan->next->next, &l_new, NULL, FALSE);
2619
2620                         if (f_new->len > f->len) {
2621                                 f = f_new;
2622                                 if      (l_new->v       == v_fan) { l_fan = l_new; }
2623                                 else if (l_new->next->v == v_fan) { l_fan = l_new->next; }
2624                                 else if (l_new->prev->v == v_fan) { l_fan = l_new->prev; }
2625                                 else { BLI_assert(0); }
2626                         }
2627                         else {
2628                                 if      (l_fan->v       == v_fan) { /* l_fan = l_fan; */ }
2629                                 else if (l_fan->next->v == v_fan) { l_fan = l_fan->next; }
2630                                 else if (l_fan->prev->v == v_fan) { l_fan = l_fan->prev; }
2631                                 else { BLI_assert(0); }
2632                         }
2633                 }
2634         }
2635 }
2636
2637 static void bevel_build_quadstrip(BMesh *bm, BevVert *bv)
2638 {
2639         BMFace *f;
2640         BLI_assert(bv->selcount == 2);
2641
2642         f = bevel_build_poly(bm, bv);
2643
2644         if (f) {
2645                 /* we have a polygon which we know starts at this vertex, make it into strips */
2646                 EdgeHalf *eh_a = bv->vmesh->boundstart->elast;
2647                 EdgeHalf *eh_b = next_bev(bv, eh_a->next);  /* since (selcount == 2) we know this is valid */
2648                 BMLoop *l_a = BM_face_vert_share_loop(f, eh_a->rightv->nv.v);
2649                 BMLoop *l_b = BM_face_vert_share_loop(f, eh_b->leftv->nv.v);
2650                 int split_count = bv->vmesh->seg + 1;  /* ensure we don't walk past the segments */
2651
2652                 while (f->len > 4 && split_count > 0) {
2653                         BMLoop *l_new;
2654                         BLI_assert(l_a->f == f);
2655                         BLI_assert(l_b->f == f);
2656
2657                         if (l_a-> v == l_b->v || l_a->next == l_b) {
2658                                 /* l_a->v and l_b->v can be the same or such that we'd make a 2-vertex poly */
2659                                 l_a = l_a->prev;
2660                                 l_b = l_b->next;
2661                         }
2662                         else {
2663                                 BM_face_split(bm, f, l_a, l_b, &l_new, NULL, FALSE);
2664                                 f = l_new->f;
2665
2666                                 /* walk around the new face to get the next verts to split */
2667                                 l_a = l_new->prev;
2668                                 l_b = l_new->next->next;
2669                         }
2670                         split_count--;
2671                 }
2672         }
2673 }
2674
2675 /* Given that the boundary is built, now make the actual BMVerts
2676  * for the boundary and the interior of the vertex mesh. */
2677 static void build_vmesh(BevelParams *bp, BMesh *bm, BevVert *bv)
2678 {
2679         MemArena *mem_arena = bp->mem_arena;
2680         VMesh *vm = bv->vmesh;
2681         BoundVert *v, *weld1, *weld2;
2682         int n, ns, ns2, i, k, weld;
2683         float *va, *vb, co[3];
2684
2685         n = vm->count;
2686         ns = vm->seg;
2687         ns2 = ns / 2;
2688
2689         vm->mesh = (NewVert *)BLI_memarena_alloc(mem_arena, n * (ns2 + 1) * (ns + 1) * sizeof(NewVert));
2690
2691         /* special case: two beveled ends welded together */
2692         weld = (bv->selcount == 2) && (vm->count == 2);
2693         weld1 = weld2 = NULL;   /* will hold two BoundVerts involved in weld */
2694
2695         /* make (i, 0, 0) mesh verts for all i */
2696         v = vm->boundstart;
2697         do {
2698                 i = v->index;
2699                 copy_v3_v3(mesh_vert(vm, i, 0, 0)->co, v->nv.co);
2700                 create_mesh_bmvert(bm, vm, i, 0, 0, bv->v);
2701                 v->nv.v = mesh_vert(vm, i, 0, 0)->v;
2702                 if (weld && v->ebev) {
2703                         if (!weld1)
2704                                 weld1 = v;
2705                         else
2706                                 weld2 = v;
2707                 }
2708         } while ((v = v->next) != vm->boundstart);
2709
2710         /* copy other ends to (i, 0, ns) for all i, and fill in profiles for edges */
2711         v = vm->boundstart;
2712         do {
2713                 i = v->index;
2714                 copy_mesh_vert(vm, i, 0, ns, v->next->index, 0, 0);
2715 #ifdef USE_ADJ_SUBDIV
2716                 for (k = 1; k < ns; k++) {
2717                         if (v->ebev && vm->mesh_kind != M_ADJ_SUBDIV) {
2718                                 get_bndv_profile_point(v, 2.0f * (float)k / (float) ns, co);
2719                                 copy_v3_v3(mesh_vert(vm, i, 0, k)->co, co);
2720                                 if (!weld)
2721                                         create_mesh_bmvert(bm, vm, i, 0, k, bv->v);
2722                         }
2723                 }
2724 #else
2725                 if (v->ebev) {
2726                         float midco[3];
2727
2728                         va = mesh_vert(vm, i, 0, 0)->co;
2729                         vb = mesh_vert(vm, i, 0, ns)->co;
2730                         if (bv->edgecount == 3 && bv->selcount == 1) {
2731                                 /* special case: profile cuts the third face, so line it up with that */
2732                                 copy_v3_v3(midco, bv->v->co);
2733                         }
2734                         else {
2735                                 project_to_edge(v->ebev->e, va, vb, midco);
2736                         }
2737                         for (k = 1; k < ns; k++) {
2738                                 get_point_on_round_edge(v->ebev, k, va, midco, vb, co);
2739                                 copy_v3_v3(mesh_vert(vm, i, 0, k)->co, co);
2740                                 if (!weld)
2741                                         create_mesh_bmvert(bm, vm, i, 0, k, bv->v);
2742                         }
2743                 }
2744 #endif
2745         } while ((v = v->next) != vm->boundstart);
2746
2747         if (weld) {
2748                 vm->mesh_kind = M_NONE;
2749                 for (k = 1; k < ns; k++) {
2750                         va = mesh_vert(vm, weld1->index, 0, k)->co;
2751                         vb = mesh_vert(vm, weld2->index, 0, ns - k)->co;
2752                         mid_v3_v3v3(co, va, vb);
2753                         copy_v3_v3(mesh_vert(vm, weld1->index, 0, k)->co, co);
2754                         create_mesh_bmvert(bm, vm, weld1->index, 0, k, bv->v);
2755                 }
2756                 for (k = 1; k < ns; k++)
2757                         copy_mesh_vert(vm, weld2->index, 0, ns - k, weld1->index, 0, k);
2758         }
2759
2760         switch (vm->mesh_kind) {
2761                 case M_NONE:
2762                         /* do nothing */
2763                         break;
2764                 case M_POLY:
2765                         bevel_build_poly(bm, bv);
2766                         break;
2767                 case M_ADJ:
2768 #ifndef USE_ADJ_SUBDIV
2769                         bevel_build_rings(bm, bv);
2770                         break;
2771 #endif
2772                 case M_ADJ_SUBDIV:
2773                         bevel_build_rings_subdiv(bp, bm, bv);
2774                         break;
2775                 case M_TRI_FAN:
2776                         bevel_build_trifan(bm, bv);
2777                         break;
2778                 case M_QUAD_STRIP:
2779                         bevel_build_quadstrip(bm, bv);
2780                         break;
2781         }
2782 }
2783
2784 /* Return the angle between the two faces adjacent to e.
2785  * If there are not two, return 0. */
2786 static float edge_face_angle(EdgeHalf *e)
2787 {
2788         if (e->fprev && e->fnext) {
2789                 /* angle between faces is supplement of angle between face normals */
2790                 return (float)M_PI - angle_normalized_v3v3(e->fprev->no, e->fnext->no);
2791         }
2792         else {
2793                 return 0.0f;
2794         }
2795 }
2796
2797 /* take care, this flag isn't cleared before use, it just so happens that its not set */
2798 #define BM_BEVEL_EDGE_TAG_ENABLE(bme)  BM_ELEM_API_FLAG_ENABLE(  (bme), _FLAG_OVERLAP)
2799 #define BM_BEVEL_EDGE_TAG_DISABLE(bme) BM_ELEM_API_FLAG_DISABLE( (bme), _FLAG_OVERLAP)
2800 #define BM_BEVEL_EDGE_TAG_TEST(bme)    BM_ELEM_API_FLAG_TEST(    (bme), _FLAG_OVERLAP)
2801
2802 /*
2803  * Construction around the vertex
2804  */
2805 static BevVert * bevel_vert_construct(BMesh *bm, BevelParams *bp, BMVert *v)
2806 {
2807         BMEdge *bme;
2808         BevVert *bv;
2809         BMEdge *bme2, *unflagged_bme, *first_bme;
2810         BMFace *f;
2811         BMVert *v1, *v2;
2812         BMIter iter, iter2;
2813         EdgeHalf *e;
2814         float weight, z;
2815         int i, found_shared_face, ccw_test_sum;
2816         int nsel = 0;
2817         int ntot = 0;
2818         int fcnt;
2819
2820         /* Gather input selected edges.
2821          * Only bevel selected edges that have exactly two incident faces.
2822          * Want edges to be ordered so that they share faces.
2823          * There may be one or more chains of shared faces broken by
2824          * gaps where there are no faces.
2825          * TODO: make following work when more than one gap.
2826          */
2827
2828         first_bme = NULL;
2829         BM_ITER_ELEM (bme, &iter, v, BM_EDGES_OF_VERT) {
2830                 fcnt = BM_edge_face_count(bme);
2831                 if (BM_elem_flag_test(bme, BM_ELEM_TAG) && !bp->vertex_only) {
2832                         BLI_assert(fcnt == 2);
2833                         nsel++;
2834                         if (!first_bme)
2835                                 first_bme = bme;
2836                 }
2837                 if (fcnt == 1) {
2838                         /* good to start face chain from this edge */
2839                         first_bme = bme;
2840                 }
2841                 ntot++;
2842
2843                 BM_BEVEL_EDGE_TAG_DISABLE(bme);
2844         }
2845         if (!first_bme)
2846                 first_bme = v->e;
2847
2848         if ((nsel == 0 && !bp->vertex_only) || (ntot < 2 && bp->vertex_only)) {
2849                 /* signal this vert isn't being beveled */
2850                 BM_elem_flag_disable(v, BM_ELEM_TAG);
2851                 return NULL;
2852         }
2853
2854         /* avoid calling BM_vert_edge_count since we loop over edges already */
2855         // ntot = BM_vert_edge_count(v);
2856         // BLI_assert(ntot == BM_vert_edge_count(v));
2857
2858         bv = (BevVert *)BLI_memarena_alloc(bp->mem_arena, (sizeof(BevVert)));
2859         bv->v = v;
2860         bv->edgecount = ntot;
2861         bv->selcount = nsel;
2862         bv->offset = bp->offset;
2863         bv->edges = (EdgeHalf *)BLI_memarena_alloc(bp->mem_arena, ntot * sizeof(EdgeHalf));
2864         bv->vmesh = (VMesh *)BLI_memarena_alloc(bp->mem_arena, sizeof(VMesh));
2865         bv->vmesh->seg = bp->seg;
2866
2867         if (bp->vertex_only) {
2868                 /* if weighted, modify offset by weight */
2869                 if (bp->dvert != NULL && bp->vertex_group != -1) {
2870                         weight = defvert_find_weight(bp->dvert + BM_elem_index_get(v), bp->vertex_group);
2871                         if (weight <= 0.0f) {
2872                                 BM_elem_flag_disable(v, BM_ELEM_TAG);
2873                                 return NULL;
2874                         }
2875                         bv->offset *= weight;
2876                 }
2877         }
2878         BLI_ghash_insert(bp->vert_hash, v, bv);
2879
2880         /* add edges to bv->edges in order that keeps adjacent edges sharing
2881          * a face, if possible */
2882         i = 0;
2883
2884         bme = first_bme;
2885         BM_BEVEL_EDGE_TAG_ENABLE(bme);
2886         e = &bv->edges[0];
2887         e->e = bme;
2888         for (i = 0; i < ntot; i++) {
2889                 if (i > 0) {
2890                         /* find an unflagged edge bme2 that shares a face f with previous bme */
2891                         found_shared_face = 0;
2892                         unflagged_bme = NULL;
2893                         BM_ITER_ELEM (bme2, &iter, v, BM_EDGES_OF_VERT) {
2894                                 if (BM_BEVEL_EDGE_TAG_TEST(bme2))
2895                                         continue;
2896                                 if (!unflagged_bme)
2897                                         unflagged_bme = bme2;
2898                                 if (!bme->l)
2899                                         continue;
2900                                 BM_ITER_ELEM (f, &iter2, bme2, BM_FACES_OF_EDGE) {
2901                                         if (BM_face_edge_share_loop(f, bme)) {
2902                                                 found_shared_face = 1;
2903                                                 break;
2904                                         }
2905                                 }
2906                                 if (found_shared_face)
2907                                         break;
2908                         }
2909                         e = &bv->edges[i];
2910                         if (found_shared_face) {
2911                                 e->e = bme2;
2912                                 e->fprev = f;
2913                                 bv->edges[i - 1].fnext = f;
2914                         }
2915                         else {
2916                                 e->e = unflagged_bme;
2917                         }
2918                 }
2919                 bme = e->e;
2920                 BM_BEVEL_EDGE_TAG_ENABLE(bme);
2921                 if (BM_elem_flag_test(bme, BM_ELEM_TAG) && !bp->vertex_only) {
2922                         e->is_bev = TRUE;
2923                         e->seg = bp->seg;
2924                 }
2925                 else {
2926                         e->is_bev = FALSE;
2927                         e->seg = 0;
2928                 }
2929                 e->is_rev = (bme->v2 == v);
2930         }
2931         /* find wrap-around shared face */
2932         BM_ITER_ELEM (f, &iter2, bme, BM_FACES_OF_EDGE) {
2933                 if (bv->edges[0].e->l && BM_face_edge_share_loop(f, bv->edges[0].e)) {
2934                         if (bv->edges[0].fnext == f)
2935                                 continue;   /* if two shared faces, want the other one now */
2936                         bv->edges[ntot - 1].fnext = f;
2937                         bv->edges[0].fprev = f;
2938                         break;
2939                 }
2940         }
2941
2942         /* if edge array doesn't go CCW around vertex from average normal side,
2943          * reverse the array, being careful to reverse face pointers too */
2944         if (ntot > 1) {
2945                 ccw_test_sum = 0;
2946                 for (i = 0; i < ntot; i++)
2947                         ccw_test_sum += bev_ccw_test(bv->edges[i].e, bv->edges[(i + 1) % ntot].e,
2948                                                      bv->edges[i].fnext);
2949                 if (ccw_test_sum < 0) {
2950                         for (i = 0; i <= (ntot / 2) - 1; i++) {
2951                                 SWAP(EdgeHalf, bv->edges[i], bv->edges[ntot - i - 1]);
2952                                 SWAP(BMFace *, bv->edges[i].fprev, bv->edges[i].fnext);
2953                                 SWAP(BMFace *, bv->edges[ntot - i - 1].fprev, bv->edges[ntot - i - 1].fnext);
2954                         }
2955                         if (ntot % 2 == 1) {
2956                                 i = ntot / 2;
2957                                 SWAP(BMFace *, bv->edges[i].fprev,  bv->edges[i].fnext);
2958                         }
2959                 }
2960         }
2961
2962         for (i = 0, e = bv->edges; i < ntot; i++, e++) {
2963                 e->next = &bv->edges[(i + 1) % ntot];
2964                 e->prev = &bv->edges[(i + ntot - 1) % ntot];
2965
2966                 /* set offsets  */
2967                 if (e->is_bev) {
2968                         /* Convert distance as specified by user into offsets along
2969                          * faces on left side and right side of this edgehalf.
2970                          * Except for percent method, offset will be same on each side. */
2971
2972                         switch (bp->offset_type) {
2973                                 case BEVEL_AMT_OFFSET:
2974                                         e->offset_l_spec = bp->offset;
2975                                         break;
2976                                 case BEVEL_AMT_WIDTH:
2977                                         z = fabsf(2.0f * sinf(edge_face_angle(e) / 2.0f));
2978                                         if (z < BEVEL_EPSILON)
2979                                                 e->offset_l_spec = 0.01f * bp->offset; /* undefined behavior, so tiny bevel */
2980                                         else
2981                                                 e->offset_l_spec = bp->offset / z;
2982                                         break;
2983                                 case BEVEL_AMT_DEPTH:
2984                                         z = fabsf(cosf(edge_face_angle(e) / 2.0f));
2985                                         if (z < BEVEL_EPSILON)
2986                                                 e->offset_l_spec = 0.01f * bp->offset; /* undefined behavior, so tiny bevel */
2987                                         else
2988                                                 e->offset_l_spec = bp->offset / z;
2989                                         break;
2990                                 case BEVEL_AMT_PERCENT:
2991                                         /* offset needs to be such that it meets adjacent edges at percentage of their lengths */
2992                                         v1 = BM_edge_other_vert(e->prev->e, v);
2993