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