Anim: implement a new curve-aware vertex to B-Bone segment mapping mode.

Currently vertices are mapped to B-Bone segments without taking the
rest pose curve into account. This is very simple and fast, but causes
poor deformations in some cases where the rest curvature is significant,
e.g. mouth corners in the new Rigify face rig.

This patch implements a new mapping mode that addresses this problem.
The general idea is to do an orthogonal projection on the curve. However,
since there is no analytical solution for bezier curves, it uses the
segment approximation:

* First, boundaries between segments are used for a binary space
  partitioning search to narrow down the mapping to one segment.
* Then, a position on the segment is chosen via linear
  interpolation between the BSP planes.
* Finally, to remove the sharp discontinuity at the evolute surface
  a smoothing pass is applied to the chosen position by blending to
  reduce the slope around the planes previously used in the BSP search.

In order to make per-vertex processing faster, a new array with the
necessary vectors converted to the pose space, as well as some
precomputed coefficients, is built.

The new mode is implemented as a per-bone option in order to ensure
backward compatibility, and also because the new mode may not be
optimal for all cases due to the difference in performance, and
complications like the smoothed but still present mapping
discontinuities around the evolute surface.

Wiki: https://wiki.blender.org/wiki/Source/Animation/B-Bone_Vertex_Mapping

Pull Request: https://projects.blender.org/blender/blender/pulls/110758

Pull Request: https://projects.blender.org/blender/blender/pulls/110758

scripts/startup/bl_ui/properties_data_bone.py

@@ -140,6 +140,8 @@ class BONE_PT_curved(BoneButtonsPanel, Panel):
         topcol = layout.column()
         topcol.active = bone.bbone_segments > 1
 
+        topcol.prop(bone, "bbone_mapping_mode", text="Vertex Mapping")
+
         col = topcol.column(align=True)
         col.prop(bbone, "bbone_curveinx", text="Curve In X")
         col.prop(bbone, "bbone_curveinz", text="Z")

source/blender/blenkernel/BKE_armature.h

@@ -81,6 +81,8 @@ typedef struct EditBone {
   /** for envelope scaling */
   float oldlength;
 
+  /** Mapping of vertices to segments. */
+  eBone_BBoneMappingMode bbone_mapping_mode;
   /** Type of next/prev bone handles */
   char bbone_prev_type;
   char bbone_next_type;

source/blender/blenkernel/intern/action.cc

@@ -1048,6 +1048,7 @@ void BKE_pose_channel_free_bbone_cache(bPoseChannel_Runtime *runtime)
   MEM_SAFE_FREE(runtime->bbone_pose_mats);
   MEM_SAFE_FREE(runtime->bbone_deform_mats);
   MEM_SAFE_FREE(runtime->bbone_dual_quats);
+  MEM_SAFE_FREE(runtime->bbone_segment_boundaries);
 }
 
 void BKE_pose_channel_free(bPoseChannel *pchan)

source/blender/blenkernel/intern/armature.cc

@@ -54,6 +54,7 @@
 #include "DEG_depsgraph_query.hh"
 
 #include "BIK_api.h"
+#include "BLI_math_base_safe.h"
 
 #include "BLO_read_write.hh"
 
@@ -1509,7 +1510,9 @@ int BKE_pchan_bbone_spline_compute(BBoneSplineParameters *param,
   return param->segments;
 }
 
-static void allocate_bbone_cache(bPoseChannel *pchan, int segments)
+static void allocate_bbone_cache(bPoseChannel *pchan,
+                                 const int segments,
+                                 const bool use_boundaries)
 {
   bPoseChannel_Runtime *runtime = &pchan->runtime;
 
@@ -1526,6 +1529,69 @@ static void allocate_bbone_cache(bPoseChannel *pchan, int segments)
     runtime->bbone_dual_quats = static_cast<DualQuat *>(MEM_malloc_arrayN(
         1 + uint(segments), sizeof(DualQuat), "bPoseChannel_Runtime::bbone_dual_quats"));
   }
+
+  /* If the segment count changed, the array was deallocated and nulled above. */
+  if (use_boundaries && !runtime->bbone_segment_boundaries) {
+    runtime->bbone_segment_boundaries = static_cast<bPoseChannel_BBoneSegmentBoundary *>(
+        MEM_malloc_arrayN(1 + uint(segments),
+                          sizeof(bPoseChannel_BBoneSegmentBoundary),
+                          "bPoseChannel_Runtime::bbone_segment_boundaries"));
+  }
+  else if (!use_boundaries) {
+    MEM_SAFE_FREE(runtime->bbone_segment_boundaries);
+  }
+}
+
+/** Computes the B-Bone segment boundary planes for the curved mapping. */
+static void compute_bbone_segment_boundaries(bPoseChannel *pchan)
+{
+  const Bone *bone = pchan->bone;
+  bPoseChannel_Runtime *runtime = &pchan->runtime;
+  const Mat4 *b_bone_rest = runtime->bbone_rest_mats;
+  bPoseChannel_BBoneSegmentBoundary *boundaries = runtime->bbone_segment_boundaries;
+
+  /* Convert joints to pose space. */
+  for (int i = 0; i <= bone->segments; i++) {
+    mul_v3_m4v3(boundaries[i].point, bone->arm_mat, b_bone_rest[i].mat[3]);
+    mul_v3_mat3_m4v3(boundaries[i].plane_normal, bone->arm_mat, b_bone_rest[i].mat[1]);
+    normalize_v3(boundaries[i].plane_normal);
+  }
+
+  /* Precompute coefficients for the mapping calculations. */
+  for (int i = 0; i <= bone->segments; i++) {
+    boundaries[i].plane_offset = dot_v3v3(boundaries[i].point, boundaries[i].plane_normal);
+  }
+
+  /* Precompute the inverted length of the curve. */
+  float arc_length = 0.0f;
+
+  for (int i = 0; i < bone->segments; i++) {
+    arc_length += len_v3v3(boundaries[i + 1].point, boundaries[i].point);
+  }
+
+  runtime->bbone_arc_length_reciprocal = 1.0f / arc_length;
+
+  /* Precompute the BSP depth based widening coefficients.
+   * The actual space partitioning includes two extra virtual segments for the ends. */
+  const int bsp_depth = int(ceilf(log2f(bone->segments + 2)));
+
+  BLI_assert(bsp_depth <= bone->segments);
+
+  /* Maximum half-width of the smoothing band at the bsp tree root plane, in segments.
+   * The tuning coefficient was chosen by trial and error (see PR #110758). */
+  const float tuning_factor = 0.222f;
+  const float straight_length = len_v3v3(boundaries[0].point, boundaries[bone->segments].point);
+  const float max_depth_scale = bone->segments * (straight_length / arc_length) * tuning_factor;
+
+  /* Per tree layer scaling factor, aiming to reduce the radius to 1 segment at the leaf level.
+   * Since depth_scale is actually a reciprocal of the width, this factor is >= 1. */
+  const float scale_factor = powf(max_ff(max_depth_scale, 1.0f), 1.0f / (bsp_depth - 1));
+
+  boundaries[0].depth_scale = bone->segments / max_depth_scale;
+
+  for (int i = 1; i < bsp_depth; i++) {
+    boundaries[i].depth_scale = boundaries[i - 1].depth_scale * scale_factor;
+  }
 }
 
 void BKE_pchan_bbone_segments_cache_compute(bPoseChannel *pchan)
@@ -1537,7 +1603,9 @@ void BKE_pchan_bbone_segments_cache_compute(bPoseChannel *pchan)
   BLI_assert(segments > 1);
 
   /* Allocate the cache if needed. */
-  allocate_bbone_cache(pchan, segments);
+  const bool use_curved_mapping = bone->bbone_mapping_mode == BBONE_MAPPING_CURVED;
+
+  allocate_bbone_cache(pchan, segments, use_curved_mapping);
 
   /* Compute the shape. */
   Mat4 *b_bone = runtime->bbone_pose_mats;
@@ -1549,6 +1617,12 @@ void BKE_pchan_bbone_segments_cache_compute(bPoseChannel *pchan)
   BKE_pchan_bbone_spline_setup(pchan, false, true, b_bone);
   BKE_pchan_bbone_spline_setup(pchan, true, true, b_bone_rest);
 
+  /* Compute segment boundaries. */
+  if (runtime->bbone_segment_boundaries) {
+    BLI_assert(use_curved_mapping);
+    compute_bbone_segment_boundaries(pchan);
+  }
+
   /* Compute deform matrices. */
   /* first matrix is the inverse arm_mat, to bring points in local bone space
    * for finding out which segment it belongs to */
@@ -1589,7 +1663,9 @@ void BKE_pchan_bbone_segments_cache_copy(bPoseChannel *pchan, bPoseChannel *pcha
     BKE_pose_channel_free_bbone_cache(&pchan->runtime);
   }
   else {
-    allocate_bbone_cache(pchan, segments);
+    const bool use_curved_mapping = runtime_from->bbone_segment_boundaries != nullptr;
+
+    allocate_bbone_cache(pchan, segments, use_curved_mapping);
 
     memcpy(runtime->bbone_rest_mats, runtime_from->bbone_rest_mats, sizeof(Mat4) * (1 + segments));
     memcpy(runtime->bbone_pose_mats, runtime_from->bbone_pose_mats, sizeof(Mat4) * (1 + segments));
@@ -1599,6 +1675,17 @@ void BKE_pchan_bbone_segments_cache_copy(bPoseChannel *pchan, bPoseChannel *pcha
     memcpy(runtime->bbone_dual_quats,
            runtime_from->bbone_dual_quats,
            sizeof(DualQuat) * (1 + segments));
+
+    if (use_curved_mapping) {
+      runtime->bbone_arc_length_reciprocal = runtime_from->bbone_arc_length_reciprocal;
+
+      memcpy(runtime->bbone_segment_boundaries,
+             runtime_from->bbone_segment_boundaries,
+             sizeof(bPoseChannel_BBoneSegmentBoundary) * (1 + segments));
+    }
+    else {
+      BLI_assert(runtime->bbone_segment_boundaries == nullptr);
+    }
   }
 }
 
@@ -1627,10 +1714,11 @@ void BKE_pchan_bbone_deform_clamp_segment_index(const bPoseChannel *pchan,
   *r_blend_next = blend;
 }
 
-void BKE_pchan_bbone_deform_segment_index(const bPoseChannel *pchan,
-                                          const float *co,
-                                          int *r_index,
-                                          float *r_blend_next)
+/** Implementation of the Straight B-Bone segment mapping. */
+static void find_bbone_segment_index_straight(const bPoseChannel *pchan,
+                                              const float *co,
+                                              int *r_index,
+                                              float *r_blend_next)
 {
   const Mat4 *mats = pchan->runtime.bbone_deform_mats;
   const float(*mat)[4] = mats[0].mat;
@@ -1643,6 +1731,147 @@ void BKE_pchan_bbone_deform_segment_index(const bPoseChannel *pchan,
       pchan, y / pchan->bone->length, r_index, r_blend_next);
 }
 
+/** Computes signed distance to the segment boundary BSP plane. */
+inline float bbone_segment_bsp_signed_distance(const bPoseChannel_BBoneSegmentBoundary &boundary,
+                                               const float *co)
+{
+  return dot_v3v3(co, boundary.plane_normal) - boundary.plane_offset;
+}
+
+/** Implementation of the Curved B-Bone segment mapping. */
+static void find_bbone_segment_index_curved(const bPoseChannel *pchan,
+                                            const float *co,
+                                            int *r_index,
+                                            float *r_blend_next)
+{
+  const bPoseChannel_BBoneSegmentBoundary *boundaries = pchan->runtime.bbone_segment_boundaries;
+  const int segments = pchan->runtime.bbone_segments;
+
+  /* Saved signed distances from co to each checked boundary plane. */
+  float boundary_dist[MAX_BBONE_SUBDIV + 1];
+
+  /* Stack of BSP plane indices that were checked in the binary search. */
+  int boundary_idx_stack[MAX_BBONE_SUBDIV];
+  int stack_top = -1;
+
+  /* Perform a BSP binary search to narrow things down to one segment.
+   * Checked BSP planes are stored for the smoothing pass later. */
+  int start = -1, end = segments + 1, bias = 0;
+
+  while (end - start > 1) {
+    const int mid = (start + end + bias) / 2;
+
+    BLI_assert(start < mid && mid < end);
+
+    const float dist = bbone_segment_bsp_signed_distance(boundaries[mid], co);
+
+    boundary_idx_stack[++stack_top] = mid;
+    boundary_dist[mid] = dist;
+
+    if (dist < 0) {
+      end = mid;
+      /* Bias division of odd numbers toward the previous split. This should produce
+       * a slightly smoother and more symmetrical evolute boundary near the ends. */
+      bias = 1;
+    }
+    else {
+      start = mid;
+      bias = 0;
+    }
+  }
+
+  /* Compute the mapping from the individual segment, or the curve ends. */
+  const float segment_size = 1.0f / segments;
+  float head_tail;
+
+  if (end <= 0) {
+    head_tail = 0;
+  }
+  else if (start >= segments) {
+    head_tail = 1;
+  }
+  else {
+    /* Linear interpolation between the innermost two planes. */
+    const float d1 = fabsf(boundary_dist[start]);
+    const float d2 = fabsf(boundary_dist[end]);
+    const float t = d1 / (d1 + d2);
+
+    head_tail = segment_size * (start + t);
+  }
+
+  /* Smooth the mapping to suppress discontinuities by using BSP boundaries up the stack.
+   *
+   * This works basically by pulling the mapped position towards the boundary in order to
+   * reduce the gradient slope to the ideal value (the one you get for points directly on
+   * the curve), using heuristic blend strength falloff coefficients based on the distances
+   * to the boundary plane before and after mapping. See PR #110758 for more details, or
+   * https://wiki.blender.org/wiki/Source/Animation/B-Bone_Vertex_Mapping#Curved_Mapping */
+  const float segment_scale = pchan->runtime.bbone_arc_length_reciprocal;
+
+  for (int i = stack_top; i >= 0; --i) {
+    const int boundary_idx = boundary_idx_stack[i];
+
+    /* Boundary in the head-tail space. */
+    const float boundary_pos = boundary_idx * segment_size;
+
+    /* Distance of the original 3d point (co) from the boundary plane,
+     * mapped to the head-tail space using the ideal slope ratio. */
+    const float point_dist = boundary_dist[boundary_idx] * segment_scale;
+    const float point_dist_abs = fabsf(point_dist);
+
+    /* Distance of the current mapped position from the boundary in the head-tail space. */
+    const float mapped_dist = fabsf(head_tail - boundary_pos);
+
+    /* Only reduce the local gradient slope, don't increase it. This basically limits
+     * smoothing to the inside of the curve, leaving outside as is. */
+    const float slope_gap = mapped_dist - point_dist_abs;
+
+    if (slope_gap <= 0) {
+      continue;
+    }
+
+    /* Only affect points close to the split line; the radius depends on the depth
+     * in the stack using precomputed coefficients. */
+    const float dist_coeff = 1.0f - point_dist_abs * boundaries[i].depth_scale;
+
+    if (dist_coeff <= 0) {
+      continue;
+    }
+
+    /* Asymptotically clamp the slope coefficient to 1. The tune coefficients here and
+     * below control the sharpness of the transition and were chosen by trial and error. */
+    const float slope_tune_coeff = 3.0f;
+    const float scaled_gap = slope_gap * slope_tune_coeff;
+    const float slope_coeff = scaled_gap / (scaled_gap + point_dist_abs);
+
+    /* Smooth the distance based coefficient around zero. */
+    const float dist_tune_coeff = 7.0f;
+    const float dist_coeff_smooth = dist_coeff * dist_coeff * (dist_tune_coeff + 1.0f) /
+                                    (dist_tune_coeff * dist_coeff + 1.0f);
+
+    /* Blend towards the point on the ideal slope. */
+    const float target_pos = boundary_pos + point_dist;
+
+    head_tail = interpf(target_pos, head_tail, slope_coeff * dist_coeff_smooth);
+  }
+
+  /* Calculate the indices of the 2 affecting b_bone segments. */
+  BKE_pchan_bbone_deform_clamp_segment_index(pchan, head_tail, r_index, r_blend_next);
+}
+
+void BKE_pchan_bbone_deform_segment_index(const bPoseChannel *pchan,
+                                          const float *co,
+                                          int *r_index,
+                                          float *r_blend_next)
+{
+  if (pchan->runtime.bbone_segment_boundaries) {
+    find_bbone_segment_index_curved(pchan, co, r_index, r_blend_next);
+  }
+  else {
+    find_bbone_segment_index_straight(pchan, co, r_index, r_blend_next);
+  }
+}
+
 /** \} */
 
 /* -------------------------------------------------------------------- */

source/blender/editors/armature/armature_add.cc

@@ -1313,6 +1313,7 @@ static int armature_symmetrize_exec(bContext *C, wmOperator *op)
         ebone->bbone_prev_type = ebone_iter->bbone_prev_type;
         ebone->bbone_next_type = ebone_iter->bbone_next_type;
 
+        ebone->bbone_mapping_mode = ebone_iter->bbone_mapping_mode;
         ebone->bbone_flag = ebone_iter->bbone_flag;
         ebone->bbone_prev_flag = ebone_iter->bbone_prev_flag;
         ebone->bbone_next_flag = ebone_iter->bbone_next_flag;

source/blender/editors/armature/armature_utils.cc

@@ -519,6 +519,7 @@ static EditBone *make_boneList_recursive(ListBase *edbo,
     eBone->bbone_prev_type = curBone->bbone_prev_type;
     eBone->bbone_next_type = curBone->bbone_next_type;
 
+    eBone->bbone_mapping_mode = eBone_BBoneMappingMode(curBone->bbone_mapping_mode);
     eBone->bbone_flag = curBone->bbone_flag;
     eBone->bbone_prev_flag = curBone->bbone_prev_flag;
     eBone->bbone_next_flag = curBone->bbone_next_flag;
@@ -733,6 +734,7 @@ void ED_armature_from_edit(Main *bmain, bArmature *arm)
     newBone->bbone_prev_type = eBone->bbone_prev_type;
     newBone->bbone_next_type = eBone->bbone_next_type;
 
+    newBone->bbone_mapping_mode = eBone->bbone_mapping_mode;
     newBone->bbone_flag = eBone->bbone_flag;
     newBone->bbone_prev_flag = eBone->bbone_prev_flag;
     newBone->bbone_next_flag = eBone->bbone_next_flag;

source/blender/makesdna/DNA_action_types.h

@@ -180,6 +180,20 @@ typedef struct bPoseChannelDrawData {
 struct DualQuat;
 struct Mat4;
 
+/* Describes a plane in pose space that delimits B-Bone segments. */
+typedef struct bPoseChannel_BBoneSegmentBoundary {
+  /* Boundary data in pose space. */
+  float point[3];
+  float plane_normal[3];
+  /* Dot product of point and plane_normal to speed up distance computation. */
+  float plane_offset;
+
+  /* Inverse width of the smoothing at this level in head-tail space.
+   * Optimization: this value is actually indexed by bsp depth (0 to bsp_depth-1), not joint
+   * index. It's put here to avoid allocating a separate array by utilizing the padding space. */
+  float depth_scale;
+} bPoseChannel_BBoneSegmentBoundary;
+
 typedef struct bPoseChannel_Runtime {
   SessionUUID session_uuid;
 
@@ -189,6 +203,10 @@ typedef struct bPoseChannel_Runtime {
   /* B-Bone shape data: copy of the segment count for validation. */
   int bbone_segments;
 
+  /* Inverse of the total length of the segment polyline. */
+  float bbone_arc_length_reciprocal;
+  char _pad1[4];
+
   /* Rest and posed matrices for segments. */
   struct Mat4 *bbone_rest_mats;
   struct Mat4 *bbone_pose_mats;
@@ -196,6 +214,10 @@ typedef struct bPoseChannel_Runtime {
   /* Delta from rest to pose in matrix and DualQuat form. */
   struct Mat4 *bbone_deform_mats;
   struct DualQuat *bbone_dual_quats;
+
+  /* Segment boundaries for curved mode. */
+  struct bPoseChannel_BBoneSegmentBoundary *bbone_segment_boundaries;
+  void *_pad;
 } bPoseChannel_Runtime;
 
 /* ************************************************ */

source/blender/makesdna/DNA_armature_types.h

@@ -122,6 +122,9 @@ typedef struct Bone {
   int layer;
   /** For B-bones. */
   short segments;
+  /** Vertex to segment mapping mode. */
+  char bbone_mapping_mode;
+  char _pad2[7];
 
   /** Type of next/prev bone handles. */
   char bbone_prev_type;
@@ -395,6 +398,12 @@ typedef enum eBone_BBoneHandleType {
   BBONE_HANDLE_TANGENT = 3,  /* Custom handle in tangent mode (use direction, not location). */
 } eBone_BBoneHandleType;
 
+/* bone->bbone_mapping_mode */
+typedef enum eBone_BBoneMappingMode {
+  BBONE_MAPPING_STRAIGHT = 0,  /* Default mode that ignores the rest pose curvature. */
+  BBONE_MAPPING_CURVED = 1,    /* Mode that takes the rest pose curvature into account. */
+} eBone_BBoneMappingMode;
+
 /* bone->bbone_flag */
 typedef enum eBone_BBoneFlag {
   /** Add the parent Out roll to the In roll. */

source/blender/makesrna/intern/rna_armature.cc

@@ -1104,6 +1104,22 @@ static void rna_def_bone_common(StructRNA *srna, int editbone)
       {0, nullptr, 0, nullptr, nullptr},
   };
 
+  static const EnumPropertyItem prop_bbone_mapping_mode[] = {
+      {BBONE_MAPPING_STRAIGHT,
+       "STRAIGHT",
+       0,
+       "Straight",
+       "Fast mapping that is good for most situations, but ignores the rest pose "
+       "curvature of the B-Bone"},
+      {BBONE_MAPPING_CURVED,
+       "CURVED",
+       0,
+       "Curved",
+       "Slower mapping that gives better deformation for B-Bones that are sharply "
+       "curved in rest pose"},
+      {0, nullptr, 0, nullptr, nullptr},
+  };
+
   static const EnumPropertyItem prop_inherit_scale_mode[] = {
       {BONE_INHERIT_SCALE_FULL, "FULL", 0, "Full", "Inherit all effects of parent scaling"},
       {BONE_INHERIT_SCALE_FIX_SHEAR,
@@ -1292,6 +1308,16 @@ static void rna_def_bone_common(StructRNA *srna, int editbone)
   RNA_def_property_ui_text(
       prop, "B-Bone Segments", "Number of subdivisions of bone (for B-Bones only)");
 
+  prop = RNA_def_property(srna, "bbone_mapping_mode", PROP_ENUM, PROP_NONE);
+  RNA_def_property_enum_sdna(prop, NULL, "bbone_mapping_mode");
+  RNA_def_property_enum_items(prop, prop_bbone_mapping_mode);
+  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  RNA_def_property_ui_text(
+      prop,
+      "B-Bone Vertex Mapping Mode",
+      "Selects how the vertices are mapped to B-Bone segments based on their position");
+  RNA_def_property_update(prop, 0, "rna_Armature_update_data");
+
   prop = RNA_def_property(srna, "bbone_x", PROP_FLOAT, PROP_NONE);
   if (editbone) {
     RNA_def_property_update(prop, 0, "rna_Armature_editbone_transform_update");