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