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