Cleanup: comment blocks
[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 /* Special case of build_boundary when a single edge is beveled.
1652   * The 'width adjust' part of build_boundary has been done already, and
1653  * efirst is the first beveled edge at vertex bv. */
1654 static void build_boundary_terminal_edge(BevelParams *bp, BevVert *bv, EdgeHalf *efirst, bool construct)
1655 {
1656         MemArena *mem_arena = bp->mem_arena;
1657         VMesh *vm = bv->vmesh;
1658         BoundVert *v;
1659         EdgeHalf *e;
1660         const float *no;
1661         float co[3], d;
1662
1663         e = efirst;
1664         if (bv->edgecount == 2) {
1665                 /* only 2 edges in, so terminate the edge with an artificial vertex on the unbeveled edge */
1666                 no = e->fprev ? e->fprev->no : (e->fnext ? e->fnext->no : NULL);
1667                 offset_in_plane(e, no, true, co);
1668                 if (construct) {
1669                         v = add_new_bound_vert(mem_arena, vm, co);
1670                         v->efirst = v->elast = v->ebev = e;
1671                         e->leftv = v;
1672                 }
1673                 else {
1674                         adjust_bound_vert(e->leftv, co);
1675                 }
1676                 no = e->fnext ? e->fnext->no : (e->fprev ? e->fprev->no : NULL);
1677                 offset_in_plane(e, no, false, co);
1678                 if (construct) {
1679                         v = add_new_bound_vert(mem_arena, vm, co);
1680                         v->efirst = v->elast = e;
1681                         e->rightv = v;
1682                 }
1683                 else {
1684                         adjust_bound_vert(e->rightv, co);
1685                 }
1686                 /* make artifical extra point along unbeveled edge, and form triangle */
1687                 slide_dist(e->next, bv->v, e->offset_l, co);
1688                 if (construct) {
1689                         v = add_new_bound_vert(mem_arena, vm, co);
1690                         v->efirst = v->elast = e->next;
1691                         e->next->leftv = e->next->rightv = v;
1692                         /* could use M_POLY too, but tri-fan looks nicer)*/
1693                         vm->mesh_kind = M_TRI_FAN;
1694                         set_bound_vert_seams(bv);
1695                 }
1696                 else {
1697                         adjust_bound_vert(e->next->leftv, co);
1698                 }
1699         }
1700         else {
1701                 /* More than 2 edges in. Put on-edge verts on all the other edges
1702                  * and join with the beveled edge to make a poly or adj mesh,
1703                  * Because e->prev has offset 0, offset_meet will put co on that edge. */
1704                 /* TODO: should do something else if angle between e and e->prev > 180 */
1705                 offset_meet(e->prev, e, bv->v, e->fprev, false, co);
1706                 if (construct) {
1707                         v = add_new_bound_vert(mem_arena, vm, co);
1708                         v->efirst = e->prev;
1709                         v->elast = v->ebev = e;
1710                         e->leftv = v;
1711                         e->prev->leftv = e->prev->rightv = v;
1712                 }
1713                 else {
1714                         adjust_bound_vert(e->leftv, co);
1715                 }
1716                 e = e->next;
1717                 offset_meet(e->prev, e, bv->v, e->fprev, false, co);
1718                 if (construct) {
1719                         v = add_new_bound_vert(mem_arena, vm, co);
1720                         v->efirst = e->prev;
1721                         v->elast = e;
1722                         e->leftv = e->rightv = v;
1723                         e->prev->rightv = v;
1724                 }
1725                 else {
1726                         adjust_bound_vert(e->leftv, co);
1727                 }
1728                 /* For the edges not adjacent to the beveled edge, slide the bevel amount along. */
1729                 d = efirst->offset_l_spec;
1730                 for (e = e->next; e->next != efirst; e = e->next) {
1731                         slide_dist(e, bv->v, d, co);
1732                         if (construct) {
1733                                 v = add_new_bound_vert(mem_arena, vm, co);
1734                                 v->efirst = v->elast = e;
1735                                 e->leftv = e->rightv = v;
1736                         }
1737                         else {
1738                                 adjust_bound_vert(e->leftv, co);
1739                         }
1740                 }
1741         }
1742         calculate_vm_profiles(bp, bv, vm);
1743
1744         if (bv->edgecount >= 3) {
1745                 /* special case: snap profile to plane of adjacent two edges */
1746                 v = vm->boundstart;
1747                 BLI_assert(v->ebev != NULL);
1748                 move_profile_plane(v, v->efirst, v->next->elast);
1749                 calculate_profile(bp, v);
1750         }
1751
1752         if (construct) {
1753                 set_bound_vert_seams(bv);
1754
1755                 if (vm->count == 2 && bv->edgecount == 3) {
1756                         vm->mesh_kind = M_NONE;
1757                 }
1758                 else if (vm->count == 3) {
1759                         vm->mesh_kind = M_TRI_FAN;
1760                 }
1761                 else {
1762                         vm->mesh_kind = M_POLY;
1763                 }
1764         }
1765 }
1766
1767 /* Return a value that is v if v is within BEVEL_MAX_ADJUST_PCT of the spec (assumed positive),
1768  * else clamp to make it at most that far away from spec */
1769 static float clamp_adjust(float v, float spec)
1770 {
1771         float allowed_delta = spec * (BEVEL_MAX_ADJUST_PCT / 100.0f);
1772
1773         if (v - spec > allowed_delta)
1774                 return spec + allowed_delta;
1775         else if (spec - v > allowed_delta)
1776                 return spec - allowed_delta;
1777         else
1778                 return v;
1779 }
1780
1781 /* Make a circular list of BoundVerts for bv, each of which has the coordinates
1782  * of a vertex on the boundary of the beveled vertex bv->v.
1783  * This may adjust some EdgeHalf widths, and there might have to be
1784  * a subsequent pass to make the widths as consistent as possible.
1785  * The first time through, construct will be true and we are making the BoundVerts
1786  * and setting up the BoundVert and EdgeHalf pointers appropriately.
1787  * For a width consistency path, we just recalculate the coordinates of the
1788  * BoundVerts. If the other ends have been (re)built already, then we
1789  * copy the offsets from there to match, else we use the ideal (user-specified)
1790  * widths.
1791  * Also, if construct, decide on the mesh pattern that will be used inside the boundary.
1792  * Doesn't make the actual BMVerts */
1793 static void build_boundary(BevelParams *bp, BevVert *bv, bool construct)
1794 {
1795         MemArena *mem_arena = bp->mem_arena;
1796         EdgeHalf *efirst, *e, *e2, *e3, *enip, *eip, *eother;
1797         BoundVert *v;
1798         BevVert *bvother;
1799         VMesh *vm;
1800         float co[3];
1801         int nip, nnip;
1802
1803         /* Current bevel does nothing if only one edge into a vertex */
1804         if (bv->edgecount <= 1)
1805                 return;
1806
1807         if (bp->vertex_only) {
1808                 build_boundary_vertex_only(bp, bv, construct);
1809                 return;
1810         }
1811
1812         vm = bv->vmesh;
1813
1814         /* Find a beveled edge to be efirst. Then for each edge, try matching widths to other end. */
1815         e = efirst = next_bev(bv, NULL);
1816         BLI_assert(e->is_bev);
1817         do {
1818                 eother = find_other_end_edge_half(bp, e, &bvother);
1819                 if (eother && bvother->visited && bp->offset_type != BEVEL_AMT_PERCENT) {
1820                         /* try to keep bevel even by matching other end offsets */
1821                         /* sometimes, adjustment can accumulate errors so use the bp->limit_offset to
1822                          * let user limit the adjustment to within a reasonable range around spec */
1823                         if (bp->limit_offset) {
1824                                 e->offset_l = clamp_adjust(eother->offset_r, e->offset_l_spec);
1825                                 e->offset_r = clamp_adjust(eother->offset_l, e->offset_r_spec);
1826                         }
1827                         else {
1828                                 e->offset_l = eother->offset_r;
1829                                 e->offset_r = eother->offset_l;
1830                         }
1831                 }
1832                 else {
1833                         /* reset to user spec */
1834                         e->offset_l = e->offset_l_spec;
1835                         e->offset_r = e->offset_r_spec;
1836                 }
1837         } while ((e = e->next) != efirst);
1838
1839         if (bv->selcount == 1) {
1840                 /* special case: only one beveled edge in */
1841                 build_boundary_terminal_edge(bp, bv, efirst, construct);
1842                 return;
1843         }
1844
1845         /* Here: there is more than one beveled edge.
1846          * We make BoundVerts to connect the sides of the beveled edges.
1847          * Non-beveled edges in between will just join to the appropriate juncture point. */
1848         e = efirst;
1849         do {
1850                 BLI_assert(e->is_bev);
1851                 /* Make the BoundVert for the right side of e; other side will be made
1852                  * when the beveled edge to the left of e is handled.
1853                  * Analyze edges until next beveled edge.
1854                  * They are either "in plane" (preceding and subsequent faces are coplanar)
1855                  * or not. The "non-in-plane" edges effect silhouette and we prefer to slide
1856                  * along one of those if possible. */
1857                 nip = nnip = 0;        /* counts of in-plane / not-in-plane */
1858                 enip = eip = NULL;     /* representatives of each */
1859                 for (e2 = e->next; !e2->is_bev; e2 = e2->next) {
1860                         if (eh_on_plane(e2)) {
1861                                 nip++;
1862                                 eip = e2;
1863                         }
1864                         else {
1865                                 nnip++;
1866                                 enip = e2;
1867                         }
1868                 }
1869                 if (nip == 0 && nnip == 0) {
1870                         offset_meet(e, e2, bv->v, e->fnext, false, co);
1871                 }
1872                 else if (nnip > 0) {
1873                         if (bp->loop_slide && nnip == 1 && good_offset_on_edge_between(e, e2, enip, bv->v)) {
1874                                 offset_on_edge_between(bp, e, e2, enip, bv->v, co);
1875                         }
1876                         else {
1877                                 offset_meet(e, e2, bv->v, NULL, true, co);
1878                         }
1879                 }
1880                 else {
1881                         /* nip > 0 and nnip == 0 */
1882                         if (bp->loop_slide && nip == 1 && good_offset_on_edge_between(e, e2, eip, bv->v)) {
1883                                 offset_on_edge_between(bp, e, e2, eip, bv->v, co);
1884                         }
1885                         else {
1886                                 offset_meet(e, e2, bv->v, e->fnext, true, co);
1887                         }
1888                 }
1889                 if (construct) {
1890                         v = add_new_bound_vert(mem_arena, vm, co);
1891                         v->efirst = e;
1892                         v->elast = e2;
1893                         v->ebev = e2;
1894                         e->rightv = v;
1895                         e2->leftv = v;
1896                         for (e3 = e->next; e3 != e2; e3 = e3->next) {
1897                                 e3->leftv = e3->rightv = v;
1898                         }
1899                 }
1900                 else {
1901                         adjust_bound_vert(e->rightv, co);
1902                 }
1903                 e = e2;
1904         } while (e != efirst);
1905
1906         calculate_vm_profiles(bp, bv, vm);
1907
1908         if (construct) {
1909                 set_bound_vert_seams(bv);
1910
1911                 if (vm->count == 2) {
1912                         vm->mesh_kind = M_NONE;
1913                 }
1914                 else if (efirst->seg == 1) {
1915                         vm->mesh_kind = M_POLY;
1916                 }
1917                 else {
1918                         vm->mesh_kind = M_ADJ;
1919                 }
1920         }
1921 }
1922
1923 /* Do a global pass to try to make offsets as even as possible.
1924  * Consider this graph:
1925  *   nodes = BevVerts
1926  *   edges = { (u,v) } where u and v are nodes such that u and v
1927  *        are connected by a mesh edge that has at least one end
1928  *        whose offset does not match the user spec.
1929  *
1930  * Do a breadth-first search on this graph, starting from nodes
1931  * that have any_adjust=true, and changing all
1932  * not-already-changed offsets on EdgeHalfs to match the
1933  * corresponding ones that changed on the other end.
1934  * The graph is dynamic in the sense that having an offset that
1935  * doesn't meet the user spec can be added as the search proceeds.
1936  * We want this search to be deterministic (not dependent
1937  * on order of processing through hash table), so as to avoid
1938  * flicker to to different decisions made if search is different
1939  * while dragging the offset number in the UI.  So look for the
1940  * lower vertex number when there is a choice of where to start.
1941  *
1942  * Note that this might not process all BevVerts, only the ones
1943  * that need adjustment.
1944  */
1945 static void adjust_offsets(BevelParams *bp)
1946 {
1947         BevVert *bv, *searchbv, *bvother;
1948         int i, searchi;
1949         GHashIterator giter;
1950         EdgeHalf *e, *efirst, *eother;
1951         GSQueue *q;
1952
1953         BLI_assert(!bp->vertex_only);
1954         GHASH_ITER(giter, bp->vert_hash) {
1955                 bv = BLI_ghashIterator_getValue(&giter);
1956                 bv->visited = false;
1957         }
1958
1959         q = BLI_gsqueue_new(sizeof(BevVert *));
1960         /* the following loop terminates because at least one node is visited each time */
1961         for (;;) {
1962                 /* look for root of a connected component in search graph */
1963                 searchbv = NULL;
1964                 searchi = -1;
1965                 GHASH_ITER(giter, bp->vert_hash) {
1966                         bv = BLI_ghashIterator_getValue(&giter);
1967                         if (!bv->visited && any_edge_half_offset_changed(bv)) {
1968                                 i = BM_elem_index_get(bv->v);
1969                                 if (!searchbv || i < searchi) {
1970                                         searchbv = bv;
1971                                         searchi = i;
1972                                 }
1973                         }
1974                 }
1975                 if (searchbv == NULL)
1976                         break;
1977
1978                 BLI_gsqueue_push(q, &searchbv);
1979                 while (!BLI_gsqueue_is_empty(q)) {
1980                         BLI_gsqueue_pop(q, &bv);
1981                         /* If do this check, don't have to check for already-on-queue before push, below */
1982                         if (bv->visited)
1983                                 continue;
1984                         bv->visited = true;
1985                         build_boundary(bp, bv, false);
1986
1987                         e = efirst = &bv->edges[0];
1988                         do {
1989                                 eother = find_other_end_edge_half(bp, e, &bvother);
1990                                 if (eother && !bvother->visited && edge_half_offset_changed(e)) {
1991                                         BLI_gsqueue_push(q, &bvother);
1992                                 }
1993                         } while ((e = e->next) != efirst);
1994                 }
1995         }
1996         BLI_gsqueue_free(q);
1997 }
1998
1999 /* Do the edges at bv form a "pipe"?
2000  * Current definition: 3 or 4 beveled edges, 2 in line with each other,
2001  * with other edges on opposite sides of the pipe if there are 4.
2002  * Also, the vertex boundary should have 3 or 4 vertices in it,
2003  * and all of the faces involved should be parallel to the pipe edges.
2004  * Return the boundary vert whose ebev is one of the pipe edges, and
2005  * whose next boundary vert has a beveled, non-pipe edge. */
2006 static BoundVert *pipe_test(BevVert *bv)
2007 {
2008         EdgeHalf *e, *epipe;
2009         VMesh *vm;
2010         BoundVert *v1, *v2, *v3;
2011         float dir1[3], dir3[3];
2012
2013         vm = bv->vmesh;
2014         if (vm->count < 3 || vm->count > 4 || bv->selcount < 3 || bv->selcount > 4)
2015                 return NULL;
2016
2017         /* find v1, v2, v3 all with beveled edges, where v1 and v3 have collinear edges */
2018         epipe = NULL;
2019         v1 = vm->boundstart;
2020         do {
2021                 v2 = v1->next;
2022                 v3 = v2->next;
2023                 if (v1->ebev && v2->ebev && v3->ebev) {
2024                         sub_v3_v3v3(dir1, bv->v->co, BM_edge_other_vert(v1->ebev->e, bv->v)->co);
2025                         sub_v3_v3v3(dir3, BM_edge_other_vert(v3->ebev->e, bv->v)->co, bv->v->co);
2026                         normalize_v3(dir1);
2027                         normalize_v3(dir3);
2028                         if (angle_normalized_v3v3(dir1, dir3) < BEVEL_EPSILON_ANG) {
2029                                 epipe =  v1->ebev;
2030                                 break;
2031                         }
2032                 }
2033         } while ((v1 = v1->next) != vm->boundstart);
2034
2035         if (!epipe)
2036                 return NULL;
2037
2038         /* check face planes: all should have normals perpendicular to epipe */
2039         for (e = &bv->edges[0]; e != &bv->edges[bv->edgecount]; e++) {
2040                 if (e->fnext) {
2041                         if (dot_v3v3(dir1, e->fnext->no) > BEVEL_EPSILON_BIG)
2042                                 return NULL;
2043                 }
2044         }
2045         return v1;
2046 }
2047
2048 static VMesh *new_adj_vmesh(MemArena *mem_arena, int count, int seg, BoundVert *bounds)
2049 {
2050         VMesh *vm;
2051
2052         vm = (VMesh *)BLI_memarena_alloc(mem_arena, sizeof(VMesh));
2053         vm->count = count;
2054         vm->seg = seg;
2055         vm->boundstart = bounds;
2056         vm->mesh = (NewVert *)BLI_memarena_alloc(mem_arena, count * (1 + seg / 2) * (1 + seg) * sizeof(NewVert));
2057         vm->mesh_kind = M_ADJ;
2058         return vm;
2059 }
2060
2061 /* VMesh verts for vertex i have data for (i, 0 <= j <= ns2, 0 <= k <= ns), where ns2 = floor(nseg / 2).
2062  * But these overlap data from previous and next i: there are some forced equivalences.
2063  * Let's call these indices the canonical ones: we will just calculate data for these
2064  *    0 <= j <= ns2, 0 <= k < ns2  (for odd ns2)
2065  *    0 <= j < ns2, 0 <= k <= ns2  (for even ns2)
2066  *        also (j=ns2, k=ns2) at i=0 (for even ns2)
2067  * This function returns the canonical one for any i, j, k in [0,n],[0,ns],[0,ns] */
2068 static NewVert *mesh_vert_canon(VMesh *vm, int i, int j, int k)
2069 {
2070         int n, ns, ns2, odd;
2071         NewVert *ans;
2072
2073         n = vm->count;
2074         ns = vm->seg;
2075         ns2 = ns / 2;
2076         odd = ns % 2;
2077         BLI_assert(0 <= i && i <= n && 0 <= j && j <= ns && 0 <= k && k <= ns);
2078
2079         if (!odd && j == ns2 && k == ns2)
2080                 ans = mesh_vert(vm, 0, j, k);
2081         else if (j <= ns2 - 1 + odd && k <= ns2)
2082                 ans = mesh_vert(vm, i, j, k);
2083         else if (k <= ns2)
2084                 ans = mesh_vert(vm, (i + n - 1) % n, k, ns - j);
2085         else
2086                 ans = mesh_vert(vm, (i + 1) % n, ns - k, j);
2087         return ans;
2088 }
2089
2090 static bool is_canon(VMesh *vm, int i, int j, int k)
2091 {
2092         int ns2 = vm->seg / 2;
2093         if (vm->seg % 2 == 1)
2094                 return (j <= ns2 && k <= ns2);
2095         else
2096                 return ((j < ns2 && k <= ns2) || (j == ns2 && k == ns2 && i == 0));
2097 }
2098
2099 /* Copy the vertex data to all of vm verts from canonical ones */
2100 static void vmesh_copy_equiv_verts(VMesh *vm)
2101 {
2102         int n, ns, ns2, i, j, k;
2103         NewVert *v0, *v1;
2104
2105         n = vm->count;
2106         ns = vm->seg;
2107         ns2 = ns / 2;
2108         for (i = 0; i < n; i++) {
2109                 for (j = 0; j <= ns2; j++) {
2110                         for (k = 0; k <= ns; k++) {
2111                                 if (is_canon(vm, i, j, k))
2112                                         continue;
2113                                 v1 = mesh_vert(vm, i, j, k);
2114                                 v0 = mesh_vert_canon(vm, i, j, k);
2115                                 copy_v3_v3(v1->co, v0->co);
2116                                 v1->v = v0->v;
2117                         }
2118                 }
2119         }
2120 }
2121
2122 /* Calculate and return in r_cent the centroid of the center poly */
2123 static void vmesh_center(VMesh *vm, float r_cent[3])
2124 {
2125         int n, ns2, i;
2126
2127         n = vm->count;
2128         ns2 = vm->seg / 2;
2129         if (vm->seg % 2) {
2130                 zero_v3(r_cent);
2131                 for (i = 0; i < n; i++) {
2132                         add_v3_v3(r_cent, mesh_vert(vm, i, ns2, ns2)->co);
2133                 }
2134                 mul_v3_fl(r_cent, 1.0f / (float) n);
2135         }
2136         else {
2137                 copy_v3_v3(r_cent, mesh_vert(vm, 0, ns2, ns2)->co);
2138         }
2139 }
2140
2141 static void avg4(
2142         float co[3],
2143         const NewVert *v0, const NewVert *v1,
2144         const NewVert *v2, const NewVert *v3)
2145 {
2146         add_v3_v3v3(co, v0->co, v1->co);
2147         add_v3_v3(co, v2->co);
2148         add_v3_v3(co, v3->co);
2149         mul_v3_fl(co, 0.25f);
2150 }
2151
2152 /* gamma needed for smooth Catmull-Clark, Sabin modification */
2153 static float sabin_gamma(int n)
2154 {
2155         double ans, k, k2, k4, k6, x, y;
2156
2157         /* precalculated for common cases of n */
2158         if (n < 3)
2159                 return 0.0f;
2160         else if (n == 3)
2161                 ans = 0.065247584f;
2162         else if (n == 4)
2163                 ans = 0.25f;
2164         else if (n == 5)
2165                 ans = 0.401983447f;
2166         else if (n == 6)
2167                 ans = 0.523423277f;
2168         else {
2169                 k = cos(M_PI / (double)n);
2170                 /* need x, real root of x^3 + (4k^2 - 3)x - 2k = 0.
2171                  * answer calculated via Wolfram Alpha */
2172                 k2 = k * k;
2173                 k4 = k2 * k2;
2174                 k6 = k4 * k2;
2175                 y = pow(M_SQRT3 * sqrt(64.0 * k6 - 144.0 * k4 + 135.0 * k2 - 27.0) + 9.0 * k,
2176                         1.0 / 3.0);
2177                 x = 0.480749856769136 * y - (0.231120424783545 * (12.0 * k2 - 9.0)) / y;
2178                 ans = (k * x + 2.0 * k2 - 1.0) / (x * x * (k * x + 1.0));
2179         }
2180         return (float)ans;
2181 }
2182
2183 /* Fill frac with fractions of way along ring 0 for vertex i, for use with interp_range function */
2184 static void fill_vmesh_fracs(VMesh *vm, float *frac, int i)
2185 {
2186         int k, ns;
2187         float total = 0.0f;
2188
2189         ns = vm->seg;
2190         frac[0] = 0.0f;
2191         for (k = 0; k < ns; k++) {
2192                 total += len_v3v3(mesh_vert(vm, i, 0, k)->co, mesh_vert(vm, i, 0, k + 1)->co);
2193                 frac[k + 1] = total;
2194         }
2195         if (total > BEVEL_EPSILON) {
2196                 for (k = 1; k <= ns; k++)
2197                         frac[k] /= total;
2198         }
2199 }
2200
2201 /* Like fill_vmesh_fracs but want fractions for profile points of bndv, with ns segments */
2202 static void fill_profile_fracs(BevelParams *bp, BoundVert *bndv, float *frac, int ns)
2203 {
2204         int k;
2205         float co[3], nextco[3];
2206         float total = 0.0f;
2207
2208         frac[0] = 0.0f;
2209         copy_v3_v3(co, bndv->nv.co);
2210         for (k = 0; k < ns; k++) {
2211                 get_profile_point(bp, &bndv->profile, k + 1, ns, nextco);
2212                 total += len_v3v3(co, nextco);
2213                 frac[k + 1] = total;
2214                 copy_v3_v3(co, nextco);
2215         }
2216         if (total > BEVEL_EPSILON) {
2217                 for (k = 1; k <= ns; k++) {
2218                         frac[k] /= total;
2219                 }
2220         }
2221 }
2222
2223 /* Return i such that frac[i] <= f <= frac[i + 1], where frac[n] == 1.0
2224  * and put fraction of rest of way between frac[i] and frac[i + 1] into r_rest */
2225 static int interp_range(const float *frac, int n, const float f, float *r_rest)
2226 {
2227         int i;
2228         float rest;
2229
2230         /* could binary search in frac, but expect n to be reasonably small */
2231         for (i = 0; i < n; i++) {
2232                 if (f <= frac[i + 1]) {
2233                         rest = f - frac[i];
2234                         if (rest == 0)
2235                                 *r_rest = 0.0f;
2236                         else
2237                                 *r_rest = rest / (frac[i + 1] - frac[i]);
2238                         if (i == n - 1 && *r_rest == 1.0f) {
2239                                 i = n;
2240                                 *r_rest = 0.0f;
2241                         }
2242                         return i;
2243                 }
2244         }
2245         *r_rest = 0.0f;
2246         return n;
2247 }
2248
2249 /* Interpolate given vmesh to make one with target nseg border vertices on the profiles */
2250 static VMesh *interp_vmesh(BevelParams *bp, VMesh *vm0, int nseg)
2251 {
2252         int n, ns0, nseg2, odd, i, j, k, j0, k0, k0prev;
2253         float *prev_frac, *frac, *new_frac, *prev_new_frac;
2254         float f, restj, restk, restkprev;
2255         float quad[4][3], co[3], center[3];
2256         VMesh *vm1;
2257         BoundVert *bndv;
2258
2259         n = vm0->count;
2260         ns0 = vm0->seg;
2261         nseg2 = nseg / 2;
2262         odd = nseg % 2;
2263         vm1 = new_adj_vmesh(bp->mem_arena, n, nseg, vm0->boundstart);
2264
2265         prev_frac = BLI_array_alloca(prev_frac, (ns0 + 1));
2266         frac = BLI_array_alloca(frac, (ns0 + 1));
2267         new_frac = BLI_array_alloca(new_frac, (nseg + 1));
2268         prev_new_frac = BLI_array_alloca(prev_new_frac, (nseg + 1));
2269
2270         fill_vmesh_fracs(vm0, prev_frac, n - 1);
2271         bndv = vm0->boundstart;
2272         fill_profile_fracs(bp, bndv->prev, prev_new_frac, nseg);
2273         for (i = 0; i < n; i++) {
2274                 fill_vmesh_fracs(vm0, frac, i);
2275                 fill_profile_fracs(bp, bndv, new_frac, nseg);
2276                 for (j = 0; j <= nseg2 - 1 + odd; j++) {
2277                         for (k = 0; k <= nseg2; k++) {
2278                                 f = new_frac[k];
2279                                 k0 = interp_range(frac, ns0, f, &restk);
2280                                 f = prev_new_frac[nseg - j];
2281                                 k0prev = interp_range(prev_frac, ns0, f, &restkprev);
2282                                 j0 = ns0 - k0prev;
2283                                 restj = -restkprev;
2284                                 if (restj > -BEVEL_EPSILON) {
2285                                         restj = 0.0f;
2286                                 }
2287                                 else {
2288                                         j0 = j0 - 1;
2289                                         restj = 1.0f + restj;
2290                                 }
2291                                 /* Use bilinear interpolation within the source quad; could be smarter here */
2292                                 if (restj < BEVEL_EPSILON && restk < BEVEL_EPSILON) {
2293                                         copy_v3_v3(co, mesh_vert_canon(vm0, i, j0, k0)->co);
2294                                 }
2295                                 else {
2296                                         copy_v3_v3(quad[0], mesh_vert_canon(vm0, i, j0, k0)->co);
2297                                         copy_v3_v3(quad[1], mesh_vert_canon(vm0, i, j0, k0 + 1)->co);
2298                                         copy_v3_v3(quad[2], mesh_vert_canon(vm0, i, j0 + 1, k0 + 1)->co);
2299                                         copy_v3_v3(quad[3], mesh_vert_canon(vm0, i, j0 + 1, k0)->co);
2300                                         interp_bilinear_quad_v3(quad, restk, restj, co);
2301                                 }
2302                                 copy_v3_v3(mesh_vert(vm1, i, j, k)->co, co);
2303                         }
2304                 }
2305                 bndv = bndv->next;
2306                 memcpy(prev_frac, frac, (ns0 + 1) * sizeof(float));
2307                 memcpy(prev_new_frac, new_frac, (nseg + 1) * sizeof(float));
2308         }
2309         if (!odd) {
2310                 vmesh_center(vm0, center);
2311                 copy_v3_v3(mesh_vert(vm1, 0, nseg2, nseg2)->co, center);
2312         }
2313         vmesh_copy_equiv_verts(vm1);
2314         return vm1;
2315 }
2316
2317 /* Do one step of cubic subdivision (Catmull-Clark), with special rules at boundaries.
2318  * For now, this is written assuming vm0->nseg is even and > 0.
2319  * We are allowed to modify vm0, as it will not be used after this call.
2320  * See Levin 1999 paper: "Filling an N-sided hole using combined subdivision schemes". */
2321 static VMesh *cubic_subdiv(BevelParams *bp, VMesh *vm0)
2322 {
2323         int n, ns0, ns20, ns1;
2324         int i, j, k, inext;
2325         float co[3], co1[3], co2[3], acc[3];
2326         float beta, gamma;
2327         VMesh *vm1;
2328         BoundVert *bndv;
2329         
2330         n = vm0->count;
2331         ns0 = vm0->seg;
2332         ns20 = ns0 / 2;
2333         BLI_assert(ns0 % 2 == 0);
2334         ns1 = 2 * ns0;
2335         vm1 = new_adj_vmesh(bp->mem_arena, n, ns1, vm0->boundstart);
2336
2337         /* First we adjust the boundary vertices of the input mesh, storing in output mesh */
2338         for (i = 0; i < n; i++) {
2339                 copy_v3_v3(mesh_vert(vm1, i, 0, 0)->co, mesh_vert(vm0, i, 0, 0)->co);
2340                 for (k = 1; k < ns0; k++) {
2341                         /* smooth boundary rule */
2342                         copy_v3_v3(co, mesh_vert(vm0, i, 0, k)->co);
2343                         copy_v3_v3(co1, mesh_vert(vm0, i, 0, k - 1)->co);
2344                         copy_v3_v3(co2, mesh_vert(vm0, i, 0, k + 1)->co);
2345
2346                         add_v3_v3v3(acc, co1, co2);
2347                         madd_v3_v3fl(acc, co, -2.0f);
2348                         madd_v3_v3fl(co, acc, -1.0f / 6.0f);
2349                         
2350                         copy_v3_v3(mesh_vert_canon(vm1, i, 0, 2 * k)->co, co);
2351                 }
2352         }
2353         /* now do odd ones in output mesh, based on even ones */
2354         bndv = vm1->boundstart;
2355         for (i = 0; i < n; i++) {
2356                 for (k = 1; k < ns1; k += 2) {
2357                         get_profile_point(bp, &bndv->profile, k, ns1, co);
2358                         copy_v3_v3(co1, mesh_vert_canon(vm1, i, 0, k - 1)->co);
2359                         copy_v3_v3(co2, mesh_vert_canon(vm1, i, 0, k + 1)->co);
2360
2361                         add_v3_v3v3(acc, co1, co2);
2362                         madd_v3_v3fl(acc, co, -2.0f);
2363                         madd_v3_v3fl(co, acc, -1.0f / 6.0f);
2364                         
2365                         copy_v3_v3(mesh_vert_canon(vm1, i, 0, k)->co, co);
2366                 }
2367                 bndv = bndv->next;
2368         }
2369         vmesh_copy_equiv_verts(vm1);
2370
2371         /* Copy adjusted verts back into vm0 */
2372         for (i = 0; i < n; i++) {
2373                 for (k = 0; k < ns0; k++) {
2374                         copy_v3_v3(mesh_vert(vm0, i, 0, k)->co,
2375                                    mesh_vert(vm1, i, 0, 2 * k)->co);
2376                 }
2377         }
2378
2379         vmesh_copy_equiv_verts(vm0);
2380
2381         /* Now we do the internal vertices, using standard Catmull-Clark
2382          * and assuming all boundary vertices have valence 4 */
2383         
2384         /* The new face vertices */
2385         for (i = 0; i < n; i++) {
2386                 for (j = 0; j < ns20; j++) {
2387                         for (k = 0; k < ns20; k++) {
2388                                 /* face up and right from (j, k) */
2389                                 avg4(co,
2390                                      mesh_vert(vm0, i, j, k),
2391                                      mesh_vert(vm0, i, j, k + 1),
2392                                      mesh_vert(vm0, i, j + 1, k),
2393                                      mesh_vert(vm0, i, j + 1, k + 1));
2394                                 copy_v3_v3(mesh_vert(vm1, i, 2 * j + 1, 2 * k + 1)->co, co);
2395                         }
2396                 }
2397         }
2398
2399         /* The new vertical edge vertices  */
2400         for (i = 0; i < n; i++) {
2401                 for (j = 0; j < ns20; j++) {
2402                         for (k = 1; k <= ns20; k++) {
2403                                 /* vertical edge between (j, k) and (j+1, k) */
2404                                 avg4(co, mesh_vert(vm0, i, j, k),
2405                                          mesh_vert(vm0, i, j + 1, k),
2406                                          mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k - 1),
2407                                          mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k + 1));
2408                                 copy_v3_v3(mesh_vert(vm1, i, 2 * j + 1, 2 * k)->co, co);
2409                         }
2410                 }
2411         }
2412
2413         /* The new horizontal edge vertices  */
2414         for (i = 0; i < n; i++) {
2415                 for (j = 1; j < ns20; j++) {
2416                         for (k = 0; k < ns20; k++) {
2417                                 /* horizontal edge between (j, k) and (j, k+1) */
2418                                 avg4(co, mesh_vert(vm0, i, j, k),
2419                                          mesh_vert(vm0, i, j, k + 1),
2420                                          mesh_vert_canon(vm1, i, 2 * j - 1, 2 * k + 1),
2421                                          mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k + 1));
2422                                 copy_v3_v3(mesh_vert(vm1, i, 2 * j, 2 * k + 1)->co, co);
2423                         }
2424                 }
2425         }
2426
2427         /* The new vertices, not on border */
2428         gamma = 0.25f;
2429         beta = -gamma;
2430         for (i = 0; i < n; i++) {
2431                 for (j = 1; j < ns20; j++) {
2432                         for (k = 1; k <= ns20; k++) {
2433                                 /* co1 = centroid of adjacent new edge verts */
2434                                 avg4(co1, mesh_vert_canon(vm1, i, 2 * j, 2 * k - 1),
2435                                           mesh_vert_canon(vm1, i, 2 * j, 2 * k + 1),
2436                                           mesh_vert_canon(vm1, i, 2 * j - 1, 2 * k),
2437                                           mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k));
2438                                 /* co2 = centroid of adjacent new face verts */
2439                                 avg4(co2, mesh_vert_canon(vm1, i, 2 * j - 1, 2 * k - 1),
2440                                           mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k - 1),
2441                                           mesh_vert_canon(vm1, i, 2 * j - 1, 2 * k + 1),
2442                                           mesh_vert_canon(vm1, i, 2 * j + 1, 2 * k + 1));
2443                                 /* combine with original vert with alpha, beta, gamma factors */
2444                                 copy_v3_v3(co, co1);  /* alpha = 1.0 */
2445                                 madd_v3_v3fl(co, co2, beta);
2446                                 madd_v3_v3fl(co, mesh_vert(vm0, i, j, k)->co, gamma);
2447                                 copy_v3_v3(mesh_vert(vm1, i, 2 * j, 2 * k)->co, co);
2448                         }
2449                 }
2450         }
2451
2452         vmesh_copy_equiv_verts(vm1);
2453
2454         /* The center vertex is special */
2455         gamma = sabin_gamma(n);
2456         beta = -gamma;
2457         /* accumulate edge verts in co1, face verts in co2 */
2458         zero_v3(co1);
2459         zero_v3(co2);
2460         for (i = 0; i < n; i++) {
2461                 add_v3_v3(co1, mesh_vert(vm1, i, ns0, ns0 - 1)->co);
2462                 add_v3_v3(co2, mesh_vert(vm1, i, ns0 - 1, ns0 - 1)->co);
2463                 add_v3_v3(co2, mesh_vert(vm1, i, ns0 - 1, ns0 + 1)->co);
2464         }
2465         copy_v3_v3(co, co1);
2466         mul_v3_fl(co, 1.0f / (float)n);
2467         madd_v3_v3fl(co, co2, beta / (2.0f * (float)n));
2468         madd_v3_v3fl(co, mesh_vert(vm0, 0, ns20, ns20)->co, gamma);
2469         for (i = 0; i < n; i++)
2470                 copy_v3_v3(mesh_vert(vm1, i, ns0, ns0)->co, co);
2471
2472         /* Final step: sample the boundary vertices at even parameter spacing */
2473         bndv = vm1->boundstart;
2474         for (i = 0; i < n; i++) {
2475                 inext = (i + 1) % n;
2476                 for (k = 0; k <= ns1; k++) {
2477                         get_profile_point(bp, &bndv->profile, k, ns1, co);
2478                         copy_v3_v3(mesh_vert(vm1, i, 0, k)->co, co);
2479                         if (k >= ns0 && k < ns1) {
2480                                 copy_v3_v3(mesh_vert(vm1, inext, ns1 - k, 0)->co, co);
2481                         }
2482                 }
2483                 bndv = bndv->next;
2484         }
2485
2486         return vm1;
2487 }
2488
2489 /* Special case for cube corner, when r is PRO_SQUARE_R,
2490  * meaning straight sides */
2491 static VMesh *make_cube_corner_straight(MemArena *mem_arena, int nseg)
2492 {
2493         VMesh *vm;
2494         float co[3];
2495         int i, j, k, ns2;
2496
2497         ns2 = nseg / 2;
2498         vm = new_adj_vmesh(mem_arena, 3, nseg, NULL);
2499         vm->count = 0;  // reset, so following loop will end up with correct count
2500         for (i = 0; i < 3; i++) {
2501                 zero_v3(co);
2502                 co[i] = 1.0f;
2503                 add_new_bound_vert(mem_arena, vm, co);
2504         }
2505         for (i = 0; i < 3; i++) {
2506                 for (j = 0; j <= ns2; j++) {
2507                         for (k = 0; k <= ns2; k++) {
2508                                 if (!is_canon(vm, i, j, k))
2509                                         continue;
2510                                 co[i] = 1.0f;
2511                                 co[(i + 1) % 3] = (float)k * 2.0f / (float)nseg;
2512                                 co[(i + 2) % 3] = (float)j * 2.0f / (float)nseg;
2513                                 copy_v3_v3(mesh_vert(vm, i, j, k)->co, co);
2514                         }
2515                 }
2516         }
2517         vmesh_copy_equiv_verts(vm);
2518         return vm;
2519 }
2520
2521 /* Make a VMesh with nseg segments that covers the unit radius sphere octant
2522  * with center at (0,0,0).
2523  * This has BoundVerts at (1,0,0), (0,1,0) and (0,0,1), with quarter circle arcs
2524  * on the faces for the orthogonal planes through the origin.
2525  */
2526 static VMesh *make_cube_corner_adj_vmesh(BevelParams *bp)
2527 {
2528         MemArena *mem_arena = bp->mem_arena;
2529         int nseg = bp->seg;
2530         float r = bp->pro_super_r;
2531         VMesh *vm0, *vm1;
2532         BoundVert *bndv;
2533         int i, j, k, ns2;
2534         float co[3], coc[3];
2535
2536         if (r == PRO_SQUARE_R)
2537                 return make_cube_corner_straight(mem_arena, nseg);
2538
2539         /* initial mesh has 3 sides, 2 segments */
2540         vm0 = new_adj_vmesh(mem_arena, 3, 2, NULL);
2541         vm0->count = 0;  // reset, so following loop will end up with correct count
2542         for (i = 0; i < 3; i++) {
2543                 zero_v3(co);
2544                 co[i] = 1.0f;
2545                 add_new_bound_vert(mem_arena, vm0, co);
2546         }
2547         bndv = vm0->boundstart;
2548         for (i = 0; i < 3; i++) {
2549                 /* Get point, 1/2 of the way around profile, on arc between this and next */
2550                 coc[i] = 1.0f;
2551                 coc[(i + 1) % 3] = 1.0f;
2552                 coc[(i + 2) % 3] = 0.0f;
2553                 bndv->profile.super_r = r;
2554                 copy_v3_v3(bndv->profile.coa, bndv->nv.co);
2555                 copy_v3_v3(bndv->profile.cob, bndv->next->nv.co);
2556                 copy_v3_v3(bndv->profile.midco, coc);
2557                 copy_v3_v3(mesh_vert(vm0, i, 0, 0)->co, bndv->profile.coa);
2558                 copy_v3_v3(bndv->profile.plane_co, bndv->profile.coa);
2559                 cross_v3_v3v3(bndv->profile.plane_no, bndv->profile.coa, bndv->profile.cob);
2560                 copy_v3_v3(bndv->profile.proj_dir, bndv->profile.plane_no);
2561                 calculate_profile(bp, bndv);
2562                 get_profile_point(bp, &bndv->profile, 1, 2, mesh_vert(vm0, i, 0, 1)->co);
2563                 
2564                 bndv = bndv->next;
2565         }
2566         /* center vertex */
2567         copy_v3_fl(co, M_SQRT1_3);
2568
2569         if (nseg > 2) {
2570                 if (r > 1.5f)
2571                         mul_v3_fl(co, 1.4f);
2572                 else if (r < 0.75f)
2573                         mul_v3_fl(co, 0.6f);
2574         }
2575         copy_v3_v3(mesh_vert(vm0, 0, 1, 1)->co, co);
2576
2577         vmesh_copy_equiv_verts(vm0);
2578
2579         vm1 = vm0;
2580         while (vm1->seg < nseg) {
2581                 vm1 = cubic_subdiv(bp, vm1);
2582         }
2583         if (vm1->seg != nseg)
2584                 vm1 = interp_vmesh(bp, vm1, nseg);
2585
2586         /* Now snap each vertex to the superellipsoid */
2587         ns2 = nseg / 2;
2588         for (i = 0; i < 3; i++) {
2589                 for (j = 0; j <= ns2; j++) {
2590                         for (k = 0; k <= nseg; k++) {
2591                                 snap_to_superellipsoid(mesh_vert(vm1, i, j, k)->co, r, false);
2592                         }
2593                 }
2594         }
2595         return vm1;
2596 }
2597
2598 /* Is this a good candidate for using tri_corner_adj_vmesh? */
2599 static bool tri_corner_test(BevelParams *bp, BevVert *bv)
2600 {
2601         float ang, totang, angdiff;
2602         EdgeHalf *e;
2603         int i;
2604
2605         if (bv->edgecount != 3 || bv->selcount != 3)
2606                 return false;
2607         totang = 0.0f;
2608         for (i = 0; i < 3; i++) {
2609                 e = &bv->edges[i];
2610                 ang = BM_edge_calc_face_angle_signed_ex(e->e, 0.0f);
2611                 if (ang <= (float) M_PI_4 || ang >= 3.0f * (float) M_PI_4)
2612                         return false;
2613                 totang += ang;
2614         }
2615         angdiff = fabsf(totang - 3.0f * (float)M_PI_2);
2616         if ((bp->pro_super_r == PRO_SQUARE_R && angdiff > (float)M_PI / 16.0f) ||
2617             (angdiff > (float)M_PI_4))
2618         {
2619                 return false;
2620         }
2621         return true;
2622 }
2623
2624 static VMesh *tri_corner_adj_vmesh(BevelParams *bp, BevVert *bv)
2625 {
2626         int i, j, k, ns, ns2;
2627         float co0[3], co1[3], co2[3];
2628         float mat[4][4], v[4];
2629         VMesh *vm;
2630         BoundVert *bndv;
2631
2632         BLI_assert(bv->edgecount == 3 && bv->selcount == 3);
2633         bndv = bv->vmesh->boundstart;
2634         copy_v3_v3(co0, bndv->nv.co);
2635         bndv = bndv->next;
2636         copy_v3_v3(co1, bndv->nv.co);
2637         bndv = bndv->next;
2638         copy_v3_v3(co2, bndv->nv.co);
2639         make_unit_cube_map(co0, co1, co2, bv->v->co, mat);
2640         ns = bp->seg;
2641         ns2 = ns / 2;
2642         vm = make_cube_corner_adj_vmesh(bp);
2643         for (i = 0; i < 3; i++) {
2644                 for (j = 0; j <= ns2; j++) {
2645                         for (k = 0; k <= ns; k++) {
2646                                 copy_v3_v3(v, mesh_vert(vm, i, j, k)->co);
2647                                 v[3] = 1.0f;
2648                                 mul_m4_v4(mat, v);
2649                                 copy_v3_v3(mesh_vert(vm, i, j, k)->co, v);
2650                         }
2651                 }
2652         }
2653
2654         return vm;
2655 }
2656
2657 static VMesh *adj_vmesh(BevelParams *bp, BevVert *bv)
2658 {
2659         int n, ns, i;
2660         VMesh *vm0, *vm1;
2661         float co[3], coa[3], cob[3], dir[3];
2662         BoundVert *bndv;
2663         MemArena *mem_arena = bp->mem_arena;
2664         float r, fac, fullness;
2665
2666         /* First construct an initial control mesh, with nseg==2 */
2667         n = bv->vmesh->count;
2668         ns = bv->vmesh->seg;
2669         vm0 = new_adj_vmesh(mem_arena, n, 2, bv->vmesh->boundstart);
2670
2671         bndv = vm0->boundstart;
2672         zero_v3(co);
2673         for (i = 0; i < n; i++) {
2674                 /* Boundaries just divide input polygon edges into 2 even segments */
2675                 copy_v3_v3(mesh_vert(vm0, i, 0, 0)->co, bndv->nv.co);
2676                 get_profile_point(bp, &bndv->profile, 1, 2, mesh_vert(vm0, i, 0, 1)->co);
2677                 add_v3_v3(co, bndv->nv.co);
2678                 bndv = bndv->next;
2679         }
2680         /* To place center vertex:
2681          * coa is original vertex
2682          * co is centroid of boundary corners
2683          * cob is reflection of coa in across co.
2684          * Calculate 'fullness' = fraction of way
2685          * from co to coa (if positive) or to cob (if negative).
2686          */
2687         copy_v3_v3(coa, bv->v->co);
2688         mul_v3_fl(co, 1.0f / (float)n);
2689         sub_v3_v3v3(cob, co, coa);
2690         add_v3_v3(cob, co);
2691         r = bp->pro_super_r;
2692         if (r == 1.0f)
2693                 fullness = 0.0f;
2694         else if (r > 1.0f) {
2695                 if (bp->vertex_only)
2696                         fac = 0.25f;
2697                 else if (r == PRO_SQUARE_R)
2698                         fac = -2.0;
2699                 else
2700                         fac = 0.5f;
2701                 fullness = 1.0f - fac / r;
2702         }
2703         else {
2704                 fullness = r - 1.0f;
2705         }
2706         sub_v3_v3v3(dir, coa, co);
2707         if (len_squared_v3(dir) > BEVEL_EPSILON_SQ)
2708                 madd_v3_v3fl(co, dir, fullness);
2709         copy_v3_v3(mesh_vert(vm0, 0, 1, 1)->co, co);
2710         vmesh_copy_equiv_verts(vm0);
2711
2712         vm1 = vm0;
2713         do {
2714                 vm1 = cubic_subdiv(bp, vm1);
2715         } while (vm1->seg < ns);
2716         if (vm1->seg != ns)
2717                 vm1 = interp_vmesh(bp, vm1, ns);
2718         return vm1;
2719 }
2720
2721 /* Snap co to the closest point on the profile for vpipe projected onto the plane
2722  * containing co with normal in the direction of edge vpipe->ebev.
2723  * For the square profiles, need to decide whether to snap to just one plane
2724  * or to the midpoint of the profile; do so if midline is true. */
2725 static void snap_to_pipe_profile(BoundVert *vpipe, bool midline, float co[3])
2726 {
2727         float va[3], vb[3], edir[3], va0[3], vb0[3], vmid0[3];
2728         float plane[4], m[4][4], minv[4][4], p[3], snap[3];
2729         Profile *pro = &vpipe->profile;
2730         EdgeHalf *e = vpipe->ebev;
2731
2732         copy_v3_v3(va, pro->coa);
2733         copy_v3_v3(vb, pro->cob);
2734
2735         sub_v3_v3v3(edir, e->e->v1->co, e->e->v2->co);
2736
2737         plane_from_point_normal_v3(plane, co, edir);
2738         closest_to_plane_v3(va0, plane, va);
2739         closest_to_plane_v3(vb0, plane, vb);
2740         closest_to_plane_v3(vmid0, plane, pro->midco);
2741         if (make_unit_square_map(va0, vmid0, vb0, m)) {
2742                 /* Transform co and project it onto superellipse */
2743                 if (!invert_m4_m4(minv, m)) {
2744                         /* shouldn't happen */
2745                         BLI_assert(!"failed inverse during pipe profile snap");
2746                         return;
2747                 }
2748                 mul_v3_m4v3(p, minv, co);
2749                 snap_to_superellipsoid(p, pro->super_r, midline);
2750                 mul_v3_m4v3(snap, m, p);
2751                 copy_v3_v3(co, snap);
2752         }
2753         else {
2754                 /* planar case: just snap to line va0--vb0 */
2755                 closest_to_line_segment_v3(p, co, va0, vb0);
2756                 copy_v3_v3(co, p);
2757         }
2758 }
2759
2760 /* See pipe_test for conditions that make 'pipe'; vpipe is the return value from that.
2761  * We want to make an ADJ mesh but then snap the vertices to the profile in a plane
2762  * perpendicular to the pipes.
2763  * A tricky case is for the 'square' profiles and an even nseg: we want certain vertices
2764  * to snap to the midline on the pipe, not just to one plane or the other. */
2765 static VMesh *pipe_adj_vmesh(BevelParams *bp, BevVert *bv, BoundVert *vpipe)
2766 {
2767         int i, j, k, n, ns, ns2, ipipe1, ipipe2;
2768         VMesh *vm;
2769         bool even, midline;
2770
2771         vm = adj_vmesh(bp, bv);
2772
2773         /* Now snap all interior coordinates to be on the epipe profile */
2774         n = bv->vmesh->count;
2775         ns = bv->vmesh->seg;
2776         ns2 = ns / 2;
2777         even = (ns % 2) == 0;
2778         ipipe1 = vpipe->index;
2779         ipipe2 = vpipe->next->next->index;
2780         for (i = 0; i < n; i++) {
2781                 for (j = 1; j <= ns2; j++) {
2782                         for (k = 0; k <= ns2; k++) {
2783                                 if (!is_canon(vm, i, j, k))
2784                                         continue;
2785                                 midline = even && k == ns2 &&
2786                                           ((i == 0 && j == ns2) || (i == ipipe1 || i == ipipe2));
2787                                 snap_to_pipe_profile(vpipe, midline, mesh_vert(vm, i, j, k)->co);
2788                         }
2789                 }
2790         }
2791
2792         return vm;
2793 }
2794
2795 static void get_incident_edges(BMFace *f, BMVert *v, BMEdge **r_e1, BMEdge **r_e2)
2796 {
2797         BMIter iter;
2798         BMEdge *e;
2799
2800         *r_e1 = NULL;
2801         *r_e2 = NULL;
2802         if (!f)
2803                 return;
2804         BM_ITER_ELEM (e, &iter, f, BM_EDGES_OF_FACE) {
2805                 if (e->v1 == v || e->v2 == v) {
2806                         if (*r_e1 == NULL)
2807                                 *r_e1 = e;
2808                         else if (*r_e2 == NULL)
2809                                 *r_e2 = e;
2810                 }
2811         }
2812 }
2813
2814 static BMEdge *find_closer_edge(float *co, BMEdge *e1, BMEdge *e2)
2815 {
2816         float dsq1, dsq2;
2817
2818         BLI_assert(e1 != NULL && e2 != NULL);
2819         dsq1 = dist_squared_to_line_segment_v3(co, e1->v1->co, e1->v2->co);
2820         dsq2 = dist_squared_to_line_segment_v3(co, e2->v1->co, e2->v2->co);
2821         if (dsq1 < dsq2)
2822                 return e1;
2823         else
2824                 return e2;
2825 }
2826
2827 /* Snap co to the closest edge of face f. Return the edge in *r_snap_e,
2828  * the coordinates of snap point in r_ snap_co,
2829  * and the distance squared to the snap point as function return */
2830 static float snap_face_dist_squared(float *co, BMFace *f, BMEdge **r_snap_e, float *r_snap_co)
2831 {
2832         BMIter iter;
2833         BMEdge *beste = NULL;
2834         float d2, beste_d2;
2835         BMEdge *e;
2836         float closest[3];
2837
2838         beste_d2 = 1e20;
2839         BM_ITER_ELEM(e, &iter, f, BM_EDGES_OF_FACE) {
2840                 closest_to_line_segment_v3(closest, co, e->v1->co, e->v2->co);
2841                 d2 = len_squared_v3v3(closest, co);
2842                 if (d2 < beste_d2) {
2843                         beste_d2 = d2;
2844                         beste = e;
2845                         copy_v3_v3(r_snap_co, closest);
2846                 }
2847         }
2848         *r_snap_e = beste;
2849         return beste_d2;
2850 }
2851
2852 /*
2853  * Given that the boundary is built and the boundary BMVerts have been made,
2854  * calculate the positions of the interior mesh points for the M_ADJ pattern,
2855  * using cubic subdivision, then make the BMVerts and the new faces. */
2856 static void bevel_build_rings(BevelParams *bp, BMesh *bm, BevVert *bv)
2857 {
2858         int n, ns, ns2, odd, i, j, k, ring;
2859         VMesh *vm1, *vm;
2860         BoundVert *v;
2861         BMVert *bmv1, *bmv2, *bmv3, *bmv4;
2862         BMFace *f, *f2;
2863         BMEdge *bme, *bme1, *bme2, *bme3;
2864         EdgeHalf *e;
2865         BoundVert *vpipe;
2866         int mat_nr = bp->mat_nr;
2867
2868         n = bv->vmesh->count;
2869         ns = bv->vmesh->seg;
2870         ns2 = ns / 2;
2871         odd = ns % 2;
2872         BLI_assert(n >= 3 && ns > 1);
2873
2874         vpipe = pipe_test(bv);
2875
2876         if (vpipe)
2877                 vm1 = pipe_adj_vmesh(bp, bv, vpipe);
2878         else if (tri_corner_test(bp, bv))
2879                 vm1 = tri_corner_adj_vmesh(bp, bv);
2880         else
2881                 vm1 = adj_vmesh(bp, bv);
2882
2883         /* copy final vmesh into bv->vmesh, make BMVerts and BMFaces */
2884         vm = bv->vmesh;
2885         for (i = 0; i < n; i++) {
2886                 for (j = 0; j <= ns2; j++) {
2887                         for (k = 0; k <= ns; k++) {
2888                                 if (j == 0 && (k == 0 || k == ns))
2889                                         continue;  /* boundary corners already made */
2890                                 if (!is_canon(vm, i, j, k))
2891                                         continue;
2892                                 copy_v3_v3(mesh_vert(vm, i, j, k)->co, mesh_vert(vm1, i, j, k)->co);
2893                                 create_mesh_bmvert(bm, vm, i, j, k, bv->v);
2894                         }
2895                 }
2896         }
2897         vmesh_copy_equiv_verts(vm);
2898         /* make the polygons */
2899         v = vm->boundstart;
2900         do {
2901                 i = v->index;
2902                 f = boundvert_rep_face(v, NULL);
2903                 f2 = boundvert_rep_face(v->next, NULL);
2904                 if (bp->vertex_only)
2905                         e = v->efirst;
2906                 else
2907                         e = v->ebev;
2908                 BLI_assert(e != NULL);
2909                 bme = e->e;
2910                 /* For odd ns, make polys with lower left corner at (i,j,k) for
2911                  *    j in [0, ns2-1], k in [0, ns2].  And then the center ngon.
2912                  * For even ns,
2913                  *    j in [0, ns2-1], k in [0, ns2-1] */
2914                 for (j = 0; j < ns2; j++) {
2915                         for (k = 0; k < ns2 + odd; k++) {
2916                                 bmv1 = mesh_vert(vm, i, j, k)->v;
2917                                 bmv2 = mesh_vert(vm, i, j, k + 1)->v;
2918                                 bmv3 = mesh_vert(vm, i, j + 1, k + 1)->v;
2919                                 bmv4 = mesh_vert(vm, i, j + 1, k)->v;
2920                                 BLI_assert(bmv1 && bmv2 && bmv3 && bmv4);
2921                                 if (bp->vertex_only) {
2922                                         if (j < k) {
2923                                                 if (k == ns2 && j == ns2 - 1) {
2924                                                         bev_create_quad_ex(bm, bmv1, bmv2, bmv3, bmv4, f2, f2, f2, f2,
2925                                                                            NULL, NULL, v->next->efirst->e, bme, mat_nr);
2926                                                 }
2927                                                 else {
2928                                                         bev_create_quad(bm, bmv1, bmv2, bmv3, bmv4, f2, f2, f2, f2, mat_nr);
2929                                                 }
2930                                         }
2931                                         else if (j > k) {
2932                                                 bev_create_quad(bm, bmv1, bmv2, bmv3, bmv4, f2, f2, f2, f2, mat_nr);
2933                                         }
2934                                         else { /* j == k */
2935                                                 /* only one edge attached to v, since vertex_only */
2936                                                 if (e->is_seam) {
2937                                                         bev_create_quad_ex(bm, bmv1, bmv2, bmv3, bmv4, f2, f2, f2, f2,
2938                                                                            bme, NULL, bme, NULL, mat_nr);
2939                                                 }
2940                                                 else {
2941                                                         bev_create_quad_ex(bm, bmv1, bmv2, bmv3, bmv4, f2, f2, f2, f,
2942                                                                            bme, NULL, bme, NULL, mat_nr);
2943                                                 }
2944                                         }
2945                                 }
2946                                 else { /* edge bevel */
2947                                         if (odd) {
2948                                                 if (k == ns2) {
2949                                                         if (e->is_seam) {
2950                                                                 bev_create_quad_ex(bm, bmv1, bmv2, bmv3, bmv4, f, f, f, f,
2951                                                                                    NULL, bme, bme, NULL, mat_nr);
2952                                                         }
2953                                                         else {
2954                                                                 bev_create_quad(bm, bmv1, bmv2, bmv3, bmv4, f, f2, f2, f, mat_nr);
2955                                                         }
2956                                                 }
2957                                                 else {
2958                                                         bev_create_quad(bm, bmv1, bmv2, bmv3, bmv4, f, f, f, f, mat_nr);
2959                                                 }
2960                                         }
2961                                         else {
2962                                                 bme1 = k == ns2 - 1 ? bme : NULL;
2963                                                 bme3 = j == ns2 - 1 ? v->prev->ebev->e : NULL;
2964                                                 bme2 = bme1 != NULL ? bme1 : bme3;
2965                                                 bev_create_quad_ex(bm, bmv1, bmv2, bmv3, bmv4, f, f, f, f,
2966                                                                    NULL, bme1, bme2, bme3, mat_nr);
2967                                         }
2968                                 }
2969                         }
2970                 }
2971         } while ((v = v->next) != vm->boundstart);
2972
2973         /* Fix UVs along center lines if even number of segments */
2974         if (!odd) {
2975                 v = vm->boundstart;
2976                 do {
2977                         i = v->index;
2978                         if (!v->any_seam) {
2979                                 for (ring = 1; ring < ns2; ring++) {
2980                                         BMVert *v_uv = mesh_vert(vm, i, ring, ns2)->v;
2981                                         if (v_uv) {
2982                                                 bev_merge_uvs(bm, v_uv);
2983                                         }
2984                                 }
2985                         }
2986                 } while ((v = v->next) != vm->boundstart);
2987                 bmv1 = mesh_vert(vm, 0, ns2, ns2)->v;
2988                 if (bp->vertex_only || count_bound_vert_seams(bv) <= 1)
2989                         bev_merge_uvs(bm, bmv1);
2990         }
2991
2992         /* center ngon */
2993         if (odd) {
2994                 BMFace *frep;
2995                 BMEdge *frep_e1, *frep_e2, *frep_e;
2996                 BMVert **vv = NULL;
2997                 BMFace **vf = NULL;
2998                 BMEdge **ve = NULL;
2999                 BLI_array_staticdeclare(vv, BM_DEFAULT_NGON_STACK_SIZE);
3000                 BLI_array_staticdeclare(vf, BM_DEFAULT_NGON_STACK_SIZE);
3001                 BLI_array_staticdeclare(ve, BM_DEFAULT_NGON_STACK_SIZE);
3002
3003                 if (bv->any_seam) {
3004                         frep = boundvert_rep_face(vm->boundstart, NULL);
3005                         get_incident_edges(frep, bv->v, &frep_e1, &frep_e2);
3006                 }
3007                 else {
3008                         frep = NULL;
3009                         frep_e1 = frep_e2 = NULL;
3010                 }
3011                 v = vm->boundstart;
3012                 do {
3013                         i = v->index;
3014                         BLI_array_append(vv, mesh_vert(vm, i, ns2, ns2)->v);
3015                         if (frep) {
3016                                 BLI_array_append(vf, frep);
3017                                 frep_e = find_closer_edge(mesh_vert(vm, i, ns2, ns2)->v->co, frep_e1, frep_e2);
3018                                 BLI_array_append(ve, v == vm->boundstart ? NULL : frep_e);
3019                         }
3020                         else {
3021                                 BLI_array_append(vf, boundvert_rep_face(v, NULL));
3022                                 BLI_array_append(ve, NULL);
3023                         }
3024                 } while ((v = v->next) != vm->boundstart);
3025                 bev_create_ngon(bm, vv, BLI_array_count(vv), vf, frep, ve, mat_nr, true);
3026
3027                 BLI_array_free(vv);
3028                 BLI_array_free(vf);
3029                 BLI_array_free(ve);
3030         }
3031 }
3032
3033 /* If we make a poly out of verts around bv, snapping to rep frep, will uv poly have zero area?
3034  * The uv poly is made by snapping all outside-of-frep vertices to the closest edge in frep.
3035  * Assume that this funciton is called when the only inside-of-frep vertex is vm->boundstart.
3036  * The poly will have zero area if the distance of that first vertex to some edge e is zero, and all
3037  * the other vertices snap to e or snap to an edge at a point that is essentially on e too.  */
3038 static bool is_bad_uv_poly(BevVert *bv, BMFace *frep)
3039 {
3040         BoundVert *v;
3041         BMEdge *snape, *firste;
3042         float co[3];
3043         VMesh *vm = bv->vmesh;
3044         float d2;
3045
3046         v = vm->boundstart;
3047         d2 = snap_face_dist_squared(v->nv.v->co, frep, &firste, co);
3048         if (d2 > BEVEL_EPSILON_BIG_SQ || firste == NULL)
3049                 return false;
3050
3051         for (v = v->next; v != vm->boundstart; v = v->next) {
3052                 snap_face_dist_squared(v->nv.v->co, frep, &snape, co);
3053                 if (snape  != firste) {
3054                         d2 = dist_to_line_v3(co, firste->v1->co, firste->v2->co);
3055                         if (d2 > BEVEL_EPSILON_BIG_SQ)
3056                                 return false;
3057                 }
3058         }
3059         return true;
3060 }
3061
3062 static BMFace *bevel_build_poly(BevelParams *bp, BMesh *bm, BevVert *bv)
3063 {
3064         BMFace *f, *frep, *frep2;
3065         int n, k;
3066         VMesh *vm = bv->vmesh;
3067         BoundVert *v;
3068         BMEdge *frep_e1, *frep_e2, *frep_e;
3069         BMVert **vv = NULL;
3070         BMFace **vf = NULL;
3071         BMEdge **ve = NULL;
3072         BLI_array_staticdeclare(vv, BM_DEFAULT_NGON_STACK_SIZE);