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