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