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