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