Curves: new Add brush

This adds a new Add brush for the new curves object type in sculpt mode.
The brush is used to insert new curves (typically hair) on the surface object.

Supported features:
* Add single curve exactly at the cursor position when `Add Amount` is 1.
* Front faces only.
* Independent interpolate shape and interpolate length settings.
* Smooth and flat shading affects curve shape interpolation.
* Spherical and projection brush.

This also adds the `surface_triangle_index` and `surface_triangle_coordinate`
attributes. Those store information about what position on the surface each
added curve is attached to:
* `surface_triangle_index` (`int`): Index of the internal triangle that a curve
  is attached to. `-1` when the curve is not attached to the surface.
* `surface_triangle_coordinate` (`float2`): First two numbers of a barycentric
  coordinate that reference a specific position within the triangle.

Ref T96444.

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

5 files changed, 806 insertions, 4 deletions

diff --git a/source/blender/editors/sculpt_paint/CMakeLists.txt b/source/blender/editors/sculpt_paint/CMakeLists.txt
index a2043c9be21..a08ec49d9c4 100644
--- a/source/blender/editors/sculpt_paint/CMakeLists.txt
+++ b/source/blender/editors/sculpt_paint/CMakeLists.txt
@@ -27,6 +27,7 @@ set(INC
 )
 
 set(SRC
+  curves_sculpt_add.cc
   curves_sculpt_comb.cc
   curves_sculpt_delete.cc
   curves_sculpt_ops.cc
diff --git a/source/blender/editors/sculpt_paint/curves_sculpt_add.cc b/source/blender/editors/sculpt_paint/curves_sculpt_add.cc
new file mode 100644
index 00000000000..e57a6e43983
--- /dev/null
+++ b/source/blender/editors/sculpt_paint/curves_sculpt_add.cc
@@ -0,0 +1,792 @@
+/* SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include <algorithm>
+
+#include "curves_sculpt_intern.hh"
+
+#include "BLI_float4x4.hh"
+#include "BLI_kdtree.h"
+#include "BLI_rand.hh"
+#include "BLI_vector.hh"
+
+#include "PIL_time.h"
+
+#include "DEG_depsgraph.h"
+
+#include "BKE_attribute_math.hh"
+#include "BKE_brush.h"
+#include "BKE_bvhutils.h"
+#include "BKE_context.h"
+#include "BKE_curves.hh"
+#include "BKE_mesh.h"
+#include "BKE_mesh_runtime.h"
+#include "BKE_paint.h"
+#include "BKE_spline.hh"
+
+#include "DNA_brush_enums.h"
+#include "DNA_brush_types.h"
+#include "DNA_curves_types.h"
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_object_types.h"
+#include "DNA_screen_types.h"
+#include "DNA_space_types.h"
+
+#include "ED_screen.h"
+#include "ED_view3d.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.
+ * re: 2D coordinates within the region.
+ */
+
+namespace blender::ed::sculpt_paint {
+
+using bke::CurvesGeometry;
+
+class AddOperation : public CurvesSculptStrokeOperation {
+ private:
+  /** Used when some data should be interpolated from existing curves. */
+  KDTree_3d *curve_roots_kdtree_ = nullptr;
+
+  friend struct AddOperationExecutor;
+
+ public:
+  ~AddOperation()
+  {
+    if (curve_roots_kdtree_ != nullptr) {
+      BLI_kdtree_3d_free(curve_roots_kdtree_);
+    }
+  }
+
+  void on_stroke_extended(bContext *C, const StrokeExtension &stroke_extension) override;
+};
+
+static void initialize_straight_curve_positions(const float3 &p1,
+                                                const float3 &p2,
+                                                MutableSpan<float3> r_positions)
+{
+  const float step = 1.0f / (float)(r_positions.size() - 1);
+  for (const int i : r_positions.index_range()) {
+    r_positions[i] = math::interpolate(p1, p2, i * step);
+  }
+}
+
+/**
+ * Utility class that actually executes the update when the stroke is updated. That's useful
+ * because it avoids passing a very large number of parameters between functions.
+ */
+struct AddOperationExecutor {
+  AddOperation *self_ = nullptr;
+  Depsgraph *depsgraph_ = nullptr;
+  Scene *scene_ = nullptr;
+  Object *object_ = nullptr;
+  ARegion *region_ = nullptr;
+  View3D *v3d_ = nullptr;
+  Curves *curves_id_ = nullptr;
+  CurvesGeometry *curves_ = nullptr;
+
+  Object *surface_ob_ = nullptr;
+  Mesh *surface_ = nullptr;
+  Span<MLoopTri> surface_looptris_;
+  Span<float3> corner_normals_su_;
+
+  CurvesSculpt *curves_sculpt_ = nullptr;
+  Brush *brush_ = nullptr;
+
+  float brush_radius_re_;
+  float2 brush_pos_re_;
+
+  bool use_front_face_;
+  bool interpolate_length_;
+  bool interpolate_shape_;
+  bool use_interpolation_;
+  float new_curve_length_;
+  int add_amount_;
+  int points_per_curve_ = 8;
+
+  /** Various matrices to convert between coordinate spaces. */
+  float4x4 curves_to_world_mat_;
+  float4x4 world_to_curves_mat_;
+  float4x4 world_to_surface_mat_;
+  float4x4 surface_to_world_mat_;
+  float4x4 surface_to_curves_mat_;
+  float4x4 surface_to_curves_normal_mat_;
+
+  BVHTreeFromMesh surface_bvh_;
+
+  int tot_old_curves_;
+  int tot_old_points_;
+
+  struct AddedPoints {
+    Vector<float3> positions_cu;
+    Vector<float3> bary_coords;
+    Vector<int> looptri_indices;
+  };
+
+  void execute(AddOperation &self, bContext *C, const StrokeExtension &stroke_extension)
+  {
+    self_ = &self;
+    depsgraph_ = CTX_data_depsgraph_pointer(C);
+    scene_ = CTX_data_scene(C);
+    object_ = CTX_data_active_object(C);
+    region_ = CTX_wm_region(C);
+    v3d_ = CTX_wm_view3d(C);
+
+    curves_id_ = static_cast<Curves *>(object_->data);
+    curves_ = &CurvesGeometry::wrap(curves_id_->geometry);
+
+    if (curves_id_->surface == nullptr || curves_id_->surface->type != OB_MESH) {
+      return;
+    }
+
+    curves_to_world_mat_ = object_->obmat;
+    world_to_curves_mat_ = curves_to_world_mat_.inverted();
+
+    surface_ob_ = curves_id_->surface;
+    surface_ = static_cast<Mesh *>(surface_ob_->data);
+    surface_to_world_mat_ = surface_ob_->obmat;
+    world_to_surface_mat_ = surface_to_world_mat_.inverted();
+    surface_to_curves_mat_ = world_to_curves_mat_ * surface_to_world_mat_;
+    surface_to_curves_normal_mat_ = surface_to_curves_mat_.inverted().transposed();
+
+    if (!CustomData_has_layer(&surface_->ldata, CD_NORMAL)) {
+      BKE_mesh_calc_normals_split(surface_);
+    }
+    corner_normals_su_ = {
+        reinterpret_cast<const float3 *>(CustomData_get_layer(&surface_->ldata, CD_NORMAL)),
+        surface_->totloop};
+
+    curves_sculpt_ = scene_->toolsettings->curves_sculpt;
+    brush_ = BKE_paint_brush(&curves_sculpt_->paint);
+    brush_radius_re_ = BKE_brush_size_get(scene_, brush_);
+    brush_pos_re_ = stroke_extension.mouse_position;
+
+    use_front_face_ = brush_->flag & BRUSH_FRONTFACE;
+    const eBrushFalloffShape falloff_shape = static_cast<eBrushFalloffShape>(
+        brush_->falloff_shape);
+    add_amount_ = std::max(0, brush_->curves_sculpt_settings->add_amount);
+    interpolate_length_ = curves_sculpt_->flag & CURVES_SCULPT_FLAG_INTERPOLATE_LENGTH;
+    interpolate_shape_ = curves_sculpt_->flag & CURVES_SCULPT_FLAG_INTERPOLATE_SHAPE;
+    use_interpolation_ = interpolate_length_ || interpolate_shape_;
+    new_curve_length_ = curves_sculpt_->curve_length;
+
+    tot_old_curves_ = curves_->curves_size();
+    tot_old_points_ = curves_->points_size();
+
+    if (add_amount_ == 0) {
+      return;
+    }
+
+    RandomNumberGenerator rng{(uint32_t)(PIL_check_seconds_timer() * 1000000.0f)};
+
+    BKE_bvhtree_from_mesh_get(&surface_bvh_, surface_, BVHTREE_FROM_LOOPTRI, 2);
+    BLI_SCOPED_DEFER([&]() { free_bvhtree_from_mesh(&surface_bvh_); });
+
+    surface_looptris_ = {BKE_mesh_runtime_looptri_ensure(surface_),
+                         BKE_mesh_runtime_looptri_len(surface_)};
+
+    /* Sample points on the surface using one of multiple strategies. */
+    AddedPoints added_points;
+    if (add_amount_ == 1) {
+      this->sample_in_center(added_points);
+    }
+    else if (falloff_shape == PAINT_FALLOFF_SHAPE_TUBE) {
+      this->sample_projected(rng, added_points);
+    }
+    else if (falloff_shape == PAINT_FALLOFF_SHAPE_SPHERE) {
+      this->sample_spherical(rng, added_points);
+    }
+    else {
+      BLI_assert_unreachable();
+    }
+
+    if (added_points.bary_coords.is_empty()) {
+      /* No new points have been added. */
+      return;
+    }
+
+    if (use_interpolation_) {
+      this->ensure_curve_roots_kdtree();
+    }
+
+    const int tot_added_curves = added_points.bary_coords.size();
+    const int tot_added_points = tot_added_curves * points_per_curve_;
+
+    curves_->resize(curves_->points_size() + tot_added_points,
+                    curves_->curves_size() + tot_added_curves);
+
+    threading::parallel_invoke([&]() { this->initialize_curve_offsets(tot_added_curves); },
+                               [&]() { this->initialize_attributes(added_points); });
+
+    DEG_id_tag_update(&curves_id_->id, ID_RECALC_GEOMETRY);
+    ED_region_tag_redraw(region_);
+  }
+
+  float3 get_bary_coords(const Mesh &mesh, const MLoopTri &looptri, const float3 position) const
+  {
+    const float3 &v0 = mesh.mvert[mesh.mloop[looptri.tri[0]].v].co;
+    const float3 &v1 = mesh.mvert[mesh.mloop[looptri.tri[1]].v].co;
+    const float3 &v2 = mesh.mvert[mesh.mloop[looptri.tri[2]].v].co;
+    float3 bary_coords;
+    interp_weights_tri_v3(bary_coords, v0, v1, v2, position);
+    return bary_coords;
+  }
+
+  /**
+   * Sample a single point exactly at the mouse position.
+   */
+  void sample_in_center(AddedPoints &r_added_points)
+  {
+    float3 ray_start_wo, ray_end_wo;
+    ED_view3d_win_to_segment_clipped(
+        depsgraph_, region_, v3d_, brush_pos_re_, ray_start_wo, ray_end_wo, true);
+    const float3 ray_start_su = world_to_surface_mat_ * ray_start_wo;
+    const float3 ray_end_su = world_to_surface_mat_ * ray_end_wo;
+    const float3 ray_direction_su = math::normalize(ray_end_su - ray_start_su);
+
+    BVHTreeRayHit ray_hit;
+    ray_hit.dist = FLT_MAX;
+    ray_hit.index = -1;
+    BLI_bvhtree_ray_cast(surface_bvh_.tree,
+                         ray_start_su,
+                         ray_direction_su,
+                         0.0f,
+                         &ray_hit,
+                         surface_bvh_.raycast_callback,
+                         &surface_bvh_);
+
+    if (ray_hit.index == -1) {
+      return;
+    }
+
+    const int looptri_index = ray_hit.index;
+    const float3 brush_pos_su = ray_hit.co;
+    const float3 bary_coords = this->get_bary_coords(
+        *surface_, surface_looptris_[looptri_index], brush_pos_su);
+
+    const float3 brush_pos_cu = surface_to_curves_mat_ * brush_pos_su;
+
+    r_added_points.positions_cu.append(brush_pos_cu);
+    r_added_points.bary_coords.append(bary_coords);
+    r_added_points.looptri_indices.append(looptri_index);
+  }
+
+  /**
+   * Sample points by shooting rays within the brush radius in the 3D view.
+   */
+  void sample_projected(RandomNumberGenerator &rng, AddedPoints &r_added_points)
+  {
+    const int max_iterations = std::max(100'000, add_amount_ * 10);
+    int current_iteration = 0;
+    while (r_added_points.bary_coords.size() < add_amount_) {
+      if (current_iteration++ >= max_iterations) {
+        break;
+      }
+
+      const float r = brush_radius_re_ * std::sqrt(rng.get_float());
+      const float angle = rng.get_float() * 2.0f * M_PI;
+      const float2 pos_re = brush_pos_re_ + r * float2(std::cos(angle), std::sin(angle));
+
+      float3 ray_start_wo, ray_end_wo;
+      ED_view3d_win_to_segment_clipped(
+          depsgraph_, region_, v3d_, pos_re, ray_start_wo, ray_end_wo, true);
+      const float3 ray_start_su = world_to_surface_mat_ * ray_start_wo;
+      const float3 ray_end_su = world_to_surface_mat_ * ray_end_wo;
+      const float3 ray_direction_su = math::normalize(ray_end_su - ray_start_su);
+
+      BVHTreeRayHit ray_hit;
+      ray_hit.dist = FLT_MAX;
+      ray_hit.index = -1;
+      BLI_bvhtree_ray_cast(surface_bvh_.tree,
+                           ray_start_su,
+                           ray_direction_su,
+                           0.0f,
+                           &ray_hit,
+                           surface_bvh_.raycast_callback,
+                           &surface_bvh_);
+
+      if (ray_hit.index == -1) {
+        continue;
+      }
+
+      if (use_front_face_) {
+        const float3 normal_su = ray_hit.no;
+        if (math::dot(ray_direction_su, normal_su) >= 0.0f) {
+          continue;
+        }
+      }
+
+      const int looptri_index = ray_hit.index;
+      const float3 pos_su = ray_hit.co;
+
+      const float3 bary_coords = this->get_bary_coords(
+          *surface_, surface_looptris_[looptri_index], pos_su);
+
+      const float3 pos_cu = surface_to_curves_mat_ * pos_su;
+
+      r_added_points.positions_cu.append(pos_cu);
+      r_added_points.bary_coords.append(bary_coords);
+      r_added_points.looptri_indices.append(looptri_index);
+    }
+  }
+
+  /**
+   * Sample points in a 3D sphere around the surface position that the mouse hovers over.
+   */
+  void sample_spherical(RandomNumberGenerator &rng, AddedPoints &r_added_points)
+  {
+    /* Find ray that starts in the center of the brush. */
+    float3 brush_ray_start_wo, brush_ray_end_wo;
+    ED_view3d_win_to_segment_clipped(
+        depsgraph_, region_, v3d_, brush_pos_re_, brush_ray_start_wo, brush_ray_end_wo, true);
+    const float3 brush_ray_start_su = world_to_surface_mat_ * brush_ray_start_wo;
+    const float3 brush_ray_end_su = world_to_surface_mat_ * brush_ray_end_wo;
+    const float3 brush_ray_direction_su = math::normalize(brush_ray_end_su - brush_ray_start_su);
+
+    /* Find ray that starts on the boundary of the brush. That is used to compute the brush radius
+     * in 3D. */
+    float3 brush_radius_ray_start_wo, brush_radius_ray_end_wo;
+    ED_view3d_win_to_segment_clipped(depsgraph_,
+                                     region_,
+                                     v3d_,
+                                     brush_pos_re_ + float2(brush_radius_re_, 0),
+                                     brush_radius_ray_start_wo,
+                                     brush_radius_ray_end_wo,
+                                     true);
+    const float3 brush_radius_ray_start_su = world_to_surface_mat_ * brush_radius_ray_start_wo;
+    const float3 brush_radius_ray_end_su = world_to_surface_mat_ * brush_radius_ray_end_wo;
+
+    BVHTreeRayHit ray_hit;
+    ray_hit.dist = FLT_MAX;
+    ray_hit.index = -1;
+    BLI_bvhtree_ray_cast(surface_bvh_.tree,
+                         brush_ray_start_su,
+                         brush_ray_direction_su,
+                         0.0f,
+                         &ray_hit,
+                         surface_bvh_.raycast_callback,
+                         &surface_bvh_);
+
+    if (ray_hit.index == -1) {
+      return;
+    }
+
+    /* Compute brush radius. */
+    const float3 brush_pos_su = ray_hit.co;
+    const float brush_radius_su = dist_to_line_v3(
+        brush_pos_su, brush_radius_ray_start_su, brush_radius_ray_end_su);
+    const float brush_radius_sq_su = pow2f(brush_radius_su);
+
+    /* Find surface triangles within brush radius. */
+    Vector<int> looptri_indices;
+    if (use_front_face_) {
+      BLI_bvhtree_range_query_cpp(
+          *surface_bvh_.tree,
+          brush_pos_su,
+          brush_radius_su,
+          [&](const int index, const float3 &UNUSED(co), const float UNUSED(dist_sq)) {
+            const MLoopTri &looptri = surface_looptris_[index];
+            const float3 v0_su = surface_->mvert[surface_->mloop[looptri.tri[0]].v].co;
+            const float3 v1_su = surface_->mvert[surface_->mloop[looptri.tri[1]].v].co;
+            const float3 v2_su = surface_->mvert[surface_->mloop[looptri.tri[2]].v].co;
+            float3 normal_su;
+            normal_tri_v3(normal_su, v0_su, v1_su, v2_su);
+            if (math::dot(normal_su, brush_ray_direction_su) >= 0.0f) {
+              return;
+            }
+            looptri_indices.append(index);
+          });
+    }
+    else {
+      BLI_bvhtree_range_query_cpp(
+          *surface_bvh_.tree,
+          brush_pos_su,
+          brush_radius_su,
+          [&](const int index, const float3 &UNUSED(co), const float UNUSED(dist_sq)) {
+            looptri_indices.append(index);
+          });
+    }
+
+    /* Density used for sampling points. This does not have to be exact, because the loop below
+     * automatically runs until enough samples have been found. If too many samples are found, some
+     * will be discarded afterwards. */
+    const float brush_plane_area_su = M_PI * brush_radius_sq_su;
+    const float approximate_density_su = add_amount_ / brush_plane_area_su;
+
+    /* Used for switching between two triangle sampling strategies. */
+    const float area_threshold = brush_plane_area_su;
+
+    /* Usually one or two iterations should be enough. */
+    const int max_iterations = 5;
+    int current_iteration = 0;
+
+    while (r_added_points.bary_coords.size() < add_amount_) {
+      if (current_iteration++ >= max_iterations) {
+        break;
+      }
+
+      for (const int looptri_index : looptri_indices) {
+        const MLoopTri &looptri = surface_looptris_[looptri_index];
+
+        const float3 v0_su = surface_->mvert[surface_->mloop[looptri.tri[0]].v].co;
+        const float3 v1_su = surface_->mvert[surface_->mloop[looptri.tri[1]].v].co;
+        const float3 v2_su = surface_->mvert[surface_->mloop[looptri.tri[2]].v].co;
+
+        const float looptri_area_su = area_tri_v3(v0_su, v1_su, v2_su);
+
+        if (looptri_area_su < area_threshold) {
+          /* The triangle is small compared to the brush radius. Sample by generating random
+           * barycentric coordinates. */
+          const int amount = rng.round_probabilistic(approximate_density_su * looptri_area_su);
+          for ([[maybe_unused]] const int i : IndexRange(amount)) {
+            const float3 bary_coord = rng.get_barycentric_coordinates();
+            const float3 point_pos_su = attribute_math::mix3(bary_coord, v0_su, v1_su, v2_su);
+            const float distance_to_brush_sq_su = math::distance_squared(point_pos_su,
+                                                                         brush_pos_su);
+            if (distance_to_brush_sq_su > brush_radius_sq_su) {
+              continue;
+            }
+
+            r_added_points.bary_coords.append(bary_coord);
+            r_added_points.looptri_indices.append(looptri_index);
+            r_added_points.positions_cu.append(surface_to_curves_mat_ * point_pos_su);
+          }
+        }
+        else {
+          /* The triangle is large compared to the brush radius. Sample by generating random points
+           * on the triangle plane within the brush radius. */
+          float3 normal_su;
+          normal_tri_v3(normal_su, v0_su, v1_su, v2_su);
+
+          float3 brush_pos_proj_su = brush_pos_su;
+          project_v3_plane(brush_pos_proj_su, normal_su, v0_su);
+
+          const float proj_distance_sq_su = math::distance_squared(brush_pos_proj_su,
+                                                                   brush_pos_su);
+          const float brush_radius_factor_sq = 1.0f -
+                                               std::min(1.0f,
+                                                        proj_distance_sq_su / brush_radius_sq_su);
+          const float radius_proj_sq_su = brush_radius_sq_su * brush_radius_factor_sq;
+          const float radius_proj_su = std::sqrt(radius_proj_sq_su);
+          const float circle_area_su = M_PI * radius_proj_su;
+
+          const int amount = rng.round_probabilistic(approximate_density_su * circle_area_su);
+
+          const float3 axis_1_su = math::normalize(v1_su - v0_su) * radius_proj_su;
+          const float3 axis_2_su = math::normalize(math::cross(
+                                       axis_1_su, math::cross(axis_1_su, v2_su - v0_su))) *
+                                   radius_proj_su;
+
+          for ([[maybe_unused]] const int i : IndexRange(amount)) {
+            const float r = std::sqrt(rng.get_float());
+            const float angle = rng.get_float() * 2.0f * M_PI;
+            const float x = r * std::cos(angle);
+            const float y = r * std::sin(angle);
+            const float3 point_pos_su = brush_pos_proj_su + axis_1_su * x + axis_2_su * y;
+            if (!isect_point_tri_prism_v3(point_pos_su, v0_su, v1_su, v2_su)) {
+              /* Sampled point is not in the triangle. */
+              continue;
+            }
+
+            float3 bary_coord;
+            interp_weights_tri_v3(bary_coord, v0_su, v1_su, v2_su, point_pos_su);
+
+            r_added_points.bary_coords.append(bary_coord);
+            r_added_points.looptri_indices.append(looptri_index);
+            r_added_points.positions_cu.append(surface_to_curves_mat_ * point_pos_su);
+          }
+        }
+      }
+    }
+
+    /* Remove samples when there are too many. */
+    while (r_added_points.bary_coords.size() > add_amount_) {
+      const int index_to_remove = rng.get_int32(r_added_points.bary_coords.size());
+      r_added_points.bary_coords.remove_and_reorder(index_to_remove);
+      r_added_points.looptri_indices.remove_and_reorder(index_to_remove);
+      r_added_points.positions_cu.remove_and_reorder(index_to_remove);
+    }
+  }
+
+  void ensure_curve_roots_kdtree()
+  {
+    if (self_->curve_roots_kdtree_ == nullptr) {
+      self_->curve_roots_kdtree_ = BLI_kdtree_3d_new(curves_->curves_size());
+      for (const int curve_i : curves_->curves_range()) {
+        const int root_point_i = curves_->offsets()[curve_i];
+        const float3 &root_pos_cu = curves_->positions()[root_point_i];
+        BLI_kdtree_3d_insert(self_->curve_roots_kdtree_, curve_i, root_pos_cu);
+      }
+      BLI_kdtree_3d_balance(self_->curve_roots_kdtree_);
+    }
+  }
+
+  void initialize_curve_offsets(const int tot_added_curves)
+  {
+    MutableSpan<int> offsets = curves_->offsets();
+    threading::parallel_for(IndexRange(tot_added_curves), 1024, [&](const IndexRange range) {
+      for (const int i : range) {
+        const int curve_i = tot_old_curves_ + i;
+        offsets[curve_i + 1] = tot_old_points_ + (i + 1) * points_per_curve_;
+      }
+    });
+  }
+
+  struct NeighborInfo {
+    /* Curve index of the neighbor. */
+    int index;
+    /* The weights of all neighbors of a new curve add up to 1. */
+    float weight;
+  };
+  static constexpr int max_neighbors = 5;
+  using NeighborsVector = Vector<NeighborInfo, max_neighbors>;
+
+  void initialize_attributes(const AddedPoints &added_points)
+  {
+    Array<NeighborsVector> neighbors_per_curve;
+    if (use_interpolation_) {
+      neighbors_per_curve = this->find_curve_neighbors(added_points);
+    }
+
+    Array<float> new_lengths_cu(added_points.bary_coords.size());
+    if (interpolate_length_) {
+      this->interpolate_lengths(neighbors_per_curve, new_lengths_cu);
+    }
+    else {
+      new_lengths_cu.fill(new_curve_length_);
+    }
+
+    Array<float3> new_normals_su = this->compute_normals_for_added_curves_su(added_points);
+    this->initialize_surface_attachment(added_points);
+
+    if (interpolate_shape_) {
+      this->initialize_position_with_interpolation(
+          added_points, neighbors_per_curve, new_normals_su, new_lengths_cu);
+    }
+    else {
+      this->initialize_position_without_interpolation(
+          added_points, new_lengths_cu, new_normals_su);
+    }
+  }
+
+  Array<NeighborsVector> find_curve_neighbors(const AddedPoints &added_points)
+  {
+    const int tot_added_curves = added_points.bary_coords.size();
+    Array<NeighborsVector> neighbors_per_curve(tot_added_curves);
+    threading::parallel_for(IndexRange(tot_added_curves), 128, [&](const IndexRange range) {
+      for (const int i : range) {
+        const float3 root_cu = added_points.positions_cu[i];
+        std::array<KDTreeNearest_3d, max_neighbors> nearest_n;
+        const int found_neighbors = BLI_kdtree_3d_find_nearest_n(
+            self_->curve_roots_kdtree_, root_cu, nearest_n.data(), max_neighbors);
+        float tot_weight = 0.0f;
+        for (const int neighbor_i : IndexRange(found_neighbors)) {
+          KDTreeNearest_3d &nearest = nearest_n[neighbor_i];
+          const float weight = 1.0f / std::max(nearest.dist, 0.00001f);
+          tot_weight += weight;
+          neighbors_per_curve[i].append({nearest.index, weight});
+        }
+        /* Normalize weights. */
+        for (NeighborInfo &neighbor : neighbors_per_curve[i]) {
+          neighbor.weight /= tot_weight;
+        }
+      }
+    });
+    return neighbors_per_curve;
+  }
+
+  void interpolate_lengths(const Span<NeighborsVector> neighbors_per_curve,
+                           MutableSpan<float> r_lengths)
+  {
+    const Span<float3> positions_cu = curves_->positions();
+
+    threading::parallel_for(r_lengths.index_range(), 128, [&](const IndexRange range) {
+      for (const int added_curve_i : range) {
+        const Span<NeighborInfo> neighbors = neighbors_per_curve[added_curve_i];
+        float length_sum = 0.0f;
+        for (const NeighborInfo &neighbor : neighbors) {
+          const IndexRange neighbor_points = curves_->range_for_curve(neighbor.index);
+          float neighbor_length = 0.0f;
+          const int tot_segments = neighbor_points.size() - 1;
+          for (const int segment_i : IndexRange(tot_segments)) {
+            const float3 &p1 = positions_cu[neighbor_points[segment_i]];
+            const float3 &p2 = positions_cu[neighbor_points[segment_i] + 1];
+            neighbor_length += math::distance(p1, p2);
+          }
+          length_sum += neighbor.weight * neighbor_length;
+        }
+        const float length = neighbors.is_empty() ? new_curve_length_ : length_sum;
+        r_lengths[added_curve_i] = length;
+      }
+    });
+  }
+
+  float3 compute_point_normal_su(const int looptri_index, const float3 &bary_coord)
+  {
+    const MLoopTri &looptri = surface_looptris_[looptri_index];
+    const int l0 = looptri.tri[0];
+    const int l1 = looptri.tri[1];
+    const int l2 = looptri.tri[2];
+
+    const float3 &l0_normal_su = corner_normals_su_[l0];
+    const float3 &l1_normal_su = corner_normals_su_[l1];
+    const float3 &l2_normal_su = corner_normals_su_[l2];
+
+    const float3 normal_su = math::normalize(
+        attribute_math::mix3(bary_coord, l0_normal_su, l1_normal_su, l2_normal_su));
+    return normal_su;
+  }
+
+  Array<float3> compute_normals_for_added_curves_su(const AddedPoints &added_points)
+  {
+    Array<float3> normals_su(added_points.bary_coords.size());
+    threading::parallel_for(normals_su.index_range(), 256, [&](const IndexRange range) {
+      for (const int i : range) {
+        const int looptri_index = added_points.looptri_indices[i];
+        const float3 &bary_coord = added_points.bary_coords[i];
+        normals_su[i] = this->compute_point_normal_su(looptri_index, bary_coord);
+      }
+    });
+    return normals_su;
+  }
+
+  void initialize_surface_attachment(const AddedPoints &added_points)
+  {
+    MutableSpan<int> surface_triangle_indices = curves_->surface_triangle_indices();
+    MutableSpan<float2> surface_triangle_coords = curves_->surface_triangle_coords();
+    threading::parallel_for(
+        added_points.bary_coords.index_range(), 1024, [&](const IndexRange range) {
+          for (const int i : range) {
+            const int curve_i = tot_old_curves_ + i;
+            surface_triangle_indices[curve_i] = added_points.looptri_indices[i];
+            surface_triangle_coords[curve_i] = float2(added_points.bary_coords[i]);
+          }
+        });
+  }
+
+  /**
+   * Initialize new curves so that they are just a straight line in the normal direction.
+   */
+  void initialize_position_without_interpolation(const AddedPoints &added_points,
+                                                 const Span<float> lengths_cu,
+                                                 const MutableSpan<float3> normals_su)
+  {
+    MutableSpan<float3> positions_cu = curves_->positions();
+
+    threading::parallel_for(
+        added_points.bary_coords.index_range(), 256, [&](const IndexRange range) {
+          for (const int i : range) {
+            const int first_point_i = tot_old_points_ + i * points_per_curve_;
+            const float3 &root_cu = added_points.positions_cu[i];
+            const float length = lengths_cu[i];
+            const float3 &normal_su = normals_su[i];
+            const float3 normal_cu = math::normalize(surface_to_curves_normal_mat_ * normal_su);
+            const float3 tip_cu = root_cu + length * normal_cu;
+
+            initialize_straight_curve_positions(
+                root_cu, tip_cu, positions_cu.slice(first_point_i, points_per_curve_));
+          }
+        });
+  }
+
+  /**
+   * Use neighboring curves to determine the shape.
+   */
+  void initialize_position_with_interpolation(const AddedPoints &added_points,
+                                              const Span<NeighborsVector> neighbors_per_curve,
+                                              const Span<float3> new_normals_su,
+                                              const Span<float> new_lengths_cu)
+  {
+    MutableSpan<float3> positions_cu = curves_->positions();
+    const Span<int> surface_triangle_indices = curves_->surface_triangle_indices();
+    const Span<float2> surface_triangle_coords = curves_->surface_triangle_coords();
+
+    threading::parallel_for(
+        added_points.bary_coords.index_range(), 256, [&](const IndexRange range) {
+          for (const int i : range) {
+            const Span<NeighborInfo> neighbors = neighbors_per_curve[i];
+
+            const float length_cu = new_lengths_cu[i];
+            const float3 &normal_su = new_normals_su[i];
+            const float3 normal_cu = math::normalize(surface_to_curves_normal_mat_ * normal_su);
+
+            const float3 &root_cu = added_points.positions_cu[i];
+            const int first_point_i = tot_old_points_ + i * points_per_curve_;
+
+            if (neighbors.is_empty()) {
+              /* If there are no neighbors, just make a straight line. */
+              const float3 tip_cu = root_cu + length_cu * normal_cu;
+              initialize_straight_curve_positions(
+                  root_cu, tip_cu, positions_cu.slice(first_point_i, points_per_curve_));
+              continue;
+            }
+
+            positions_cu.slice(first_point_i, points_per_curve_).fill(root_cu);
+
+            for (const NeighborInfo &neighbor : neighbors) {
+              const int neighbor_curve_i = neighbor.index;
+              const int neighbor_looptri_index = surface_triangle_indices[neighbor_curve_i];
+
+              float3 neighbor_bary_coord{surface_triangle_coords[neighbor_curve_i]};
+              neighbor_bary_coord.z = 1.0f - neighbor_bary_coord.x - neighbor_bary_coord.y;
+
+              const float3 neighbor_normal_su = this->compute_point_normal_su(
+                  neighbor_looptri_index, neighbor_bary_coord);
+              const float3 neighbor_normal_cu = math::normalize(surface_to_curves_normal_mat_ *
+                                                                neighbor_normal_su);
+
+              /* The rotation matrix used to transform relative coordinates of the neighbor curve
+               * to the new curve. */
+              float normal_rotation_cu[3][3];
+              rotation_between_vecs_to_mat3(normal_rotation_cu, neighbor_normal_cu, normal_cu);
+
+              const IndexRange neighbor_points = curves_->range_for_curve(neighbor_curve_i);
+              const float3 &neighbor_root_cu = positions_cu[neighbor_points[0]];
+
+              /* Use a temporary #PolySpline, because that's the easiest way to resample an
+               * existing curve right now. Resampling is necessary if the length of the new curve
+               * does not match the length of the neighbors or the number of handle points is
+               * different. */
+              PolySpline neighbor_spline;
+              neighbor_spline.resize(neighbor_points.size());
+              neighbor_spline.positions().copy_from(positions_cu.slice(neighbor_points));
+              neighbor_spline.mark_cache_invalid();
+
+              const float neighbor_length_cu = neighbor_spline.length();
+              const float length_factor = std::min(1.0f, length_cu / neighbor_length_cu);
+
+              const float resample_factor = (1.0f / (points_per_curve_ - 1.0f)) * length_factor;
+              for (const int j : IndexRange(points_per_curve_)) {
+                const Spline::LookupResult lookup = neighbor_spline.lookup_evaluated_factor(
+                    j * resample_factor);
+                const float index_factor = lookup.evaluated_index + lookup.factor;
+                float3 p;
+                neighbor_spline.sample_with_index_factors<float3>(
+                    neighbor_spline.positions(), {&index_factor, 1}, {&p, 1});
+                const float3 relative_coord = p - neighbor_root_cu;
+                float3 rotated_relative_coord = relative_coord;
+                mul_m3_v3(normal_rotation_cu, rotated_relative_coord);
+                positions_cu[first_point_i + j] += neighbor.weight * rotated_relative_coord;
+              }
+            }
+          }
+        });
+  }
+};
+
+void AddOperation::on_stroke_extended(bContext *C, const StrokeExtension &stroke_extension)
+{
+  AddOperationExecutor executor;
+  executor.execute(*this, C, stroke_extension);
+}
+
+std::unique_ptr<CurvesSculptStrokeOperation> new_add_operation()
+{
+  return std::make_unique<AddOperation>();
+}
+
+}  // namespace blender::ed::sculpt_paint
diff --git a/source/blender/editors/sculpt_paint/curves_sculpt_intern.hh b/source/blender/editors/sculpt_paint/curves_sculpt_intern.hh
index 14d7ec01251..d1de69f36b6 100644
--- a/source/blender/editors/sculpt_paint/curves_sculpt_intern.hh
+++ b/source/blender/editors/sculpt_paint/curves_sculpt_intern.hh
@@ -22,6 +22,7 @@ class CurvesSculptStrokeOperation {
   virtual void on_stroke_extended(bContext *C, const StrokeExtension &stroke_extension) = 0;
 };
 
+std::unique_ptr<CurvesSculptStrokeOperation> new_add_operation();
 std::unique_ptr<CurvesSculptStrokeOperation> new_comb_operation();
 std::unique_ptr<CurvesSculptStrokeOperation> new_delete_operation();
 std::unique_ptr<CurvesSculptStrokeOperation> new_snake_hook_operation();
diff --git a/source/blender/editors/sculpt_paint/curves_sculpt_ops.cc b/source/blender/editors/sculpt_paint/curves_sculpt_ops.cc
index 06fdeabf3e1..b9a019a012c 100644
--- a/source/blender/editors/sculpt_paint/curves_sculpt_ops.cc
+++ b/source/blender/editors/sculpt_paint/curves_sculpt_ops.cc
@@ -186,7 +186,7 @@ class ShrinkOperation : public CurvesSculptStrokeOperation {
   }
 };
 
-class AddOperation : public CurvesSculptStrokeOperation {
+class DensityAddOperation : public CurvesSculptStrokeOperation {
  private:
   /** Contains the root points of the curves that existed before this operation started. */
   KDTree_3d *old_kdtree_ = nullptr;
@@ -207,7 +207,7 @@ class AddOperation : public CurvesSculptStrokeOperation {
   };
 
  public:
-  ~AddOperation() override
+  ~DensityAddOperation() override
   {
     if (old_kdtree_ != nullptr) {
       BLI_kdtree_3d_free(old_kdtree_);
@@ -624,8 +624,10 @@ static std::unique_ptr<CurvesSculptStrokeOperation> start_brush_operation(bConte
       return new_delete_operation();
     case CURVES_SCULPT_TOOL_SNAKE_HOOK:
       return new_snake_hook_operation();
+    case CURVES_SCULPT_TOOL_ADD:
+      return new_add_operation();
     case CURVES_SCULPT_TOOL_TEST1:
-      return std::make_unique<AddOperation>();
+      return std::make_unique<DensityAddOperation>();
     case CURVES_SCULPT_TOOL_TEST2:
       return std::make_unique<ShrinkOperation>();
   }
diff --git a/source/blender/editors/sculpt_paint/paint_stroke.c b/source/blender/editors/sculpt_paint/paint_stroke.c
index c0bf89107e0..0f7b8ad1f3d 100644
--- a/source/blender/editors/sculpt_paint/paint_stroke.c
+++ b/source/blender/editors/sculpt_paint/paint_stroke.c
@@ -986,6 +986,11 @@ static void stroke_done(bContext *C, wmOperator *op, PaintStroke *stroke)
   paint_stroke_free(C, op, stroke);
 }
 
+static bool curves_sculpt_brush_uses_spacing(const eBrushCurvesSculptTool tool)
+{
+  return ELEM(tool, CURVES_SCULPT_TOOL_ADD);
+}
+
 bool paint_space_stroke_enabled(Brush *br, ePaintMode mode)
 {
   if ((br->flag & BRUSH_SPACE) == 0) {
@@ -1000,7 +1005,8 @@ bool paint_space_stroke_enabled(Brush *br, ePaintMode mode)
     return true;
   }
 
-  if (mode == PAINT_MODE_SCULPT_CURVES) {
+  if (mode == PAINT_MODE_SCULPT_CURVES &&
+      !curves_sculpt_brush_uses_spacing(br->curves_sculpt_tool)) {
     return false;
   }