Curves: support deforming curves on surface

Curves that are attached to a surface can now follow the surface when
it is modified using shape keys or modifiers (but not when the original
surface is deformed in edit or sculpt mode).

The surface is allowed to be changed in any way that keeps uv maps
intact. So deformation is allowed, but also some topology changes like
subdivision.

The following features are added:
* A new `Deform Curves on Surface` node, which deforms curves with
  attachment information based on the surface object and uv map set
  in the properties panel.
* A new `Add Rest Position` checkbox in the shape keys panel. When checked,
  a new `rest_position` vector attribute is added to the mesh before shape
  keys and modifiers are applied. This is necessary to support proper
  deformation of the curves, but can also be used for other purposes.
* The `Add > Curve > Empty Hair` operator now sets up a simple geometry
  nodes setup that deforms the hair. It also makes sure that the rest
  position attribute is added to the surface.
* A new `Object (Attach Curves to Surface)` operator in the `Set Parent To`
  (ctrl+P) menu, which attaches existing curves to the surface and sets the
  surface object as parent.

Limitations:
* Sculpting the procedurally deformed curves will be implemented separately.
* The `Deform Curves on Surface` node is not generic and can only be used
  for one specific purpose currently. We plan to generalize this more in the
  future by adding support by exposing more inputs and/or by turning it into
  a node group.

Differential Revision: https://developer.blender.org/D14864

4 files changed, 367 insertions, 0 deletions

diff --git a/source/blender/nodes/NOD_geometry.h b/source/blender/nodes/NOD_geometry.h
index 8f15add33fd..86c276fbd6f 100644
--- a/source/blender/nodes/NOD_geometry.h
+++ b/source/blender/nodes/NOD_geometry.h
@@ -47,6 +47,7 @@ void register_node_type_geo_curve_subdivide(void);
 void register_node_type_geo_curve_to_mesh(void);
 void register_node_type_geo_curve_to_points(void);
 void register_node_type_geo_curve_trim(void);
+void register_node_type_geo_deform_curves_on_surface(void);
 void register_node_type_geo_delete_geometry(void);
 void register_node_type_geo_duplicate_elements(void);
 void register_node_type_geo_distribute_points_on_faces(void);
diff --git a/source/blender/nodes/NOD_static_types.h b/source/blender/nodes/NOD_static_types.h
index 302656b8a7c..7c7f114bb78 100644
--- a/source/blender/nodes/NOD_static_types.h
+++ b/source/blender/nodes/NOD_static_types.h
@@ -301,6 +301,7 @@ DefNode(GeometryNode, GEO_NODE_CURVE_SPLINE_PARAMETER, 0, "SPLINE_PARAMETER", Sp
 DefNode(GeometryNode, GEO_NODE_CURVE_SPLINE_TYPE, def_geo_curve_spline_type, "CURVE_SPLINE_TYPE", CurveSplineType, "Set Spline Type", "")
 DefNode(GeometryNode, GEO_NODE_CURVE_TO_MESH, 0, "CURVE_TO_MESH", CurveToMesh, "Curve to Mesh", "")
 DefNode(GeometryNode, GEO_NODE_CURVE_TO_POINTS, def_geo_curve_to_points, "CURVE_TO_POINTS", CurveToPoints, "Curve to Points", "")
+DefNode(GeometryNode, GEO_NODE_DEFORM_CURVES_ON_SURFACE, 0, "DEFORM_CURVES_ON_SURFACE", DeformCurvesOnSurface, "Deform Curves on Surface", "")
 DefNode(GeometryNode, GEO_NODE_DELETE_GEOMETRY, def_geo_delete_geometry, "DELETE_GEOMETRY", DeleteGeometry, "Delete Geometry", "")
 DefNode(GeometryNode, GEO_NODE_DUPLICATE_ELEMENTS, def_geo_duplicate_elements, "DUPLICATE_ELEMENTS", DuplicateElements, "Duplicate Elements", "")
 DefNode(GeometryNode, GEO_NODE_DISTRIBUTE_POINTS_ON_FACES, def_geo_distribute_points_on_faces, "DISTRIBUTE_POINTS_ON_FACES", DistributePointsOnFaces, "Distribute Points on Faces", "")
diff --git a/source/blender/nodes/geometry/CMakeLists.txt b/source/blender/nodes/geometry/CMakeLists.txt
index 950124f75d0..d87f0312958 100644
--- a/source/blender/nodes/geometry/CMakeLists.txt
+++ b/source/blender/nodes/geometry/CMakeLists.txt
@@ -57,6 +57,7 @@ set(SRC
   nodes/node_geo_curve_to_mesh.cc
   nodes/node_geo_curve_to_points.cc
   nodes/node_geo_curve_trim.cc
+  nodes/node_geo_deform_curves_on_surface.cc
   nodes/node_geo_delete_geometry.cc
   nodes/node_geo_distribute_points_on_faces.cc
   nodes/node_geo_dual_mesh.cc
diff --git a/source/blender/nodes/geometry/nodes/node_geo_deform_curves_on_surface.cc b/source/blender/nodes/geometry/nodes/node_geo_deform_curves_on_surface.cc
new file mode 100644
index 00000000000..60b5f0383ca
--- /dev/null
+++ b/source/blender/nodes/geometry/nodes/node_geo_deform_curves_on_surface.cc
@@ -0,0 +1,364 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "BKE_attribute_math.hh"
+#include "BKE_curves.hh"
+#include "BKE_editmesh.h"
+#include "BKE_lib_id.h"
+#include "BKE_mesh.h"
+#include "BKE_mesh_runtime.h"
+#include "BKE_mesh_wrapper.h"
+#include "BKE_modifier.h"
+#include "BKE_type_conversions.hh"
+
+#include "BLI_float3x3.hh"
+#include "BLI_task.hh"
+
+#include "UI_interface.h"
+#include "UI_resources.h"
+
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+
+#include "NOD_socket_search_link.hh"
+
+#include "GEO_reverse_uv_sampler.hh"
+
+#include "DEG_depsgraph_query.h"
+
+#include "node_geometry_util.hh"
+
+namespace blender::nodes::node_geo_deform_curves_on_surface_cc {
+
+using attribute_math::mix3;
+using bke::CurvesGeometry;
+using geometry::ReverseUVSampler;
+
+NODE_STORAGE_FUNCS(NodeGeometryCurveTrim)
+
+static void node_declare(NodeDeclarationBuilder &b)
+{
+  b.add_input<decl::Geometry>(N_("Curves")).supported_type(GEO_COMPONENT_TYPE_CURVE);
+  b.add_output<decl::Geometry>(N_("Curves"));
+}
+
+static void deform_curves(CurvesGeometry &curves,
+                          const Mesh &surface_mesh_old,
+                          const Mesh &surface_mesh_new,
+                          const Span<float2> curve_attachment_uvs,
+                          const ReverseUVSampler &reverse_uv_sampler_old,
+                          const ReverseUVSampler &reverse_uv_sampler_new,
+                          const Span<float3> corner_normals_old,
+                          const Span<float3> corner_normals_new,
+                          const Span<float3> rest_positions,
+                          const float4x4 &surface_to_curves,
+                          std::atomic<int> &r_invalid_uv_count)
+{
+  /* Find attachment points on old and new mesh. */
+  const int curves_num = curves.curves_num();
+  Array<ReverseUVSampler::Result> surface_samples_old(curves_num);
+  Array<ReverseUVSampler::Result> surface_samples_new(curves_num);
+  threading::parallel_invoke(
+      [&]() { reverse_uv_sampler_old.sample_many(curve_attachment_uvs, surface_samples_old); },
+      [&]() { reverse_uv_sampler_new.sample_many(curve_attachment_uvs, surface_samples_new); });
+
+  MutableSpan<float3> positions = curves.positions_for_write();
+
+  const float4x4 curves_to_surface = surface_to_curves.inverted();
+
+  threading::parallel_for(curves.curves_range(), 256, [&](const IndexRange range) {
+    for (const int curve_i : range) {
+      const ReverseUVSampler::Result &surface_sample_old = surface_samples_old[curve_i];
+      if (surface_sample_old.type != ReverseUVSampler::ResultType::Ok) {
+        r_invalid_uv_count++;
+        continue;
+      }
+      const ReverseUVSampler::Result &surface_sample_new = surface_samples_new[curve_i];
+      if (surface_sample_new.type != ReverseUVSampler::ResultType::Ok) {
+        r_invalid_uv_count++;
+        continue;
+      }
+
+      const MLoopTri &looptri_old = *surface_sample_old.looptri;
+      const MLoopTri &looptri_new = *surface_sample_new.looptri;
+      const float3 &bary_weights_old = surface_sample_old.bary_weights;
+      const float3 &bary_weights_new = surface_sample_new.bary_weights;
+
+      const int corner_0_old = looptri_old.tri[0];
+      const int corner_1_old = looptri_old.tri[1];
+      const int corner_2_old = looptri_old.tri[2];
+
+      const int corner_0_new = looptri_new.tri[0];
+      const int corner_1_new = looptri_new.tri[1];
+      const int corner_2_new = looptri_new.tri[2];
+
+      const int vert_0_old = surface_mesh_old.mloop[corner_0_old].v;
+      const int vert_1_old = surface_mesh_old.mloop[corner_1_old].v;
+      const int vert_2_old = surface_mesh_old.mloop[corner_2_old].v;
+
+      const int vert_0_new = surface_mesh_new.mloop[corner_0_new].v;
+      const int vert_1_new = surface_mesh_new.mloop[corner_1_new].v;
+      const int vert_2_new = surface_mesh_new.mloop[corner_2_new].v;
+
+      const float3 &normal_0_old = corner_normals_old[corner_0_old];
+      const float3 &normal_1_old = corner_normals_old[corner_1_old];
+      const float3 &normal_2_old = corner_normals_old[corner_2_old];
+      const float3 normal_old = math::normalize(
+          mix3(bary_weights_old, normal_0_old, normal_1_old, normal_2_old));
+
+      const float3 &normal_0_new = corner_normals_new[corner_0_new];
+      const float3 &normal_1_new = corner_normals_new[corner_1_new];
+      const float3 &normal_2_new = corner_normals_new[corner_2_new];
+      const float3 normal_new = math::normalize(
+          mix3(bary_weights_new, normal_0_new, normal_1_new, normal_2_new));
+
+      const float3 &pos_0_old = surface_mesh_old.mvert[vert_0_old].co;
+      const float3 &pos_1_old = surface_mesh_old.mvert[vert_1_old].co;
+      const float3 &pos_2_old = surface_mesh_old.mvert[vert_2_old].co;
+      const float3 pos_old = mix3(bary_weights_old, pos_0_old, pos_1_old, pos_2_old);
+
+      const float3 &pos_0_new = surface_mesh_new.mvert[vert_0_new].co;
+      const float3 &pos_1_new = surface_mesh_new.mvert[vert_1_new].co;
+      const float3 &pos_2_new = surface_mesh_new.mvert[vert_2_new].co;
+      const float3 pos_new = mix3(bary_weights_new, pos_0_new, pos_1_new, pos_2_new);
+
+      /* The translation is just the difference between the old and new position on the surface. */
+      const float3 translation = pos_new - pos_old;
+
+      const float3 &rest_pos_0 = rest_positions[vert_0_new];
+      const float3 &rest_pos_1 = rest_positions[vert_1_new];
+
+      /* The tangent reference direction is used to determine the rotation of the surface point
+       * around its normal axis. It's important that the old and new tangent reference are computed
+       * in a consistent way. If the surface has not been rotated, the old and new tangent
+       * reference have to have the same direction. For that reason, the old tangent reference is
+       * computed based on the rest position attribute instead of positions on the old mesh. This
+       * way the old and new tangent reference use the same topology.
+       *
+       * TODO: Figure out if this can be smoothly interpolated across the surface as well.
+       * Currently, this is a source of discontinuity in the deformation, because the vector
+       * changes intantly from one triangle to the next. */
+      const float3 tangent_reference_dir_old = rest_pos_1 - rest_pos_0;
+      const float3 tangent_reference_dir_new = pos_1_new - pos_0_new;
+
+      /* Compute first local tangent based on the (potentially smoothed) normal and the tangent
+       * reference. */
+      const float3 tangent_x_old = math::normalize(
+          math::cross(normal_old, tangent_reference_dir_old));
+      const float3 tangent_x_new = math::normalize(
+          math::cross(normal_new, tangent_reference_dir_new));
+
+      /* The second tangent defined by the normal and first tangent. */
+      const float3 tangent_y_old = math::normalize(math::cross(normal_old, tangent_x_old));
+      const float3 tangent_y_new = math::normalize(math::cross(normal_new, tangent_x_new));
+
+      /* Construct rotation matrix that encodes the orientation of the old surface position. */
+      float3x3 rotation_old;
+      copy_v3_v3(rotation_old.values[0], tangent_x_old);
+      copy_v3_v3(rotation_old.values[1], tangent_y_old);
+      copy_v3_v3(rotation_old.values[2], normal_old);
+
+      /* Construct rotation matrix that encodes the orientation of the new surface position. */
+      float3x3 rotation_new;
+      copy_v3_v3(rotation_new.values[0], tangent_x_new);
+      copy_v3_v3(rotation_new.values[1], tangent_y_new);
+      copy_v3_v3(rotation_new.values[2], normal_new);
+
+      /* Can use transpose instead of inverse because the matrix is orthonormal. In the case of
+       * zero-area triangles, the matrix would not be orthonormal, but in this case, none of this
+       * works anyway. */
+      const float3x3 rotation_old_inv = rotation_old.transposed();
+
+      /* Compute a rotation matrix that rotates points from the old to the new surface
+       * orientation. */
+      const float3x3 rotation = rotation_new * rotation_old_inv;
+      float4x4 rotation_4x4;
+      copy_m4_m3(rotation_4x4.values, rotation.values);
+
+      /* Construction transformation matrix for this surface position that includes rotation and
+       * translation. */
+      float4x4 surface_transform = float4x4::identity();
+      /* Subtract and add #pos_old, so that the rotation origin is the position on the surface. */
+      sub_v3_v3(surface_transform.values[3], pos_old);
+      mul_m4_m4_pre(surface_transform.values, rotation_4x4.values);
+      add_v3_v3(surface_transform.values[3], pos_old);
+      add_v3_v3(surface_transform.values[3], translation);
+
+      /* Change the basis of the transformation so to that it can be applied in the local space of
+       * the curves. */
+      const float4x4 curve_transform = surface_to_curves * surface_transform * curves_to_surface;
+
+      /* Actually transform all points. */
+      const IndexRange points = curves.points_for_curve(curve_i);
+      for (const int point_i : points) {
+        const float3 old_point_pos = positions[point_i];
+        const float3 new_point_pos = curve_transform * old_point_pos;
+        positions[point_i] = new_point_pos;
+      }
+    }
+  });
+}
+
+static void node_geo_exec(GeoNodeExecParams params)
+{
+  GeometrySet curves_geometry = params.extract_input<GeometrySet>("Curves");
+
+  Mesh *surface_mesh_orig = nullptr;
+  bool free_suface_mesh_orig = false;
+  BLI_SCOPED_DEFER([&]() {
+    if (free_suface_mesh_orig) {
+      BKE_id_free(nullptr, surface_mesh_orig);
+    }
+  });
+
+  auto pass_through_input = [&]() { params.set_output("Curves", std::move(curves_geometry)); };
+
+  const Object *self_ob_eval = params.self_object();
+  if (self_ob_eval == nullptr || self_ob_eval->type != OB_CURVES) {
+    pass_through_input();
+    return;
+  }
+  const Curves *self_curves_eval = static_cast<const Curves *>(self_ob_eval->data);
+  /* Take surface information from self-object. */
+  Object *surface_ob_eval = self_curves_eval->surface;
+  const StringRefNull uv_map_name = self_curves_eval->surface_uv_map;
+  const StringRefNull rest_position_name = "rest_position";
+
+  if (!curves_geometry.has_curves()) {
+    pass_through_input();
+    return;
+  }
+  if (surface_ob_eval == nullptr || surface_ob_eval->type != OB_MESH) {
+    pass_through_input();
+    params.error_message_add(NodeWarningType::Error, "Curves not attached to a surface.");
+    return;
+  }
+  Object *surface_ob_orig = DEG_get_original_object(surface_ob_eval);
+  Mesh &surface_object_data = *static_cast<Mesh *>(surface_ob_orig->data);
+
+  if (BMEditMesh *em = surface_object_data.edit_mesh) {
+    surface_mesh_orig = BKE_mesh_from_bmesh_for_eval_nomain(em->bm, NULL, &surface_object_data);
+    free_suface_mesh_orig = true;
+  }
+  else {
+    surface_mesh_orig = &surface_object_data;
+  }
+  Mesh *surface_mesh_eval = BKE_modifier_get_evaluated_mesh_from_evaluated_object(surface_ob_eval,
+                                                                                  false);
+  if (surface_mesh_eval == nullptr) {
+    pass_through_input();
+    params.error_message_add(NodeWarningType::Error, "Surface has no mesh.");
+    return;
+  }
+
+  BKE_mesh_wrapper_ensure_mdata(surface_mesh_eval);
+
+  MeshComponent mesh_eval;
+  mesh_eval.replace(surface_mesh_eval, GeometryOwnershipType::ReadOnly);
+  MeshComponent mesh_orig;
+  mesh_orig.replace(surface_mesh_orig, GeometryOwnershipType::ReadOnly);
+
+  Curves &curves_id = *curves_geometry.get_curves_for_write();
+  CurvesGeometry &curves = CurvesGeometry::wrap(curves_id.geometry);
+
+  if (uv_map_name.is_empty()) {
+    pass_through_input();
+    const char *message = TIP_("Surface UV map not defined.");
+    params.error_message_add(NodeWarningType::Error, message);
+    return;
+  }
+  if (!mesh_eval.attribute_exists(uv_map_name)) {
+    pass_through_input();
+    char *message = BLI_sprintfN(TIP_("Evaluated surface missing UV map: %s."),
+                                 uv_map_name.c_str());
+    params.error_message_add(NodeWarningType::Error, message);
+    MEM_freeN(message);
+    return;
+  }
+  if (!mesh_orig.attribute_exists(uv_map_name)) {
+    pass_through_input();
+    char *message = BLI_sprintfN(TIP_("Original surface missing UV map: %s."),
+                                 uv_map_name.c_str());
+    params.error_message_add(NodeWarningType::Error, message);
+    MEM_freeN(message);
+    return;
+  }
+  if (!mesh_eval.attribute_exists(rest_position_name)) {
+    pass_through_input();
+    params.error_message_add(NodeWarningType::Error,
+                             TIP_("Evaluated surface missing attribute: rest_position."));
+    return;
+  }
+  if (curves.surface_uv_coords().is_empty()) {
+    pass_through_input();
+    params.error_message_add(NodeWarningType::Error,
+                             TIP_("Curves are not attached to any UV map."));
+    return;
+  }
+  const VArraySpan<float2> uv_map_orig = mesh_orig.attribute_get_for_read<float2>(
+      uv_map_name, ATTR_DOMAIN_CORNER, {0.0f, 0.0f});
+  const VArraySpan<float2> uv_map_eval = mesh_eval.attribute_get_for_read<float2>(
+      uv_map_name, ATTR_DOMAIN_CORNER, {0.0f, 0.0f});
+  const VArraySpan<float3> rest_positions = mesh_eval.attribute_get_for_read<float3>(
+      rest_position_name, ATTR_DOMAIN_POINT, {0.0f, 0.0f, 0.0f});
+  const Span<float2> surface_uv_coords = curves.surface_uv_coords();
+
+  const Span<MLoopTri> looptris_orig{BKE_mesh_runtime_looptri_ensure(surface_mesh_orig),
+                                     BKE_mesh_runtime_looptri_len(surface_mesh_orig)};
+  const Span<MLoopTri> looptris_eval{BKE_mesh_runtime_looptri_ensure(surface_mesh_eval),
+                                     BKE_mesh_runtime_looptri_len(surface_mesh_eval)};
+  const ReverseUVSampler reverse_uv_sampler_orig{uv_map_orig, looptris_orig};
+  const ReverseUVSampler reverse_uv_sampler_eval{uv_map_eval, looptris_eval};
+
+  /* Retrieve face corner normals from each mesh. It's necessary to use face corner normals
+   * because face normals or vertex normals may lose information (custom normals, auto smooth) in
+   * some cases. It isn't yet possible to retrieve lazily calculated face corner normals from a
+   * const mesh, so they are calculated here every time. */
+  Array<float3> corner_normals_orig(surface_mesh_orig->totloop);
+  Array<float3> corner_normals_eval(surface_mesh_eval->totloop);
+  BKE_mesh_calc_normals_split_ex(
+      surface_mesh_orig, nullptr, reinterpret_cast<float(*)[3]>(corner_normals_orig.data()));
+  BKE_mesh_calc_normals_split_ex(
+      surface_mesh_eval, nullptr, reinterpret_cast<float(*)[3]>(corner_normals_eval.data()));
+
+  std::atomic<int> invalid_uv_count = 0;
+
+  const bke::CurvesSurfaceTransforms transforms{*self_ob_eval, surface_ob_eval};
+
+  deform_curves(curves,
+                *surface_mesh_orig,
+                *surface_mesh_eval,
+                surface_uv_coords,
+                reverse_uv_sampler_orig,
+                reverse_uv_sampler_eval,
+                corner_normals_orig,
+                corner_normals_eval,
+                rest_positions,
+                transforms.surface_to_curves,
+                invalid_uv_count);
+
+  curves.tag_positions_changed();
+
+  if (invalid_uv_count) {
+    char *message = BLI_sprintfN(TIP_("Invalid surface UVs on %d curves."),
+                                 invalid_uv_count.load());
+    params.error_message_add(NodeWarningType::Warning, message);
+    MEM_freeN(message);
+  }
+
+  params.set_output("Curves", curves_geometry);
+}
+
+}  // namespace blender::nodes::node_geo_deform_curves_on_surface_cc
+
+void register_node_type_geo_deform_curves_on_surface()
+{
+  namespace file_ns = blender::nodes::node_geo_deform_curves_on_surface_cc;
+
+  static bNodeType ntype;
+  geo_node_type_base(
+      &ntype, GEO_NODE_DEFORM_CURVES_ON_SURFACE, "Deform Curves on Surface", NODE_CLASS_GEOMETRY);
+  ntype.geometry_node_execute = file_ns::node_geo_exec;
+  ntype.declare = file_ns::node_declare;
+  node_type_size(&ntype, 170, 120, 700);
+  nodeRegisterType(&ntype);
+}