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