Curves: initial surface collision for curves sculpt mode

During hair grooming in curves sculpt mode, it is very useful when hair strands
are prevented from intersecting with the surface mesh. Unfortunately, it also
decreases performance significantly so we don't want it to be turned on all the time.

The surface collision is used by the Comb, Pinch and Puff brushes currently.
It can be turned on or off on a per-geometry basis.

The intersection prevention quality of this patch is not perfect yet. This can
be improved over time using a better solver. Overall, perfect collision detection
at the cost of bad performance is not necessary for interactive sculpting,
because the user can fix small mistakes very quickly. Nevertheless, the quality
can probably still be improved significantly without too big slow-downs depending
on the use case. This can be done separately from this patch.

Pull Request #104469

release/scripts/startup/bl_ui/space_view3d.py

@@ -161,6 +161,8 @@ class VIEW3D_HT_tool_header(Header):
             sub.prop(context.object.data, "use_mirror_y", text="Y", toggle=True)
             sub.prop(context.object.data, "use_mirror_z", text="Z", toggle=True)
 
+            layout.prop(context.object.data, "use_sculpt_collision", icon='MOD_PHYSICS', icon_only=True, toggle=True)
+
         # Expand panels from the side-bar as popovers.
         popover_kw = {"space_type": 'VIEW_3D', "region_type": 'UI', "category": "Tool"}
 

source/blender/blenlib/BLI_index_mask_ops.hh

@@ -71,4 +71,9 @@ IndexMask find_indices_from_virtual_array(IndexMask indices_to_check,
                                           int64_t parallel_grain_size,
                                           Vector<int64_t> &r_indices);
 
+/**
+ * Find the true indices in a boolean span.
+ */
+IndexMask find_indices_from_array(Span<bool> array, Vector<int64_t> &r_indices);
+
 }  // namespace blender::index_mask_ops

source/blender/blenlib/intern/index_mask.cc

@@ -208,8 +208,7 @@ IndexMask find_indices_from_virtual_array(const IndexMask indices_to_check,
   }
   if (virtual_array.is_span()) {
     const Span<bool> span = virtual_array.get_internal_span();
-    return find_indices_based_on_predicate(
-        indices_to_check, 4096, r_indices, [&](const int64_t i) { return span[i]; });
+    return find_indices_from_array(span, r_indices);
   }
 
   threading::EnumerableThreadSpecific<Vector<bool>> materialize_buffers;
@@ -241,4 +240,10 @@ IndexMask find_indices_from_virtual_array(const IndexMask indices_to_check,
   return detail::find_indices_based_on_predicate__merge(indices_to_check, sub_masks, r_indices);
 }
 
+IndexMask find_indices_from_array(const Span<bool> array, Vector<int64_t> &r_indices)
+{
+  return find_indices_based_on_predicate(
+      array.index_range(), 4096, r_indices, [array](const int64_t i) { return array[i]; });
+}
+
 }  // namespace blender::index_mask_ops

source/blender/editors/sculpt_paint/curves_sculpt_brush.cc

@@ -22,6 +22,8 @@
 
 #include "BLT_translation.h"
 
+#include "GEO_curve_constraints.hh"
+
 /**
  * The code below uses a prefix naming convention to indicate the coordinate space:
  * cu: Local space of the curves object that is being edited.
@@ -431,4 +433,40 @@ void report_invalid_uv_map(ReportList *reports)
   BKE_report(reports, RPT_WARNING, TIP_("Invalid UV map: UV islands must not overlap"));
 }
 
+void CurvesConstraintSolver::initialize(const bke::CurvesGeometry &curves,
+                                        const IndexMask curve_selection,
+                                        const bool use_surface_collision)
+{
+  use_surface_collision_ = use_surface_collision;
+  segment_lengths_.reinitialize(curves.points_num());
+  geometry::curve_constraints::compute_segment_lengths(
+      curves.points_by_curve(), curves.positions(), curve_selection, segment_lengths_);
+  if (use_surface_collision_) {
+    start_positions_ = curves.positions();
+  }
+}
+
+void CurvesConstraintSolver::solve_step(bke::CurvesGeometry &curves,
+                                        const IndexMask curve_selection,
+                                        const Mesh *surface,
+                                        const CurvesSurfaceTransforms &transforms)
+{
+  if (use_surface_collision_ && surface != nullptr) {
+    geometry::curve_constraints::solve_length_and_collision_constraints(
+        curves.points_by_curve(),
+        curve_selection,
+        segment_lengths_,
+        start_positions_,
+        *surface,
+        transforms,
+        curves.positions_for_write());
+    start_positions_ = curves.positions();
+  }
+  else {
+    geometry::curve_constraints::solve_length_constraints(
+        curves.points_by_curve(), curve_selection, segment_lengths_, curves.positions_for_write());
+  }
+  curves.tag_positions_changed();
+}
+
 }  // namespace blender::ed::sculpt_paint

source/blender/editors/sculpt_paint/curves_sculpt_comb.cc

@@ -66,8 +66,8 @@ class CombOperation : public CurvesSculptStrokeOperation {
   /** Only used when a 3D brush is used. */
   CurvesBrush3D brush_3d_;
 
-  /** Length of each segment indexed by the index of the first point in the segment. */
-  Array<float> segment_lengths_cu_;
+  /** Solver for length and collision constraints. */
+  CurvesConstraintSolver constraint_solver_;
 
   friend struct CombOperationExecutor;
 
@@ -144,12 +144,13 @@ struct CombOperationExecutor {
       if (falloff_shape_ == PAINT_FALLOFF_SHAPE_SPHERE) {
         this->initialize_spherical_brush_reference_point();
       }
-      this->initialize_segment_lengths();
+      self_->constraint_solver_.initialize(
+          *curves_orig_, curve_selection_, curves_id_orig_->flag & CV_SCULPT_COLLISION_ENABLED);
       /* Combing does nothing when there is no mouse movement, so return directly. */
       return;
     }
 
-    EnumerableThreadSpecific<Vector<int>> changed_curves;
+    Array<bool> changed_curves(curves_orig_->curves_num(), false);
 
     if (falloff_shape_ == PAINT_FALLOFF_SHAPE_TUBE) {
       this->comb_projected_with_symmetry(changed_curves);
@@ -161,7 +162,14 @@ struct CombOperationExecutor {
       BLI_assert_unreachable();
     }
 
-    this->restore_segment_lengths(changed_curves);
+    const Mesh *surface = curves_id_orig_->surface && curves_id_orig_->surface->type == OB_MESH ?
+                              static_cast<Mesh *>(curves_id_orig_->surface->data) :
+                              nullptr;
+
+    Vector<int64_t> indices;
+    const IndexMask changed_curves_mask = index_mask_ops::find_indices_from_array(changed_curves,
+                                                                                  indices);
+    self_->constraint_solver_.solve_step(*curves_orig_, changed_curves_mask, surface, transforms_);
 
     curves_orig_->tag_positions_changed();
     DEG_id_tag_update(&curves_id_orig_->id, ID_RECALC_GEOMETRY);
@@ -172,7 +180,7 @@ struct CombOperationExecutor {
   /**
    * Do combing in screen space.
    */
-  void comb_projected_with_symmetry(EnumerableThreadSpecific<Vector<int>> &r_changed_curves)
+  void comb_projected_with_symmetry(MutableSpan<bool> r_changed_curves)
   {
     const Vector<float4x4> symmetry_brush_transforms = get_symmetry_brush_transforms(
         eCurvesSymmetryType(curves_id_orig_->symmetry));
@@ -181,8 +189,7 @@ struct CombOperationExecutor {
     }
   }
 
-  void comb_projected(EnumerableThreadSpecific<Vector<int>> &r_changed_curves,
-                      const float4x4 &brush_transform)
+  void comb_projected(MutableSpan<bool> r_changed_curves, const float4x4 &brush_transform)
   {
     const float4x4 brush_transform_inv = math::invert(brush_transform);
 
@@ -198,7 +205,6 @@ struct CombOperationExecutor {
     const float brush_radius_sq_re = pow2f(brush_radius_re);
 
     threading::parallel_for(curve_selection_.index_range(), 256, [&](const IndexRange range) {
-      Vector<int> &local_changed_curves = r_changed_curves.local();
       for (const int curve_i : curve_selection_.slice(range)) {
         bool curve_changed = false;
         const IndexRange points = points_by_curve[curve_i];
@@ -246,7 +252,7 @@ struct CombOperationExecutor {
           curve_changed = true;
         }
         if (curve_changed) {
-          local_changed_curves.append(curve_i);
+          r_changed_curves[curve_i] = true;
         }
       }
     });
@@ -255,7 +261,7 @@ struct CombOperationExecutor {
   /**
    * Do combing in 3D space.
    */
-  void comb_spherical_with_symmetry(EnumerableThreadSpecific<Vector<int>> &r_changed_curves)
+  void comb_spherical_with_symmetry(MutableSpan<bool> r_changed_curves)
   {
     float4x4 projection;
     ED_view3d_ob_project_mat_get(ctx_.rv3d, curves_ob_orig_, projection.ptr());
@@ -289,7 +295,7 @@ struct CombOperationExecutor {
     }
   }
 
-  void comb_spherical(EnumerableThreadSpecific<Vector<int>> &r_changed_curves,
+  void comb_spherical(MutableSpan<bool> r_changed_curves,
                       const float3 &brush_start_cu,
                       const float3 &brush_end_cu,
                       const float brush_radius_cu)
@@ -303,7 +309,6 @@ struct CombOperationExecutor {
     const OffsetIndices points_by_curve = curves_orig_->points_by_curve();
 
     threading::parallel_for(curve_selection_.index_range(), 256, [&](const IndexRange range) {
-      Vector<int> &local_changed_curves = r_changed_curves.local();
       for (const int curve_i : curve_selection_.slice(range)) {
         bool curve_changed = false;
         const IndexRange points = points_by_curve[curve_i];
@@ -335,7 +340,7 @@ struct CombOperationExecutor {
           curve_changed = true;
         }
         if (curve_changed) {
-          local_changed_curves.append(curve_i);
+          r_changed_curves[curve_i] = true;
         }
       }
     });
@@ -357,53 +362,6 @@ struct CombOperationExecutor {
       self_->brush_3d_ = *brush_3d;
     }
   }
-
-  /**
-   * Remember the initial length of all curve segments. This allows restoring the length after
-   * combing.
-   */
-  void initialize_segment_lengths()
-  {
-    const Span<float3> positions_cu = curves_orig_->positions();
-    const OffsetIndices points_by_curve = curves_orig_->points_by_curve();
-    self_->segment_lengths_cu_.reinitialize(curves_orig_->points_num());
-    threading::parallel_for(curves_orig_->curves_range(), 128, [&](const IndexRange range) {
-      for (const int curve_i : range) {
-        const IndexRange points = points_by_curve[curve_i];
-        for (const int point_i : points.drop_back(1)) {
-          const float3 &p1_cu = positions_cu[point_i];
-          const float3 &p2_cu = positions_cu[point_i + 1];
-          const float length_cu = math::distance(p1_cu, p2_cu);
-          self_->segment_lengths_cu_[point_i] = length_cu;
-        }
-      }
-    });
-  }
-
-  /**
-   * Restore previously stored length for each segment in the changed curves.
-   */
-  void restore_segment_lengths(EnumerableThreadSpecific<Vector<int>> &changed_curves)
-  {
-    const Span<float> expected_lengths_cu = self_->segment_lengths_cu_;
-    const OffsetIndices points_by_curve = curves_orig_->points_by_curve();
-    MutableSpan<float3> positions_cu = curves_orig_->positions_for_write();
-
-    threading::parallel_for_each(changed_curves, [&](const Vector<int> &changed_curves) {
-      threading::parallel_for(changed_curves.index_range(), 256, [&](const IndexRange range) {
-        for (const int curve_i : changed_curves.as_span().slice(range)) {
-          const IndexRange points = points_by_curve[curve_i];
-          for (const int segment_i : points.drop_back(1)) {
-            const float3 &p1_cu = positions_cu[segment_i];
-            float3 &p2_cu = positions_cu[segment_i + 1];
-            const float3 direction = math::normalize(p2_cu - p1_cu);
-            const float expected_length_cu = expected_lengths_cu[segment_i];
-            p2_cu = p1_cu + direction * expected_length_cu;
-          }
-        }
-      });
-    });
-  }
 };
 
 void CombOperation::on_stroke_extended(const bContext &C, const StrokeExtension &stroke_extension)

source/blender/editors/sculpt_paint/curves_sculpt_intern.hh

@@ -144,4 +144,25 @@ void report_missing_uv_map_on_original_surface(ReportList *reports);
 void report_missing_uv_map_on_evaluated_surface(ReportList *reports);
 void report_invalid_uv_map(ReportList *reports);
 
+/**
+ * Utility class to make it easy for brushes to implement length preservation and surface
+ * collision.
+ */
+struct CurvesConstraintSolver {
+ private:
+  bool use_surface_collision_;
+  Array<float3> start_positions_;
+  Array<float> segment_lengths_;
+
+ public:
+  void initialize(const bke::CurvesGeometry &curves,
+                  const IndexMask curve_selection,
+                  const bool use_surface_collision);
+
+  void solve_step(bke::CurvesGeometry &curves,
+                  const IndexMask curve_selection,
+                  const Mesh *surface,
+                  const CurvesSurfaceTransforms &transforms);
+};
+
 }  // namespace blender::ed::sculpt_paint

source/blender/editors/sculpt_paint/curves_sculpt_pinch.cc

@@ -4,6 +4,7 @@
 
 #include "curves_sculpt_intern.hh"
 
+#include "BLI_index_mask_ops.hh"
 #include "BLI_math_matrix_types.hh"
 #include "BLI_task.hh"
 #include "BLI_vector.hh"
@@ -42,7 +43,9 @@ namespace blender::ed::sculpt_paint {
 class PinchOperation : public CurvesSculptStrokeOperation {
  private:
   bool invert_pinch_;
-  Array<float> segment_lengths_cu_;
+
+  /** Solver for length and collision constraints. */
+  CurvesConstraintSolver constraint_solver_;
 
   /** Only used when a 3D brush is used. */
   CurvesBrush3D brush_3d_;
@@ -115,8 +118,6 @@ struct PinchOperationExecutor {
         brush_->falloff_shape);
 
     if (stroke_extension.is_first) {
-      this->initialize_segment_lengths();
-
       if (falloff_shape == PAINT_FALLOFF_SHAPE_SPHERE) {
         self_->brush_3d_ = *sample_curves_3d_brush(*ctx_.depsgraph,
                                                    *ctx_.region,
@@ -126,6 +127,9 @@ struct PinchOperationExecutor {
                                                    brush_pos_re_,
                                                    brush_radius_base_re_);
       }
+
+      self_->constraint_solver_.initialize(
+          *curves_, curve_selection_, curves_id_->flag & CV_SCULPT_COLLISION_ENABLED);
     }
 
     Array<bool> changed_curves(curves_->curves_num(), false);
@@ -139,7 +143,14 @@ struct PinchOperationExecutor {
       BLI_assert_unreachable();
     }
 
-    this->restore_segment_lengths(changed_curves);
+    Vector<int64_t> indices;
+    const IndexMask changed_curves_mask = index_mask_ops::find_indices_from_array(changed_curves,
+                                                                                  indices);
+    const Mesh *surface = curves_id_->surface && curves_id_->surface->type == OB_MESH ?
+                              static_cast<const Mesh *>(curves_id_->surface->data) :
+                              nullptr;
+    self_->constraint_solver_.solve_step(*curves_, changed_curves_mask, surface, transforms_);
+
     curves_->tag_positions_changed();
     DEG_id_tag_update(&curves_id_->id, ID_RECALC_GEOMETRY);
     WM_main_add_notifier(NC_GEOM | ND_DATA, &curves_id_->id);
@@ -270,47 +281,6 @@ struct PinchOperationExecutor {
       }
     });
   }
-
-  void initialize_segment_lengths()
-  {
-    const Span<float3> positions_cu = curves_->positions();
-    const OffsetIndices points_by_curve = curves_->points_by_curve();
-    self_->segment_lengths_cu_.reinitialize(curves_->points_num());
-    threading::parallel_for(curve_selection_.index_range(), 256, [&](const IndexRange range) {
-      for (const int curve_i : curve_selection_.slice(range)) {
-        const IndexRange points = points_by_curve[curve_i];
-        for (const int point_i : points.drop_back(1)) {
-          const float3 &p1_cu = positions_cu[point_i];
-          const float3 &p2_cu = positions_cu[point_i + 1];
-          const float length_cu = math::distance(p1_cu, p2_cu);
-          self_->segment_lengths_cu_[point_i] = length_cu;
-        }
-      }
-    });
-  }
-
-  void restore_segment_lengths(const Span<bool> changed_curves)
-  {
-    const Span<float> expected_lengths_cu = self_->segment_lengths_cu_;
-    const OffsetIndices points_by_curve = curves_->points_by_curve();
-    MutableSpan<float3> positions_cu = curves_->positions_for_write();
-
-    threading::parallel_for(changed_curves.index_range(), 256, [&](const IndexRange range) {
-      for (const int curve_i : range) {
-        if (!changed_curves[curve_i]) {
-          continue;
-        }
-        const IndexRange points = points_by_curve[curve_i];
-        for (const int segment_i : IndexRange(points.size() - 1)) {
-          const float3 &p1_cu = positions_cu[points[segment_i]];
-          float3 &p2_cu = positions_cu[points[segment_i] + 1];
-          const float3 direction = math::normalize(p2_cu - p1_cu);
-          const float expected_length_cu = expected_lengths_cu[points[segment_i]];
-          p2_cu = p1_cu + direction * expected_length_cu;
-        }
-      }
-    });
-  }
 };
 
 void PinchOperation::on_stroke_extended(const bContext &C, const StrokeExtension &stroke_extension)

source/blender/editors/sculpt_paint/curves_sculpt_puff.cc

@@ -19,6 +19,7 @@
 
 #include "WM_api.h"
 
+#include "BLI_index_mask_ops.hh"
 #include "BLI_length_parameterize.hh"
 #include "BLI_math_matrix.hh"
 #include "BLI_task.hh"
@@ -34,8 +35,8 @@ class PuffOperation : public CurvesSculptStrokeOperation {
   /** Only used when a 3D brush is used. */
   CurvesBrush3D brush_3d_;
 
-  /** Length of each segment indexed by the index of the first point in the segment. */
-  Array<float> segment_lengths_cu_;
+  /** Solver for length and collision constraints. */
+  CurvesConstraintSolver constraint_solver_;
 
   friend struct PuffOperationExecutor;
 
@@ -130,7 +131,6 @@ struct PuffOperationExecutor {
     BLI_SCOPED_DEFER([&]() { free_bvhtree_from_mesh(&surface_bvh_); });
 
     if (stroke_extension.is_first) {
-      this->initialize_segment_lengths();
       if (falloff_shape_ == PAINT_FALLOFF_SHAPE_SPHERE) {
         self.brush_3d_ = *sample_curves_3d_brush(*ctx_.depsgraph,
                                                  *ctx_.region,
@@ -140,6 +140,9 @@ struct PuffOperationExecutor {
                                                  brush_pos_re_,
                                                  brush_radius_base_re_);
       }
+
+      self_->constraint_solver_.initialize(
+          *curves_, curve_selection_, curves_id_->flag & CV_SCULPT_COLLISION_ENABLED);
     }
 
     Array<float> curve_weights(curve_selection_.size(), 0.0f);
@@ -155,7 +158,17 @@ struct PuffOperationExecutor {
     }
 
     this->puff(curve_weights);
-    this->restore_segment_lengths();
+
+    Vector<int64_t> changed_curves_indices;
+    changed_curves_indices.reserve(curve_selection_.size());
+    for (int64_t select_i : curve_selection_.index_range()) {
+      if (curve_weights[select_i] > 0.0f) {
+        changed_curves_indices.append(curve_selection_[select_i]);
+      }
+    }
+
+    self_->constraint_solver_.solve_step(
+        *curves_, IndexMask(changed_curves_indices), surface_, transforms_);
 
     curves_->tag_positions_changed();
     DEG_id_tag_update(&curves_id_->id, ID_RECALC_GEOMETRY);
@@ -344,44 +357,6 @@ struct PuffOperationExecutor {
       }
     });
   }
-
-  void initialize_segment_lengths()
-  {
-    const OffsetIndices points_by_curve = curves_->points_by_curve();
-    const Span<float3> positions_cu = curves_->positions();
-    self_->segment_lengths_cu_.reinitialize(curves_->points_num());
-    threading::parallel_for(curves_->curves_range(), 128, [&](const IndexRange range) {
-      for (const int curve_i : range) {
-        const IndexRange points = points_by_curve[curve_i];
-        for (const int point_i : points.drop_back(1)) {
-          const float3 &p1_cu = positions_cu[point_i];
-          const float3 &p2_cu = positions_cu[point_i + 1];
-          const float length_cu = math::distance(p1_cu, p2_cu);
-          self_->segment_lengths_cu_[point_i] = length_cu;
-        }
-      }
-    });
-  }
-
-  void restore_segment_lengths()
-  {
-    const Span<float> expected_lengths_cu = self_->segment_lengths_cu_;
-    const OffsetIndices points_by_curve = curves_->points_by_curve();
-    MutableSpan<float3> positions_cu = curves_->positions_for_write();
-
-    threading::parallel_for(curves_->curves_range(), 256, [&](const IndexRange range) {
-      for (const int curve_i : range) {
-        const IndexRange points = points_by_curve[curve_i];
-        for (const int segment_i : points.drop_back(1)) {
-          const float3 &p1_cu = positions_cu[segment_i];
-          float3 &p2_cu = positions_cu[segment_i + 1];
-          const float3 direction = math::normalize(p2_cu - p1_cu);
-          const float expected_length_cu = expected_lengths_cu[segment_i];
-          p2_cu = p1_cu + direction * expected_length_cu;
-        }
-      }
-    });
-  }
 };
 
 void PuffOperation::on_stroke_extended(const bContext &C, const StrokeExtension &stroke_extension)

source/blender/geometry/CMakeLists.txt

@@ -16,6 +16,7 @@ set(INC
 
 set(SRC
   intern/add_curves_on_mesh.cc
+  intern/curve_constraints.cc
   intern/fillet_curves.cc
   intern/mesh_merge_by_distance.cc
   intern/mesh_primitive_cuboid.cc
@@ -32,6 +33,7 @@ set(SRC
   intern/uv_parametrizer.cc
 
   GEO_add_curves_on_mesh.hh
+  GEO_curve_constraints.hh
   GEO_fillet_curves.hh
   GEO_mesh_merge_by_distance.hh
   GEO_mesh_primitive_cuboid.hh

source/blender/geometry/GEO_curve_constraints.hh

@@ -0,0 +1,27 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#pragma once
+
+#include "BKE_curves.hh"
+
+namespace blender::geometry::curve_constraints {
+
+void compute_segment_lengths(OffsetIndices<int> points_by_curve,
+                             Span<float3> positions,
+                             IndexMask curve_selection,
+                             MutableSpan<float> r_segment_lengths);
+
+void solve_length_constraints(OffsetIndices<int> points_by_curve,
+                              IndexMask curve_selection,
+                              Span<float> segment_lenghts,
+                              MutableSpan<float3> positions);
+
+void solve_length_and_collision_constraints(OffsetIndices<int> points_by_curve,
+                                            IndexMask curve_selection,
+                                            Span<float> segment_lengths,
+                                            Span<float3> start_positions,
+                                            const Mesh &surface,
+                                            const bke::CurvesSurfaceTransforms &transforms,
+                                            MutableSpan<float3> positions);
+
+}  // namespace blender::geometry::curve_constraints

source/blender/geometry/intern/curve_constraints.cc

@@ -0,0 +1,168 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "BLI_math_matrix.hh"
+#include "BLI_task.hh"
+
+#include "GEO_curve_constraints.hh"
+
+#include "BKE_bvhutils.h"
+
+/**
+ * The code below uses a prefix naming convention to indicate the coordinate space:
+ * `cu`: Local space of the curves object that is being edited.
+ * `su`: Local space of the surface object.
+ * `wo`: World space.
+ */
+
+namespace blender::geometry::curve_constraints {
+
+void compute_segment_lengths(const OffsetIndices<int> points_by_curve,
+                             const Span<float3> positions,
+                             const IndexMask curve_selection,
+                             MutableSpan<float> r_segment_lengths)
+{
+  BLI_assert(r_segment_lengths.size() == points_by_curve.total_size());
+
+  threading::parallel_for(curve_selection.index_range(), 256, [&](const IndexRange range) {
+    for (const int curve_i : curve_selection.slice(range)) {
+      const IndexRange points = points_by_curve[curve_i].drop_back(1);
+      for (const int point_i : points) {
+        const float3 &p1 = positions[point_i];
+        const float3 &p2 = positions[point_i + 1];
+        const float length = math::distance(p1, p2);
+        r_segment_lengths[point_i] = length;
+      }
+    }
+  });
+}
+
+void solve_length_constraints(const OffsetIndices<int> points_by_curve,
+                              const IndexMask curve_selection,
+                              const Span<float> segment_lenghts,
+                              MutableSpan<float3> positions)
+{
+  BLI_assert(segment_lenghts.size() == points_by_curve.total_size());
+
+  threading::parallel_for(curve_selection.index_range(), 256, [&](const IndexRange range) {
+    for (const int curve_i : curve_selection.slice(range)) {
+      const IndexRange points = points_by_curve[curve_i].drop_back(1);
+      for (const int point_i : points) {
+        const float3 &p1 = positions[point_i];
+        float3 &p2 = positions[point_i + 1];
+        const float3 direction = math::normalize(p2 - p1);
+        const float goal_length = segment_lenghts[point_i];
+        p2 = p1 + direction * goal_length;
+      }
+    }
+  });
+}
+
+void solve_length_and_collision_constraints(const OffsetIndices<int> points_by_curve,
+                                            const IndexMask curve_selection,
+                                            const Span<float> segment_lengths_cu,
+                                            const Span<float3> start_positions_cu,
+                                            const Mesh &surface,
+                                            const bke::CurvesSurfaceTransforms &transforms,
+                                            MutableSpan<float3> positions_cu)
+{
+  solve_length_constraints(points_by_curve, curve_selection, segment_lengths_cu, positions_cu);
+
+  BVHTreeFromMesh surface_bvh;
+  BKE_bvhtree_from_mesh_get(&surface_bvh, &surface, BVHTREE_FROM_LOOPTRI, 2);
+  BLI_SCOPED_DEFER([&]() { free_bvhtree_from_mesh(&surface_bvh); });
+
+  const float radius = 0.001f;
+  const int max_collisions = 5;
+
+  threading::parallel_for(curve_selection.index_range(), 64, [&](const IndexRange range) {
+    for (const int curve_i : curve_selection.slice(range)) {
+      const IndexRange points = points_by_curve[curve_i];
+
+      /* Sometimes not all collisions can be handled. This happens relatively rarely, but if it
+       * happens it's better to just not to move the curve instead of going into the surface. */
+      bool revert_curve = false;
+      for (const int point_i : points.drop_front(1)) {
+        const float goal_segment_length_cu = segment_lengths_cu[point_i - 1];
+        const float3 &prev_pos_cu = positions_cu[point_i - 1];
+        const float3 &start_pos_cu = start_positions_cu[point_i];
+
+        int used_iterations = 0;
+        for ([[maybe_unused]] const int iteration : IndexRange(max_collisions)) {
+          used_iterations++;
+          const float3 &old_pos_cu = positions_cu[point_i];
+          if (start_pos_cu == old_pos_cu) {
+            /* The point did not move, done. */
+            break;
+          }
+
+          /* Check if the point moved through a surface. */
+          const float3 start_pos_su = math::transform_point(transforms.curves_to_surface,
+                                                            start_pos_cu);
+          const float3 old_pos_su = math::transform_point(transforms.curves_to_surface,
+                                                          old_pos_cu);
+          const float3 pos_diff_su = old_pos_su - start_pos_su;
+          float max_ray_length_su;
+          const float3 ray_direction_su = math::normalize_and_get_length(pos_diff_su,
+                                                                         max_ray_length_su);
+          BVHTreeRayHit hit;
+          hit.index = -1;
+          hit.dist = max_ray_length_su + radius;
+          BLI_bvhtree_ray_cast(surface_bvh.tree,
+                               start_pos_su,
+                               ray_direction_su,
+                               radius,
+                               &hit,
+                               surface_bvh.raycast_callback,
+                               &surface_bvh);
+          if (hit.index == -1) {
+            break;
+          }
+          const float3 hit_pos_su = hit.co;
+          const float3 hit_normal_su = hit.no;
+          if (math::dot(hit_normal_su, ray_direction_su) > 0.0f) {
+            /* Moving from the inside to the outside is ok. */
+            break;
+          }
+
+          /* The point was moved through a surface. Now put it back on the correct side of the
+           * surface and slide it on the surface to keep the length the same. */
+
+          const float3 hit_pos_cu = math::transform_point(transforms.surface_to_curves,
+                                                          hit_pos_su);
+          const float3 hit_normal_cu = math::normalize(
+              math::transform_direction(transforms.surface_to_curves_normal, hit_normal_su));
+
+          /* Slide on a plane that is slightly above the surface. */
+          const float3 plane_pos_cu = hit_pos_cu + hit_normal_cu * radius;
+          const float3 plane_normal_cu = hit_normal_cu;
+
+          /* Decompose the current segment into the part normal and tangent to the collision
+           * surface. */
+          const float3 collided_segment_cu = plane_pos_cu - prev_pos_cu;
+          const float3 slide_normal_cu = plane_normal_cu *
+                                         math::dot(collided_segment_cu, plane_normal_cu);
+          const float3 slide_direction_cu = collided_segment_cu - slide_normal_cu;
+
+          float slide_direction_length_cu;
+          const float3 normalized_slide_direction_cu = math::normalize_and_get_length(
+              slide_direction_cu, slide_direction_length_cu);
+
+          /* Use pythagorian theorem to determine how far to slide. */
+          const float slide_distance_cu = std::sqrt(pow2f(goal_segment_length_cu) -
+                                                    math::length_squared(slide_normal_cu)) -
+                                          slide_direction_length_cu;
+          positions_cu[point_i] = plane_pos_cu + normalized_slide_direction_cu * slide_distance_cu;
+        }
+        if (used_iterations == max_collisions) {
+          revert_curve = true;
+          break;
+        }
+      }
+      if (revert_curve) {
+        positions_cu.slice(points).copy_from(start_positions_cu.slice(points));
+      }
+    }
+  });
+}
+
+}  // namespace blender::geometry::curve_constraints

source/blender/makesdna/DNA_curves_types.h

@@ -197,6 +197,7 @@ typedef struct Curves {
 /** #Curves.flag */
 enum {
   HA_DS_EXPAND = (1 << 0),
+  CV_SCULPT_COLLISION_ENABLED = (1 << 1),
 };
 
 /** #Curves.symmetry */

source/blender/makesrna/intern/rna_curves.c

@@ -453,6 +453,13 @@ static void rna_def_curves(BlenderRNA *brna)
   RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
   RNA_def_property_update(prop, 0, "rna_Curves_update_data");
 
+  prop = RNA_def_property(srna, "use_sculpt_collision", PROP_BOOLEAN, PROP_NONE);
+  RNA_def_property_boolean_sdna(prop, NULL, "flag", CV_SCULPT_COLLISION_ENABLED);
+  RNA_def_property_ui_text(
+      prop, "Use Sculpt Collision", "Enable collision with the surface while sculpting");
+  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
+  RNA_def_property_update(prop, 0, "rna_Curves_update_draw");
+
   /* attributes */
   rna_def_attributes_common(srna);