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(
+    return find_indices_from_array(span, r_indices);
-        indices_to_check, 4096, r_indices, [&](const int64_t i) { return span[i]; });
   }
 
   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. */
+  /** Solver for length and collision constraints. */
-  Array<float> segment_lengths_cu_;
+  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,
+  void comb_projected(MutableSpan<bool> r_changed_curves, const float4x4 &brush_transform)
-                      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. */
+  /** Solver for length and collision constraints. */
-  Array<float> segment_lengths_cu_;
+  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);