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