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