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