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