Make multisegment bevel profiles even for all parameters.

The method for calculating points on the profile for non-circles
and non-lines meant that the segments making up an edge had
uneven widths.
Use a numeric search technique to find superellipse evaluation
points that lead to equal-sized chords across the profile.
Also calculate the actual profile points sooner, so that they
don't have to be recalculated again and again.
This also sets up for a possible later feature of arbitrary
profile shapes, set by user.

1 files changed, 339 insertions, 67 deletions

diff --git a/source/blender/bmesh/tools/bmesh_bevel.c b/source/blender/bmesh/tools/bmesh_bevel.c
index 98e8c793bbc..1071f2522ca 100644
--- a/source/blender/bmesh/tools/bmesh_bevel.c
+++ b/source/blender/bmesh/tools/bmesh_bevel.c
@@ -98,6 +98,10 @@ typedef struct EdgeHalf {
  * The profile is an arc with control points coa, midco,
  * projected onto a plane (plane_no is normal, plane_co is a point on it)
  * via lines in a given direction (proj_dir).
+ * After the parameters are all set, the actual profile points are calculated
+ * and point in prof_co. We also may need profile points for a higher resolution
+ * number of segments, in order to make the vertex mesh pattern, and that goes
+ * in prof_co_2.
  */
 typedef struct Profile {
 	float super_r;       /* superellipse r parameter */
@@ -107,12 +111,23 @@ typedef struct Profile {
 	float plane_no[3];   /* normal of plane to project to */
 	float plane_co[3];   /* coordinate on plane to project to */
 	float proj_dir[3];   /* direction of projection line */
+	float *prof_co;      /* seg+1 profile coordinates (triples of floats) */
+	float *prof_co_2;    /* like prof_co, but for seg power of 2 >= seg */
 } Profile;
 #define PRO_SQUARE_R 4.0f
 #define PRO_CIRCLE_R 2.0f
 #define PRO_LINE_R 1.0f
 #define PRO_SQUARE_IN_R 0.0f
 
+/* Cache result of expensive calculation of u parameter values to
+ * get even spacing on superellipse for current BevelParams seg
+ * and pro_super_r. */
+typedef struct ProfileSpacing {
+	float *uvals;       /* seg+1 u values */
+	float *uvals_2;     /* seg_2+1 u values, seg_2 = power of 2 >= seg */
+	int seg_2;          /* the seg_2 value */
+} ProfileSpacing;
+
 /* An element in a cyclic boundary of a Vertex Mesh (VMesh) */
 typedef struct BoundVert {
 	struct BoundVert *next, *prev;  /* in CCW order */
@@ -161,6 +176,7 @@ typedef struct BevelParams {
 	 * GHash: (key=(BMVert *), value=(BevVert *)) */
 	GHash    *vert_hash;
 	MemArena *mem_arena;    /* use for all allocs while bevel runs, if we need to free we can switch to mempool */
+	ProfileSpacing pro_spacing; /* parameter values for evenly spaced profiles */
 
 	float offset;           /* blender units to offset each side of a beveled edge */
 	int offset_type;        /* how offset is measured; enum defined in bmesh_operators.h */
@@ -1096,60 +1112,166 @@ static void get_point_on_round_edge(EdgeHalf *e, int k,
 }
 #endif
 
-/* Find the point on given profile at parameter u which goes from 0 to 2 as
- * the profile is moved from pro->coa to pro->cob. */
-static void get_profile_point(const Profile *pro, float u, float r_co[3])
+/* Get the coordinate on the superellipse (exponent r),
+ * at parameter value u.  u goes from u to 2 as the
+ * superellipse moves on the quadrant (0,1) to (1,0). */
+static void superellipse_co(float u, float r, float r_co[2])
 {
-	float co[3], co2[3], p[3], vo[3], angle, r, w;
-	float m[4][4];
-
-	if (u <= 0.0f)
-		copy_v3_v3(co, pro->coa);
-	else if (u >= 2.0f)
-		copy_v3_v3(co, pro->cob);
-	else {
-		r = pro->super_r;
-		if (r == 1.0f || !make_unit_square_map(pro->coa, pro->midco, pro->cob, m)) {
-			interp_v3_v3v3(co, pro->coa, pro->cob, u / 2.0f);
-		}
-		else if (r == PRO_SQUARE_IN_R) {
-			/* square inward concave */
-			zero_v3(p);
-			mul_v3_m4v3(vo, m, p);
-			if (u <= 1.0f)
-				interp_v3_v3v3(co, pro->coa, vo, u);
-			else
-				interp_v3_v3v3(co, vo, pro->cob, u - 1.0f);
+	float t;
+	
+	if (u <= 0.0f) {
+		r_co[0] = 0.0f;
+		r_co[1] = 1.0f;
+ 	}
+	else if (u >= 2.0f) {
+		r_co[0] = 1.0f;
+		r_co[1] = 0.0f;
+	}
+	else if (r == PRO_LINE_R) {
+		t = u / 2.0f;
+		r_co[0] = t;
+		r_co[1] = 1.0f - t;
+		
+	}
+	else if (r == PRO_SQUARE_IN_R) {
+		if (u < 1.0f) {
+			r_co[0] = 0.0f;
+			r_co[1] = 1.0f - u;
 		}
-		else if (r >= PRO_SQUARE_R) {
-			/* square outward convex */
-			if (u <= 1.0f)
-				interp_v3_v3v3(co, pro->coa, pro->midco, u);
-			else
-				interp_v3_v3v3(co, pro->midco, pro->cob, u - 1.0f);
+		else {
+			r_co[0] = u - 1.0f;
+			r_co[1] = 0.0f;
+		}
+	}
+	else if (r == PRO_SQUARE_R) {
+		if (u < 1.0f) {
+			r_co[0] = u;
+			r_co[1] = 1.0f;
 		}
 		else {
-			angle = u * (float)M_PI / 4.0f;  /* angle from y axis */
-			p[0] = sinf(angle);
-			p[1] = cosf(angle);
-			p[2] = 0.0f;
-			if (r != PRO_CIRCLE_R) {
-				w = powf(powf(p[0], r) + powf(p[1], r), -1.0f / r);
-				mul_v2_fl(p, w);
-			}
-			mul_v3_m4v3(co, m, p);
+			r_co[0] = 1.0f;
+			r_co[1] = 2.0f - u;
 		}
+		
 	}
-	/* project co onto final profile plane */
-	if (!is_zero_v3(pro->proj_dir)) {
-		add_v3_v3v3(co2, co, pro->proj_dir);
-		if (!isect_line_plane_v3(r_co, co, co2, pro->plane_co, pro->plane_no)) {
-			/* shouldn't happen */
-			copy_v3_v3(r_co, co);
+	else {
+		t = u * (float)M_PI / 4.0f;  /* angle from y axis */
+		r_co[0] = sinf(t);
+		r_co[1] = cosf(t);
+		if (r != PRO_SQUARE_R) {
+			r_co[0] = pow(r_co[0], 2.0f / r);
+			r_co[1] = pow(r_co[1], 2.0f / r);
 		}
 	}
+}
+
+/* Find the point on given profile at parameter i which goes from 0 to n as
+ * the profile is moved from pro->coa to pro->cob.
+ * We assume that n is either the global seg number or a power of 2 less than
+ * or equal to the power of 2 >= seg.
+ * In the latter case, we subsample the profile for seg_2, which will not necessarily
+ * give equal spaced chords, but is in fact more what is desired by the cubic subdivision
+ * method used to make the vmesh pattern. */
+static void get_profile_point(BevelParams *bp, const Profile *pro, int i, int n, float r_co[3])
+{
+	int d;
+
+	if (bp->seg == 1) {
+		if (i == 0)
+			copy_v3_v3(r_co, pro->coa);
+		else
+			copy_v3_v3(r_co, pro->cob);
+	}
+	
 	else {
-		copy_v3_v3(r_co, co);
+		if (n == bp->seg) {
+			BLI_assert(pro->prof_co != NULL);
+			copy_v3_v3(r_co, pro->prof_co + 3 * i);
+		}
+		else {
+			BLI_assert(is_power_of_2_i(n) && n <= bp->pro_spacing.seg_2);
+			/* set d to spacing in prof_co_2 between subsamples */
+			d = bp->pro_spacing.seg_2 / n;
+			copy_v3_v3(r_co, pro->prof_co_2 + 3 * i * d);
+		}
+	}
+}
+
+/* Calculcate the actual coordinate values for bndv's profile.
+ * This is only needed if bp->seg > 1.
+ * Allocate the space for them if that hasn't been done already.
+ * If bp->seg is not a power of 2, also need to calculate
+ * the coordinate values for the power of 2 >= bp->seg,
+ * because the ADJ_SUBDIV pattern needs power-of-2 boundaries
+ * during construction. */
+static void calculate_profile(BevelParams *bp, BoundVert *bndv)
+{
+	int i, k, ns;
+	float *uvals;
+	float co[3], co2[3], p[3], m[4][4];
+	float *prof_co, *prof_co_k;
+	float r;
+	bool need_2, map_ok;
+	Profile *pro = &bndv->profile;
+
+	if (bp->seg == 1)
+		return;
+
+	need_2 = bp->seg != bp->pro_spacing.seg_2;
+	if (!pro->prof_co) {
+		pro->prof_co = (float *)BLI_memarena_alloc(bp->mem_arena, (bp->seg + 1) * 3 * sizeof(float));
+		if (need_2)
+			pro->prof_co_2 = (float *)BLI_memarena_alloc(bp->mem_arena, (bp->pro_spacing.seg_2 + 1) * 3 *sizeof(float));
+		else
+			pro->prof_co_2 = pro->prof_co;
+	}
+	r = pro->super_r;
+	if (r == PRO_LINE_R)
+		map_ok = false;
+	else
+		map_ok = make_unit_square_map(pro->coa, pro->midco, pro->cob, m);
+	for (i = 0; i < 2; i++) {
+		if (i == 0) {
+			ns = bp->seg;
+			uvals = bp->pro_spacing.uvals;
+			prof_co = pro->prof_co;
+		}
+		else {
+			if (!need_2)
+				break;  /* shares coords with pro->prof_co */
+			ns = bp->pro_spacing.seg_2;
+			uvals = bp->pro_spacing.uvals_2;
+			prof_co = pro->prof_co_2;
+		}
+		BLI_assert((r == PRO_LINE_R || uvals != NULL) && prof_co != NULL);
+		for (k = 0; k <= ns; k++) {
+			if (k == 0)
+				copy_v3_v3(co, pro->coa);
+			else if (k == ns)
+				copy_v3_v3(co, pro->cob);
+			else {
+				if (map_ok) {
+					superellipse_co(uvals[k], r, p);
+					p[2] = 0.0f;
+					mul_v3_m4v3(co, m, p);
+				}
+				else {
+					interp_v3_v3v3(co, pro->coa, pro->cob, (float)k / (float)ns);
+				}
+			}
+			/* project co onto final profile plane */
+			prof_co_k = prof_co + 3 * k;
+			if (!is_zero_v3(pro->proj_dir)) {
+				add_v3_v3v3(co2, co, pro->proj_dir);
+				if (!isect_line_plane_v3(prof_co_k, co, co2, pro->plane_co, pro->plane_no)) {
+					/* shouldn't happen */
+					copy_v3_v3(prof_co_k, co);
+				}
+			}
+			else {
+				copy_v3_v3(prof_co_k, co);
+			}
+		}
 	}
 }
 
@@ -1367,6 +1489,7 @@ static void build_boundary(BevelParams *bp, BevVert *bv, bool construct)
 			adjust_bound_vert(e->next->leftv, co);
 		}
 		set_profile_params(bp, vm->boundstart);
+		calculate_profile(bp, vm->boundstart);
 		return;
 	}
 
@@ -1471,6 +1594,7 @@ static void build_boundary(BevelParams *bp, BevVert *bv, bool construct)
 	v = vm->boundstart;
 	do {
 		set_profile_params(bp, v);
+		calculate_profile(bp, v);
 	} while ((v = v->next) != vm->boundstart);
 
 	if (bv->selcount == 1 && bv->edgecount == 3) {
@@ -1478,6 +1602,7 @@ static void build_boundary(BevelParams *bp, BevVert *bv, bool construct)
 		v = vm->boundstart;
 		BLI_assert(v->ebev != NULL);
 		move_profile_plane(v, v->efirst, v->next->elast);
+		calculate_profile(bp, v);
 	}
 
 	if (construct) {
@@ -2130,7 +2255,7 @@ static void fill_vmesh_fracs(VMesh *vm, float *frac, int i)
 }
 
 /* Like fill_vmesh_fracs but want fractions for profile points of bndv, with ns segments */
-static void fill_profile_fracs(BoundVert *bndv, float *frac, int ns)
+static void fill_profile_fracs(BevelParams *bp, BoundVert *bndv, float *frac, int ns)
 {
 	int k;
 	float co[3], nextco[3];
@@ -2139,7 +2264,7 @@ static void fill_profile_fracs(BoundVert *bndv, float *frac, int ns)
 	frac[0] = 0.0f;
 	copy_v3_v3(co, bndv->nv.co);
 	for (k = 0; k < ns; k++) {
-		get_profile_point(&bndv->profile, 2.0f * (float) (k + 1) / (float) ns, nextco);
+		get_profile_point(bp, &bndv->profile, k + 1, ns, nextco);
 		total += len_v3v3(co, nextco);
 		frac[k + 1] = total;
 		copy_v3_v3(co, nextco);
@@ -2178,7 +2303,7 @@ static int interp_range(const float *frac, int n, const float f, float *r_rest)
 }
 
 /* Interpolate given vmesh to make one with target nseg border vertices on the profiles */
-static VMesh *interp_vmesh(MemArena *mem_arena, VMesh *vm0, int nseg)
+static VMesh *interp_vmesh(BevelParams *bp, VMesh *vm0, int nseg)
 {
 	int n, ns0, nseg2, odd, i, j, k, j0, k0, k0prev;
 	float *prev_frac, *frac, *new_frac, *prev_new_frac;
@@ -2191,7 +2316,7 @@ static VMesh *interp_vmesh(MemArena *mem_arena, VMesh *vm0, int nseg)
 	ns0 = vm0->seg;
 	nseg2 = nseg / 2;
 	odd = nseg % 2;
-	vm1 = new_adj_subdiv_vmesh(mem_arena, n, nseg, vm0->boundstart);
+	vm1 = new_adj_subdiv_vmesh(bp->mem_arena, n, nseg, vm0->boundstart);
 
 	prev_frac = BLI_array_alloca(prev_frac, (ns0 + 1));
 	frac = BLI_array_alloca(frac, (ns0 + 1));
@@ -2200,10 +2325,10 @@ static VMesh *interp_vmesh(MemArena *mem_arena, VMesh *vm0, int nseg)
 
 	fill_vmesh_fracs(vm0, prev_frac, n - 1);
 	bndv = vm0->boundstart;
-	fill_profile_fracs(bndv->prev, prev_new_frac, nseg);
+	fill_profile_fracs(bp, bndv->prev, prev_new_frac, nseg);
 	for (i = 0; i < n; i++) {
 		fill_vmesh_fracs(vm0, frac, i);
-		fill_profile_fracs(bndv, new_frac, nseg);
+		fill_profile_fracs(bp, bndv, new_frac, nseg);
 		for (j = 0; j <= nseg2 - 1 + odd; j++) {
 			for (k = 0; k <= nseg2; k++) {
 				f = new_frac[k];
@@ -2249,12 +2374,12 @@ static VMesh *interp_vmesh(MemArena *mem_arena, VMesh *vm0, int nseg)
  * For now, this is written assuming vm0->nseg is even and > 0.
  * We are allowed to modify vm0, as it will not be used after this call.
  * See Levin 1999 paper: "Filling an N-sided hole using combined subdivision schemes". */
-static VMesh *cubic_subdiv(MemArena *mem_arena, VMesh *vm0)
+static VMesh *cubic_subdiv(BevelParams *bp, VMesh *vm0)
 {
 	int n, ns0, ns20, ns1;
 	int i, j, k, inext;
 	float co[3], co1[3], co2[3], acc[3];
-	float beta, gamma, u;
+	float beta, gamma;
 	VMesh *vm1;
 	BoundVert *bndv;
 	
@@ -2263,7 +2388,7 @@ static VMesh *cubic_subdiv(MemArena *mem_arena, VMesh *vm0)
 	ns20 = ns0 / 2;
 	BLI_assert(ns0 % 2 == 0);
 	ns1 = 2 * ns0;
-	vm1 = new_adj_subdiv_vmesh(mem_arena, n, ns1, vm0->boundstart);
+	vm1 = new_adj_subdiv_vmesh(bp->mem_arena, n, ns1, vm0->boundstart);
 
 	/* First we adjust the boundary vertices of the input mesh, storing in output mesh */
 	for (i = 0; i < n; i++) {
@@ -2285,7 +2410,7 @@ static VMesh *cubic_subdiv(MemArena *mem_arena, VMesh *vm0)
 	bndv = vm1->boundstart;
 	for (i = 0; i < n; i++) {
 		for (k = 1; k < ns1; k += 2) {
-			get_profile_point(&bndv->profile, 2.0f * (float) k / (float) ns1, co);
+			get_profile_point(bp, &bndv->profile, k, ns1, co);
 			copy_v3_v3(co1, mesh_vert_canon(vm1, i, 0, k - 1)->co);
 			copy_v3_v3(co2, mesh_vert_canon(vm1, i, 0, k + 1)->co);
 
@@ -2405,8 +2530,7 @@ static VMesh *cubic_subdiv(MemArena *mem_arena, VMesh *vm0)
 	for (i = 0; i < n; i++) {
 		inext = (i + 1) % n;
 		for (k = 0; k <= ns1; k++) {
-			u = 2.0f * (float)k / (float)ns1;
-			get_profile_point(&bndv->profile, u, co);
+			get_profile_point(bp, &bndv->profile, k, ns1, co);
 			copy_v3_v3(mesh_vert(vm1, i, 0, k)->co, co);
 			if (k >= ns0 && k < ns1) {
 				copy_v3_v3(mesh_vert(vm1, inext, ns1 - k, 0)->co, co);
@@ -2455,8 +2579,11 @@ static VMesh *make_cube_corner_straight(MemArena *mem_arena, int nseg)
  * This has BoundVerts at (1,0,0), (0,1,0) and (0,0,1), with quarter circle arcs
  * on the faces for the orthogonal planes through the origin.
  */
-static VMesh *make_cube_corner_adj_vmesh(MemArena *mem_arena, int nseg, float r)
+static VMesh *make_cube_corner_adj_vmesh(BevelParams *bp)
 {
+	MemArena *mem_arena = bp->mem_arena;
+	int nseg = bp->seg;
+	float r = bp->pro_super_r;
 	VMesh *vm0, *vm1;
 	BoundVert *bndv;
 	int i, j, k, ns2;
@@ -2488,7 +2615,8 @@ static VMesh *make_cube_corner_adj_vmesh(MemArena *mem_arena, int nseg, float r)
 		copy_v3_v3(bndv->profile.plane_co, bndv->profile.coa);
 		cross_v3_v3v3(bndv->profile.plane_no, bndv->profile.coa, bndv->profile.cob);
 		copy_v3_v3(bndv->profile.proj_dir, bndv->profile.plane_no);
-		get_profile_point(&bndv->profile, 1.0f, mesh_vert(vm0, i, 0, 1)->co);
+		calculate_profile(bp, bndv);
+		get_profile_point(bp, &bndv->profile, 1, 2, mesh_vert(vm0, i, 0, 1)->co);
 		
 		bndv = bndv->next;
 	}
@@ -2509,10 +2637,10 @@ static VMesh *make_cube_corner_adj_vmesh(MemArena *mem_arena, int nseg, float r)
 
 	vm1 = vm0;
 	while (vm1->seg < nseg) {
-		vm1 = cubic_subdiv(mem_arena, vm1);
+		vm1 = cubic_subdiv(bp, vm1);
 	}
 	if (vm1->seg != nseg)
-		vm1 = interp_vmesh(mem_arena, vm1, nseg);
+		vm1 = interp_vmesh(bp, vm1, nseg);
 
 	/* Now snap each vertex to the superellipsoid */
 	ns2 = nseg / 2;
@@ -2570,7 +2698,7 @@ static VMesh *tri_corner_adj_vmesh(BevelParams *bp, BevVert *bv)
 	make_unit_cube_map(co0, co1, co2, bv->v->co, mat);
 	ns = bp->seg;
 	ns2 = ns / 2;
-	vm = make_cube_corner_adj_vmesh(bp->mem_arena, bp->seg, bp->pro_super_r);
+	vm = make_cube_corner_adj_vmesh(bp);
 	for (i = 0; i < 3; i++) {
 		for (j = 0; j <= ns2; j++) {
 			for (k = 0; k <= ns; k++) {
@@ -2604,7 +2732,7 @@ static VMesh *adj_vmesh(BevelParams *bp, BevVert *bv)
 	for (i = 0; i < n; i++) {
 		/* Boundaries just divide input polygon edges into 2 even segments */
 		copy_v3_v3(mesh_vert(vm0, i, 0, 0)->co, bndv->nv.co);
-		get_profile_point(&bndv->profile, 1.0f, mesh_vert(vm0, i, 0, 1)->co);
+		get_profile_point(bp, &bndv->profile, 1, 2, mesh_vert(vm0, i, 0, 1)->co);
 		add_v3_v3(co, bndv->nv.co);
 		bndv = bndv->next;
 	}
@@ -2642,10 +2770,10 @@ static VMesh *adj_vmesh(BevelParams *bp, BevVert *bv)
 
 	vm1 = vm0;
 	do {
-		vm1 = cubic_subdiv(mem_arena, vm1);
+		vm1 = cubic_subdiv(bp, vm1);
 	} while (vm1->seg < ns);
 	if (vm1->seg != ns)
-		vm1 = interp_vmesh(mem_arena, vm1, ns);
+		vm1 = interp_vmesh(bp, vm1, ns);
 	return vm1;
 }
 
@@ -2934,6 +3062,8 @@ static void build_vmesh(BevelParams *bp, BMesh *bm, BevVert *bv)
 			else {
 				weld2 = v;
 				move_weld_profile_planes(bv, weld1, weld2);
+				calculate_profile(bp, weld1);
+				calculate_profile(bp, weld2);
 			}
 		}
 	} while ((v = v->next) != vm->boundstart);
@@ -2946,7 +3076,7 @@ static void build_vmesh(BevelParams *bp, BMesh *bm, BevVert *bv)
 #ifdef USE_ADJ_SUBDIV
 		for (k = 1; k < ns; k++) {
 			if (v->ebev && vm->mesh_kind != M_ADJ_SUBDIV) {
-				get_profile_point(&v->profile, 2.0f * (float)k / (float) ns, co);
+				get_profile_point(bp, &v->profile, k, ns, co);
 				copy_v3_v3(mesh_vert(vm, i, 0, k)->co, co);
 				if (!weld)
 					create_mesh_bmvert(bm, vm, i, 0, k, bv->v);
@@ -2988,7 +3118,7 @@ static void build_vmesh(BevelParams *bp, BMesh *bm, BevVert *bv)
 				copy_v3_v3(co, vb);
 			}
 			else if (weld2->profile.super_r == PRO_LINE_R &&
-			         weld1->profile.super_r != PRO_LINE_R)
+				 weld1->profile.super_r != PRO_LINE_R)
 			{
 				copy_v3_v3(co, va);
 			}
@@ -3476,6 +3606,147 @@ static void bevel_build_edge_polygons(BMesh *bm, BevelParams *bp, BMEdge *bme)
 	BM_elem_attrs_copy(bm, bm, bme, bme2);
 }
 
+/* Returns the square of the length of the chord from parameter u0 to parameter u1
+ * of superellipse_co. */
+static float superellipse_chord_length_squared(float u0, float u1, float r)
+{
+	float a[2], b[2];
+
+	BLI_assert(u0 >= 0.0f && u1 >= u0 && u1 <= 2.0f);
+	superellipse_co(u0, r, a);
+	superellipse_co(u1, r, b);
+	return len_squared_v2v2(a, b);
+}
+
+/* Find parameter u >= u0 to make chord of squared length d2goal,
+ * from u0 to u on superellipse with parameter r.
+ * If it cannot be found, return -1.0f. */
+static float find_superellipse_chord_u(float u0, float d2goal, float r)
+{
+	float ulow, uhigh, u, d2, d2max;
+	const float dtol = 1e-4f;
+	const float utol = 1e-6f;
+	const float umax = 2.0f;
+
+	if (d2goal == 0.0f)
+		return u0;
+	d2max = superellipse_chord_length_squared(u0, umax, r);
+	if (fabsf(d2goal - d2max) <= dtol)
+		return umax;
+	if (d2goal - d2max > dtol)
+		return -1.0f;
+
+	/* binary search for good u value */
+	ulow = u0;
+	uhigh = umax;
+	do {
+		u = 0.5f * (ulow + uhigh);
+		d2 = superellipse_chord_length_squared(u0, u, r);
+		if (fabsf(d2goal - d2) <= dtol)
+			break;
+		if (d2 < d2goal)
+			ulow = u;
+		else
+			uhigh = u;
+	} while (fabsf(uhigh - ulow) > utol);
+	return u;
+}
+
+/* Find parameters u0, u1, ..., un that divide the quarter-arc superellipse
+ * with parameter r into n even chords.
+ * There is no closed form way of doing this except for a few special
+ * values of r, so this uses binary search to find a chord length that works.
+ * Return the u's in *r_params, which should point to an array of size n+1. */
+static void find_even_superellipse_params(int n, float r, float *r_params)
+{
+	float d2low, d2high, d2, d2final, u;
+	int i, j, n2;
+	const int maxiters = 40;
+	const float d2tol = 1e-6f;
+	const float umax = 2.0f;
+
+	if (r == PRO_CIRCLE_R || r == PRO_LINE_R ||
+	    ((n % 2 == 0) && (r == PRO_SQUARE_IN_R || r == PRO_SQUARE_R)))
+	{
+		/* even parameter spacing works for these cases */
+		for (i = 0; i <= n; i++)
+			r_params[i] = i * 2.0f / (float) n;
+		return;
+	}
+	if (r == PRO_SQUARE_IN_R || r == PRO_SQUARE_R) {
+		/* n is odd, so get one corner-cut chord.
+		 * Solve u == sqrt(2*(1-n2*u)^2) where n2 = floor(n/2) */
+		n2 = n / 2;
+		u = (2.0f * n2 - M_SQRT2) / (2.0f * n2 * n2 - 1);
+		for (i = 0; i < n; i++)
+			r_params[i] = i * u;
+		r_params[n] = umax;
+	}
+	d2low = 2.0f / (n * n);  /* (sqrt(2)/n)**2 */
+	d2high = 2 * d2low;      /* (2/n)**2 */
+	for (i = 0; i < maxiters && fabsf(d2high - d2low) > d2tol; i++) {
+		d2 = 0.5f * (d2low + d2high);
+
+		/* find where we are after n-1 chords of squared length d2 */
+		u = 0.0f;
+		for (j = 0; j < n - 1; j++) {
+			u = find_superellipse_chord_u(u, d2, r);
+			if (u == -1.0f)
+				break;  /* d2 is too big to go n-1 chords */
+		}
+		if (u == -1.0f) {
+			d2high = d2;
+			continue;
+		}
+		d2final = superellipse_chord_length_squared(u, umax, r);
+		if (fabsf(d2final - d2) <= d2tol)
+			break;
+		if (d2final < d2)
+			d2high = d2;
+		else
+			d2low = d2;
+	}
+	u = 0.0f;
+	for (i = 0; i < n; i++) {
+		r_params[i] = u;
+		u = find_superellipse_chord_u(u, d2, r);
+	}
+	r_params[n] = umax;
+}
+
+/* The superellipse used for multisegment profiles does not
+ * have a closed-form way to generate evenly spaced points
+ * along an arc. We use an expensive search procedure to find
+ * the parameter values that lead to bp->seg even chords.
+ * We also want spacing for a number of segments that is
+ * a power of 2 >= bp->seg (but at least 4). */
+static void set_profile_spacing(BevelParams *bp)
+{
+	int seg, seg_2;
+
+	seg = bp->seg;
+	if (seg > 1) {
+		bp->pro_spacing.uvals = (float *)BLI_memarena_alloc(bp->mem_arena, (seg + 1) * sizeof(float));
+		find_even_superellipse_params(seg, bp->pro_super_r, bp->pro_spacing.uvals);
+		seg_2 = power_of_2_max_i(bp->seg);
+		if (seg_2 == 2)
+			seg_2 = 4;
+		bp->pro_spacing.seg_2 = seg_2;
+		if (seg_2 == seg) {
+			bp->pro_spacing.uvals_2 = bp->pro_spacing.uvals;
+		}
+		else {
+			bp->pro_spacing.uvals_2 = (float *)BLI_memarena_alloc(bp->mem_arena, (seg_2 + 1) * sizeof(float));
+			find_even_superellipse_params(seg_2, bp->pro_super_r, bp->pro_spacing.uvals_2);
+		}
+	}
+	else {
+		bp->pro_spacing.uvals = NULL;
+		bp->pro_spacing.uvals_2 = NULL;
+		bp->pro_spacing.seg_2 = 0;
+	}
+}
+
 /*
  * Calculate and return an offset that is the lesser of the current
  * bp.offset and the maximum possible offset before geometry
@@ -3568,6 +3839,7 @@ void BM_mesh_bevel(BMesh *bm, const float offset, const int offset_type,
 		bp.vert_hash = BLI_ghash_ptr_new(__func__);
 		bp.mem_arena = BLI_memarena_new(MEM_SIZE_OPTIMAL(1 << 16), __func__);
 		BLI_memarena_use_calloc(bp.mem_arena);
+		set_profile_spacing(&bp);
 
 		if (limit_offset)
 			bp.offset = bevel_limit_offset(bm, &bp);