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