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