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