Fix Bevel intersection continuity.

This patch from Henrik Dick improves the continuity between the
grid forming corners and the edge polyons on multisegment bevels.
For details, see patch D13867.

source/blender/bmesh/tools/bmesh_bevel.c

@@ -451,6 +451,16 @@ static bool nearly_parallel_normalized(const float d1[3], const float d2[3])
   return compare_ff(fabsf(direction_dot), 1.0f, BEVEL_EPSILON_ANG_DOT);
 }
 
+/**
+ * calculate the determinant of a matrix formed by three vectors
+ * \return dot(a, cross(b, c)) = determinant(a, b, c)
+ */
+static float determinant_v3v3v3(const float a[3], const float b[3], const float c[3])
+{
+  return a[0] * b[1] * c[2] + a[1] * b[2] * c[0] + a[2] * b[0] * c[1] - a[0] * b[2] * c[1] -
+         a[1] * b[0] * c[2] - a[2] * b[1] * c[0];
+}
+
 /* Make a new BoundVert of the given kind, inserting it at the end of the circular linked
  * list with entry point bv->boundstart, and return it. */
 static BoundVert *add_new_bound_vert(MemArena *mem_arena, VMesh *vm, const float co[3])
@@ -4118,44 +4128,114 @@ static VMesh *cubic_subdiv(BevelParams *bp, VMesh *vm_in)
   VMesh *vm_out = new_adj_vmesh(bp->mem_arena, n_boundary, ns_out, vm_in->boundstart);
 
   /* First we adjust the boundary vertices of the input mesh, storing in output mesh. */
+  BoundVert *bndv = vm_in->boundstart;
   for (int i = 0; i < n_boundary; i++) {
+    float co1[3], co2[3], acc[3];
+    EdgeHalf *e = bndv->elast;
+    /* Generate tangents. This is hacked together and would ideally be done elsewere and then only
+     * used here. */
+    float tangent[3], tangent2[3], normal[3];
+    bool convex = true;
+    bool orthogonal = false;
+    float stretch = 0.0f;
+    if (e) {
+      /* Projection direction is direction of the edge. */
+      sub_v3_v3v3(tangent, e->e->v1->co, e->e->v2->co);
+      if (e->is_rev) {
+        negate_v3(tangent);
+      }
+      normalize_v3(tangent);
+      if (bndv->is_arc_start || bndv->is_patch_start) {
+        BMFace *face = e->fnext;
+        if (face) {
+          copy_v3_v3(normal, face->no);
+        }
+        else {
+          zero_v3(normal);
+        }
+        madd_v3_v3v3fl(co2, bndv->profile.middle, normal, 0.1f);
+      }
+      if (bndv->is_arc_start || bp->affect_type == BEVEL_AFFECT_VERTICES) {
+        EdgeHalf *e1 = bndv->next->elast;
+        BLI_assert(e1);
+        sub_v3_v3v3(tangent2, e1->e->v1->co, e1->e->v2->co);
+        if (e1->is_rev) {
+          negate_v3(tangent2);
+        }
+        normalize_v3(tangent2);
+
+        convex = determinant_v3v3v3(tangent2, tangent, normal) < 0;
+
+        add_v3_v3(tangent2, tangent);
+        normalize_v3(tangent2);
+        copy_v3_v3(tangent, tangent2);
+      }
+      /* Calculate a factor which determines how much the interpolated mesh is
+       * going to be stretched out into the direction of the tangent.
+       * It is currently using the difference along the tangent of the
+       * central point on the profile and the current center vertex position. */
+      get_profile_point(bp, &bndv->profile, ns_in2, ns_in, co);
+      stretch = dot_v3v3(tangent, mesh_vert(vm_in, i, ns_in2, ns_in2)->co) - dot_v3v3(tangent, co);
+      stretch = fabsf(stretch);
+      /* Scale the tangent by stretch. The divide by ns_in2 comes from the Levin Paper. */
+      mul_v3_fl(tangent, stretch / ns_in2);
+      orthogonal = bndv->is_patch_start;
+    }
+    else if (bndv->prev->is_patch_start) {
+      /* If this is the second edge of a patch and therefore #e is NULL,
+       * then e->fprev has to be used/not NULL. */
+      BLI_assert(bndv->prev->elast);
+      BMFace *face = bndv->prev->elast->fnext;
+      if (face) {
+        copy_v3_v3(normal, face->no);
+      }
+      else {
+        zero_v3(normal);
+      }
+      orthogonal = true;
+    }
+    else {
+      /** Should only come here from make_cube_corner_adj_vmesh. */
+      sub_v3_v3v3(co1, mesh_vert(vm_in, i, 0, 0)->co, mesh_vert(vm_in, i, 0, 1)->co);
+      sub_v3_v3v3(co2, mesh_vert(vm_in, i, 0, 1)->co, mesh_vert(vm_in, i, 0, 2)->co);
+      cross_v3_v3v3(tangent, co1, co2);
+      /** The following constant is choosen to best match the old results. */
+      normalize_v3_length(tangent, 1.5f / ns_out);
+    }
+    /** Copy corner vertex. */
     copy_v3_v3(mesh_vert(vm_out, i, 0, 0)->co, mesh_vert(vm_in, i, 0, 0)->co);
+    /** Copy the rest of the boundary vertices. */
     for (int k = 1; k < ns_in; k++) {
       copy_v3_v3(co, mesh_vert(vm_in, i, 0, k)->co);
 
-      /* Smooth boundary rule. Custom profiles shouldn't be smoothed. */
-      if (bp->profile_type != BEVEL_PROFILE_CUSTOM) {
-        float co1[3], co2[3], acc[3];
-        copy_v3_v3(co1, mesh_vert(vm_in, i, 0, k - 1)->co);
-        copy_v3_v3(co2, mesh_vert(vm_in, i, 0, k + 1)->co);
+      copy_v3_v3(co1, mesh_vert(vm_in, i, 0, k - 1)->co);
+      copy_v3_v3(co2, mesh_vert(vm_in, i, 0, k + 1)->co);
 
-        add_v3_v3v3(acc, co1, co2);
-        madd_v3_v3fl(acc, co, -2.0f);
-        madd_v3_v3fl(co, acc, -1.0f / 6.0f);
+      add_v3_v3v3(acc, co1, co2);
+      if (bndv->is_arc_start) {
+        sub_v3_v3(co1, co);
+        sub_v3_v3(co2, co);
+        normalize_v3(co1);
+        normalize_v3(co2);
+        add_v3_v3v3(tangent, co1, co2);
+        /* This is an empirical formula to make the result look good. */
+        normalize_v3(tangent);
+        float dot = convex ? fminf(0, dot_v3v3(tangent2, tangent)) : 1.0f;
+        mul_v3_fl(tangent, stretch / ns_in * dot);
       }
+      else if (orthogonal) {
+        sub_v3_v3(co1, co);
+        cross_v3_v3v3(tangent, normal, co1);
+        /* This is an empirical formula to make the result look good. */
+        normalize_v3_length(tangent, -bp->offset * 0.7071f / ns_in);
+      }
+      mul_v3_fl(co, 2.0f);
+      madd_v3_v3fl(co, acc, -0.25f);
+      madd_v3_v3fl(co, mesh_vert(vm_in, i, 1, k)->co, -0.5f);
+      add_v3_v3(co, tangent);
 
       copy_v3_v3(mesh_vert_canon(vm_out, i, 0, 2 * k)->co, co);
     }
-  }
-  /* Now adjust odd boundary vertices in output mesh, based on even ones. */
-  BoundVert *bndv = vm_out->boundstart;
-  for (int i = 0; i < n_boundary; i++) {
-    for (int k = 1; k < ns_out; k += 2) {
-      get_profile_point(bp, &bndv->profile, k, ns_out, co);
-
-      /* Smooth if using a non-custom profile. */
-      if (bp->profile_type != BEVEL_PROFILE_CUSTOM) {
-        float co1[3], co2[3], acc[3];
-        copy_v3_v3(co1, mesh_vert_canon(vm_out, i, 0, k - 1)->co);
-        copy_v3_v3(co2, mesh_vert_canon(vm_out, i, 0, k + 1)->co);
-
-        add_v3_v3v3(acc, co1, co2);
-        madd_v3_v3fl(acc, co, -2.0f);
-        madd_v3_v3fl(co, acc, -1.0f / 6.0f);
-      }
-
-      copy_v3_v3(mesh_vert_canon(vm_out, i, 0, k)->co, co);
-    }
     bndv = bndv->next;
   }
   vmesh_copy_equiv_verts(vm_out);
@@ -4163,7 +4243,7 @@ static VMesh *cubic_subdiv(BevelParams *bp, VMesh *vm_in)
   /* Copy adjusted verts back into vm_in. */
   for (int i = 0; i < n_boundary; i++) {
     for (int k = 0; k < ns_in; k++) {
-      copy_v3_v3(mesh_vert(vm_in, i, 0, k)->co, mesh_vert(vm_out, i, 0, 2 * k)->co);
+      copy_v3_v3(mesh_vert_canon(vm_in, i, 0, k)->co, mesh_vert_canon(vm_out, i, 0, 2 * k)->co);
     }
   }
 
@@ -4248,7 +4328,7 @@ static VMesh *cubic_subdiv(BevelParams *bp, VMesh *vm_in)
   vmesh_copy_equiv_verts(vm_out);
 
   /* The center vertex is special. */
-  gamma = sabin_gamma(n_boundary);
+  gamma = sabin_gamma(n_boundary) * 0.5f;
   beta = -gamma;
   /* Accumulate edge verts in co1, face verts in co2. */
   float co1[3], co2[3];