Curves: add brush to add curves on surface

This adds a prototype for the first brush that can add new curves by
painting on a surface. Note that this can only be used when the curves
object has a surface object set in the properties panel.

The brush can take minimum distance into account. This allows
distributing curves with a somewhat consistent density.

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

release/scripts/startup/bl_ui/space_view3d.py

@@ -473,7 +473,8 @@ class _draw_tool_settings_context_mode:
         if (tool is None) or (not tool.has_datablock):
             return False
 
-        paint = context.tool_settings.curves_sculpt
+        tool_settings = context.tool_settings
+        paint = tool_settings.curves_sculpt
         layout.template_ID_preview(paint, "brush", rows=3, cols=8, hide_buttons=True)
 
         brush = paint.brush
@@ -500,6 +501,10 @@ class _draw_tool_settings_context_mode:
             header=True
         )
 
+        if brush.curves_sculpt_tool == "TEST3":
+            layout.prop(tool_settings.curves_sculpt, "distance")
+
+
 
 class VIEW3D_HT_header(Header):
     bl_space_type = 'VIEW_3D'

source/blender/blenlib/BLI_hash.hh

@@ -148,6 +148,13 @@ template<> struct DefaultHash<float> {
   }
 };
 
+template<> struct DefaultHash<double> {
+  uint64_t operator()(double value) const
+  {
+    return *reinterpret_cast<uint64_t *>(&value);
+  }
+};
+
 template<> struct DefaultHash<bool> {
   uint64_t operator()(bool value) const
   {

source/blender/editors/sculpt_paint/curves_sculpt_ops.cc

@@ -2,9 +2,14 @@
 
 #include "BLI_utildefines.h"
 
+#include "BKE_attribute_math.hh"
 #include "BKE_brush.h"
+#include "BKE_bvhutils.h"
 #include "BKE_context.h"
 #include "BKE_curves.hh"
+#include "BKE_lib_id.h"
+#include "BKE_mesh.h"
+#include "BKE_mesh_runtime.h"
 #include "BKE_paint.h"
 
 #include "WM_api.h"
@@ -19,12 +24,18 @@
 
 #include "DNA_brush_types.h"
 #include "DNA_curves_types.h"
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
 #include "DNA_screen_types.h"
 
 #include "RNA_access.h"
 
 #include "BLI_index_mask_ops.hh"
+#include "BLI_kdtree.h"
 #include "BLI_math_vector.hh"
+#include "BLI_rand.hh"
+
+#include "PIL_time.h"
 
 #include "curves_sculpt_intern.h"
 #include "paint_intern.h"
@@ -79,7 +90,7 @@ class DeleteOperation : public CurvesSculptStrokeOperation {
   float2 last_mouse_position_;
 
  public:
-  void on_stroke_extended(bContext *C, const StrokeExtension &stroke_extension)
+  void on_stroke_extended(bContext *C, const StrokeExtension &stroke_extension) override
   {
     Scene &scene = *CTX_data_scene(C);
     Object &object = *CTX_data_active_object(C);
@@ -146,7 +157,7 @@ class MoveOperation : public CurvesSculptStrokeOperation {
   float2 last_mouse_position_;
 
  public:
-  void on_stroke_extended(bContext *C, const StrokeExtension &stroke_extension)
+  void on_stroke_extended(bContext *C, const StrokeExtension &stroke_extension) override
   {
     Scene &scene = *CTX_data_scene(C);
     Object &object = *CTX_data_active_object(C);
@@ -201,6 +212,438 @@ class MoveOperation : public CurvesSculptStrokeOperation {
   }
 };
 
+class AddOperation : public CurvesSculptStrokeOperation {
+ private:
+  /** Contains the root points of the curves that existed before this operation started. */
+  KDTree_3d *old_kdtree_ = nullptr;
+  /** Number of points in the kdtree above. */
+  int old_kdtree_size_ = 0;
+
+  /**
+   * Indicates that the corresponding curve has already been created and can't be changed by this
+   * operation anymore.
+   */
+  static constexpr int ExistsAlreadyIndex = INT32_MAX;
+
+  struct NewPointsData {
+    Vector<float3> bary_coords;
+    Vector<int> looptri_indices;
+    Vector<float3> positions;
+    Vector<float3> normals;
+  };
+
+ public:
+  ~AddOperation()
+  {
+    if (old_kdtree_ != nullptr) {
+      BLI_kdtree_3d_free(old_kdtree_);
+    }
+  }
+
+  void on_stroke_extended(bContext *C, const StrokeExtension &stroke_extension) override
+  {
+    Depsgraph &depsgraph = *CTX_data_depsgraph_pointer(C);
+    Scene &scene = *CTX_data_scene(C);
+    Object &object = *CTX_data_active_object(C);
+    ARegion *region = CTX_wm_region(C);
+    View3D *v3d = CTX_wm_view3d(C);
+
+    Curves &curves_id = *static_cast<Curves *>(object.data);
+    CurvesGeometry &curves = CurvesGeometry::wrap(curves_id.geometry);
+
+    if (curves_id.surface == nullptr || curves_id.surface->type != OB_MESH) {
+      return;
+    }
+
+    const Object &surface_ob = *curves_id.surface;
+    const Mesh &surface = *static_cast<const Mesh *>(surface_ob.data);
+    const float4x4 surface_ob_mat = surface_ob.obmat;
+    const float4x4 surface_ob_imat = surface_ob_mat.inverted();
+
+    ToolSettings &tool_settings = *scene.toolsettings;
+    CurvesSculpt &curves_sculpt = *tool_settings.curves_sculpt;
+    Brush &brush = *BKE_paint_brush(&curves_sculpt.paint);
+    const float brush_radius_screen = BKE_brush_size_get(&scene, &brush);
+    const float strength = BKE_brush_alpha_get(&scene, &brush);
+    const float minimum_distance = curves_sculpt.distance;
+
+    /* This is the main ray that is used to determine the brush position in 3D space. */
+    float3 ray_start, ray_end;
+    ED_view3d_win_to_segment_clipped(
+        &depsgraph, region, v3d, stroke_extension.mouse_position, ray_start, ray_end, true);
+    ray_start = surface_ob_imat * ray_start;
+    ray_end = surface_ob_imat * ray_end;
+    const float3 ray_direction = math::normalize(ray_end - ray_start);
+
+    /* This ray is used to determine the brush radius in 3d space. */
+    float3 offset_ray_start, offset_ray_end;
+    ED_view3d_win_to_segment_clipped(&depsgraph,
+                                     region,
+                                     v3d,
+                                     stroke_extension.mouse_position +
+                                         float2(0, brush_radius_screen),
+                                     offset_ray_start,
+                                     offset_ray_end,
+                                     true);
+    offset_ray_start = surface_ob_imat * offset_ray_start;
+    offset_ray_end = surface_ob_imat * offset_ray_end;
+
+    float4x4 ob_imat;
+    invert_m4_m4(ob_imat.values, object.obmat);
+
+    const float4x4 transform = ob_imat * surface_ob_mat;
+
+    BVHTreeFromMesh bvhtree;
+    BKE_bvhtree_from_mesh_get(&bvhtree, &surface, BVHTREE_FROM_LOOPTRI, 2);
+
+    /* Do a raycast against the surface object to find the brush position. */
+    BVHTreeRayHit ray_hit;
+    ray_hit.dist = FLT_MAX;
+    ray_hit.index = -1;
+    BLI_bvhtree_ray_cast(bvhtree.tree,
+                         ray_start,
+                         ray_direction,
+                         0.0f,
+                         &ray_hit,
+                         bvhtree.raycast_callback,
+                         &bvhtree);
+
+    if (ray_hit.index == -1) {
+      /* The ray did not hit the surface. */
+      free_bvhtree_from_mesh(&bvhtree);
+      return;
+    }
+    /* Brush position in the space of the surface object. */
+    const float3 brush_pos_3d_surface = ray_hit.co;
+    const float brush_radius_3d_surface = dist_to_line_v3(
+        brush_pos_3d_surface, offset_ray_start, offset_ray_end);
+
+    /* Brush position in the space of the curves object. */
+    const float3 brush_pos_3d_curves = transform * brush_pos_3d_surface;
+    const float brush_radius_3d_curves = dist_to_line_v3(
+        brush_pos_3d_curves, transform * offset_ray_start, transform * offset_ray_end);
+
+    Vector<int> looptri_indices = this->find_looptri_indices_to_consider(
+        bvhtree, brush_pos_3d_surface, brush_radius_3d_surface);
+
+    free_bvhtree_from_mesh(&bvhtree);
+
+    if (old_kdtree_ == nullptr && minimum_distance > 0.0f) {
+      old_kdtree_ = this->kdtree_from_curve_roots_and_positions(curves, curves.curves_range(), {});
+      old_kdtree_size_ = curves.curves_size();
+    }
+
+    float density;
+    if (minimum_distance > 0.0f) {
+      /* Estimate the sampling density based on the target minimum distance. */
+      density = strength * pow2f(1.0f / minimum_distance);
+    }
+    else {
+      /* Sample a somewhat constant amount of points based on the strength. */
+      const float brush_circle_area_3d = M_PI * pow2f(brush_radius_3d_curves);
+      density = strength * 100.0f / brush_circle_area_3d;
+    }
+
+    NewPointsData new_points = this->sample_new_points(density,
+                                                       minimum_distance,
+                                                       brush_radius_3d_curves,
+                                                       brush_pos_3d_curves,
+                                                       looptri_indices,
+                                                       transform,
+                                                       surface);
+    if (minimum_distance > 0.0f) {
+      this->eliminate_too_close_points(new_points, curves, minimum_distance);
+    }
+    this->insert_new_curves(new_points, curves);
+
+    DEG_id_tag_update(&curves_id.id, ID_RECALC_GEOMETRY);
+    ED_region_tag_redraw(region);
+  }
+
+ private:
+  Vector<int> find_looptri_indices_to_consider(BVHTreeFromMesh &bvhtree,
+                                               const float3 &brush_pos,
+                                               const float brush_radius_3d)
+  {
+    Vector<int> looptri_indices;
+
+    struct RangeQueryUserData {
+      Vector<int> &indices;
+    } range_query_user_data = {looptri_indices};
+
+    BLI_bvhtree_range_query(
+        bvhtree.tree,
+        brush_pos,
+        brush_radius_3d,
+        [](void *userdata, int index, const float co[3], float dist_sq) {
+          UNUSED_VARS(co, dist_sq);
+          RangeQueryUserData &data = *static_cast<RangeQueryUserData *>(userdata);
+          data.indices.append(index);
+        },
+        &range_query_user_data);
+
+    return looptri_indices;
+  }
+
+  KDTree_3d *kdtree_from_curve_roots_and_positions(const CurvesGeometry &curves,
+                                                   const IndexRange curves_range,
+                                                   Span<float3> extra_positions)
+  {
+    const int tot_points = curves_range.size() + extra_positions.size();
+    KDTree_3d *kdtree = BLI_kdtree_3d_new(tot_points);
+    for (const int curve_i : curves_range) {
+      const int first_point_i = curves.offsets()[curve_i];
+      const float3 root_position = curves.positions()[first_point_i];
+      BLI_kdtree_3d_insert(kdtree, ExistsAlreadyIndex, root_position);
+    }
+    for (const int i : extra_positions.index_range()) {
+      BLI_kdtree_3d_insert(kdtree, i, extra_positions[i]);
+    }
+    BLI_kdtree_3d_balance(kdtree);
+    return kdtree;
+  }
+
+  int float_to_int_amount(float amount_f, RandomNumberGenerator &rng)
+  {
+    const float add_probability = fractf(amount_f);
+    const bool add_point = add_probability > rng.get_float();
+    return (int)amount_f + (int)add_point;
+  }
+
+  bool is_too_close_to_existing_point(const float3 position, const float minimum_distance) const
+  {
+    if (old_kdtree_ == nullptr) {
+      return false;
+    }
+    KDTreeNearest_3d nearest;
+    nearest.index = -1;
+    BLI_kdtree_3d_find_nearest(old_kdtree_, position, &nearest);
+    if (nearest.index >= 0 && nearest.dist < minimum_distance) {
+      return true;
+    }
+    return false;
+  }
+
+  NewPointsData sample_new_points(const float density,
+                                  const float minimum_distance,
+                                  const float brush_radius_3d,
+                                  const float3 &brush_pos,
+                                  const Span<int> looptri_indices,
+                                  const float4x4 &transform,
+                                  const Mesh &surface)
+  {
+    const float brush_radius_3d_sq = brush_radius_3d * brush_radius_3d;
+    const float area_threshold = M_PI * brush_radius_3d_sq;
+
+    const Span<MLoopTri> looptris{BKE_mesh_runtime_looptri_ensure(&surface),
+                                  BKE_mesh_runtime_looptri_len(&surface)};
+
+    threading::EnumerableThreadSpecific<NewPointsData> new_points_per_thread;
+
+    const double time = PIL_check_seconds_timer();
+    const uint64_t time_as_int = *reinterpret_cast<const uint64_t *>(&time);
+    const uint32_t rng_base_seed = time_as_int ^ (time_as_int >> 32);
+
+    RandomNumberGenerator rng{rng_base_seed};
+
+    threading::parallel_for(looptri_indices.index_range(), 512, [&](const IndexRange range) {
+      RandomNumberGenerator looptri_rng{rng_base_seed + (uint32_t)range.start()};
+
+      for (const int looptri_index : looptri_indices.slice(range)) {
+        const MLoopTri &looptri = looptris[looptri_index];
+        const float3 &v0 = transform * float3(surface.mvert[surface.mloop[looptri.tri[0]].v].co);
+        const float3 &v1 = transform * float3(surface.mvert[surface.mloop[looptri.tri[1]].v].co);
+        const float3 &v2 = transform * float3(surface.mvert[surface.mloop[looptri.tri[2]].v].co);
+        const float looptri_area = area_tri_v3(v0, v1, v2);
+
+        float3 normal;
+        normal_tri_v3(normal, v0, v1, v2);
+
+        /* Use a different sampling strategy depending on whether the triangle is large or small
+         * compared to the brush size. When the triangle is small, points are distributed within
+         * the triangle directly. If the triangle is larger than the brush, distribute new points
+         * in a circle on the triangle plane. */
+        if (looptri_area < area_threshold) {
+          const int amount = this->float_to_int_amount(looptri_area * density, looptri_rng);
+
+          threading::parallel_for(IndexRange(amount), 512, [&](const IndexRange amount_range) {
+            RandomNumberGenerator point_rng{rng_base_seed + looptri_index * 1000 +
+                                            (uint32_t)amount_range.start()};
+            NewPointsData &new_points = new_points_per_thread.local();
+
+            for ([[maybe_unused]] const int i : amount_range) {
+              const float3 bary_coord = point_rng.get_barycentric_coordinates();
+              const float3 point_pos = attribute_math::mix3(bary_coord, v0, v1, v2);
+
+              if (math::distance(point_pos, brush_pos) > brush_radius_3d) {
+                continue;
+              }
+              if (minimum_distance > 0.0f &&
+                  this->is_too_close_to_existing_point(point_pos, minimum_distance)) {
+                continue;
+              }
+
+              new_points.bary_coords.append(bary_coord);
+              new_points.looptri_indices.append(looptri_index);
+              new_points.positions.append(point_pos);
+              new_points.normals.append(normal);
+            }
+          });
+        }
+        else {
+          float3 hit_pos_proj = brush_pos;
+          project_v3_plane(hit_pos_proj, normal, v0);
+          const float proj_distance_sq = math::distance_squared(hit_pos_proj, brush_pos);
+          const float brush_radius_factor_sq = 1.0f -
+                                               std::min(1.0f,
+                                                        proj_distance_sq / brush_radius_3d_sq);
+          const float radius_proj_sq = brush_radius_3d_sq * brush_radius_factor_sq;
+          const float radius_proj = std::sqrt(radius_proj_sq);
+          const float circle_area = M_PI * radius_proj_sq;
+
+          const int amount = this->float_to_int_amount(circle_area * density, rng);
+
+          const float3 axis_1 = math::normalize(v1 - v0) * radius_proj;
+          const float3 axis_2 = math::normalize(
+                                    math::cross(axis_1, math::cross(axis_1, v2 - v0))) *
+                                radius_proj;
+
+          threading::parallel_for(IndexRange(amount), 512, [&](const IndexRange amount_range) {
+            RandomNumberGenerator point_rng{rng_base_seed + looptri_index * 1000 +
+                                            (uint32_t)amount_range.start()};
+            NewPointsData &new_points = new_points_per_thread.local();
+
+            for ([[maybe_unused]] const int i : amount_range) {
+              const float r = std::sqrt(rng.get_float());
+              const float angle = rng.get_float() * 2 * M_PI;
+              const float x = r * std::cos(angle);
+              const float y = r * std::sin(angle);
+
+              const float3 point_pos = hit_pos_proj + axis_1 * x + axis_2 * y;
+
+              if (!isect_point_tri_prism_v3(point_pos, v0, v1, v2)) {
+                continue;
+              }
+              if (minimum_distance > 0.0f &&
+                  this->is_too_close_to_existing_point(point_pos, minimum_distance)) {
+                continue;
+              }
+
+              float3 bary_coord;
+              interp_weights_tri_v3(bary_coord, v0, v1, v2, point_pos);
+
+              new_points.bary_coords.append(bary_coord);
+              new_points.looptri_indices.append(looptri_index);
+              new_points.positions.append(point_pos);
+              new_points.normals.append(normal);
+            }
+          });
+        }
+      }
+    });
+
+    NewPointsData new_points;
+    for (const NewPointsData &local_new_points : new_points_per_thread) {
+      new_points.bary_coords.extend(local_new_points.bary_coords);
+      new_points.looptri_indices.extend(local_new_points.looptri_indices);
+      new_points.positions.extend(local_new_points.positions);
+      new_points.normals.extend(local_new_points.normals);
+    }
+    return new_points;
+  }
+
+  void eliminate_too_close_points(NewPointsData &points,
+                                  const CurvesGeometry &curves,
+                                  const float minimum_distance)
+  {
+    Array<bool> elimination_mask(points.positions.size(), false);
+
+    const int curves_added_previously = curves.curves_size() - old_kdtree_size_;
+    KDTree_3d *new_points_kdtree = this->kdtree_from_curve_roots_and_positions(
+        curves, IndexRange(old_kdtree_size_, curves_added_previously), points.positions);
+
+    Array<Vector<int>> points_in_range(points.positions.size());
+    threading::parallel_for(points.positions.index_range(), 256, [&](const IndexRange range) {
+      for (const int point_i : range) {
+        const float3 query_position = points.positions[point_i];
+
+        struct CallbackData {
+          int point_i;
+          Vector<int> &found_indices;
+          MutableSpan<bool> elimination_mask;
+        } callback_data = {point_i, points_in_range[point_i], elimination_mask};
+
+        BLI_kdtree_3d_range_search_cb(
+            new_points_kdtree,
+            query_position,
+            minimum_distance,
+            [](void *user_data, int index, const float *UNUSED(co), float UNUSED(dist_sq)) {
+              CallbackData &data = *static_cast<CallbackData *>(user_data);
+              if (index == data.point_i) {
+                /* Ignore self. */
+                return true;
+              }
+              if (index == ExistsAlreadyIndex) {
+                /* An already existing point is too close, so this new point will be eliminated. */
+                data.elimination_mask[data.point_i] = true;
+                return false;
+              }
+              data.found_indices.append(index);
+              return true;
+            },
+            &callback_data);
+      }
+    });
+
+    for (const int point_i : points.positions.index_range()) {
+      if (elimination_mask[point_i]) {
+        /* Point is eliminated already. */
+        continue;
+      }
+
+      for (const int other_point_i : points_in_range[point_i]) {
+        elimination_mask[other_point_i] = true;
+      }
+    }
+
+    BLI_kdtree_3d_free(new_points_kdtree);
+    for (int i = points.positions.size() - 1; i >= 0; i--) {
+      if (elimination_mask[i]) {
+        points.positions.remove_and_reorder(i);
+        points.bary_coords.remove_and_reorder(i);
+        points.looptri_indices.remove_and_reorder(i);
+        points.normals.remove_and_reorder(i);
+      }
+    }
+  }
+
+  void insert_new_curves(const NewPointsData &new_points, CurvesGeometry &curves)
+  {
+    const int tot_new_curves = new_points.positions.size();
+
+    const int points_per_curve = 8;
+    curves.resize(curves.points_size() + tot_new_curves * points_per_curve,
+                  curves.curves_size() + tot_new_curves);
+
+    MutableSpan<int> offsets = curves.offsets();
+    MutableSpan<float3> positions = curves.positions();
+
+    for (const int i : IndexRange(tot_new_curves)) {
+      const int curve_i = curves.curves_size() - tot_new_curves + i;
+      const int first_point_i = offsets[curve_i];
+      offsets[curve_i + 1] = offsets[curve_i] + points_per_curve;
+
+      const float3 root = new_points.positions[i];
+      const float3 tip = root + 0.1f * new_points.normals[i];
+
+      for (const int j : IndexRange(points_per_curve)) {
+        positions[first_point_i + j] = math::interpolate(
+            root, tip, j / (float)(points_per_curve - 1));
+      }
+    }
+  }
+};
+
 static std::unique_ptr<CurvesSculptStrokeOperation> start_brush_operation(bContext *C,
                                                                           wmOperator *UNUSED(op))
 {
@@ -212,6 +655,8 @@ static std::unique_ptr<CurvesSculptStrokeOperation> start_brush_operation(bConte
       return std::make_unique<MoveOperation>();
     case CURVES_SCULPT_TOOL_TEST2:
       return std::make_unique<DeleteOperation>();
+    case CURVES_SCULPT_TOOL_TEST3:
+      return std::make_unique<AddOperation>();
   }
   BLI_assert_unreachable();
   return {};

source/blender/makesdna/DNA_brush_enums.h

@@ -459,6 +459,7 @@ typedef enum eBrushUVSculptTool {
 typedef enum eBrushCurvesSculptTool {
   CURVES_SCULPT_TOOL_TEST1 = 0,
   CURVES_SCULPT_TOOL_TEST2 = 1,
+  CURVES_SCULPT_TOOL_TEST3 = 2,
 } eBrushCurvesSculptTool;
 
 /** When #BRUSH_ACCUMULATE is used */

source/blender/makesdna/DNA_scene_types.h

@@ -994,6 +994,9 @@ typedef struct Sculpt {
 
 typedef struct CurvesSculpt {
   Paint paint;
+  /** Minimum distance between newly added curves on a surface. */
+  float distance;
+  char _pad1[4];
 } CurvesSculpt;
 
 typedef struct UvSculpt {

source/blender/makesrna/intern/rna_brush.c

@@ -246,6 +246,7 @@ const EnumPropertyItem rna_enum_brush_gpencil_weight_types_items[] = {
 const EnumPropertyItem rna_enum_brush_curves_sculpt_tool_items[] = {
     {CURVES_SCULPT_TOOL_TEST1, "TEST1", ICON_NONE, "Test 1", ""},
     {CURVES_SCULPT_TOOL_TEST2, "TEST2", ICON_NONE, "Test 2", ""},
+    {CURVES_SCULPT_TOOL_TEST3, "TEST3", ICON_NONE, "Test 3", ""},
     {0, NULL, 0, NULL, NULL},
 };
 

source/blender/makesrna/intern/rna_sculpt_paint.c

@@ -1506,10 +1506,18 @@ static void rna_def_gpencil_sculpt(BlenderRNA *brna)
 static void rna_def_curves_sculpt(BlenderRNA *brna)
 {
   StructRNA *srna;
+  PropertyRNA *prop;
 
   srna = RNA_def_struct(brna, "CurvesSculpt", "Paint");
   RNA_def_struct_path_func(srna, "rna_CurvesSculpt_path");
   RNA_def_struct_ui_text(srna, "Curves Sculpt Paint", "");
+
+  prop = RNA_def_property(srna, "distance", PROP_FLOAT, PROP_DISTANCE);
+  RNA_def_property_range(prop, 0.0f, FLT_MAX);
+  RNA_def_property_ui_range(prop, 0.0f, FLT_MAX, 1, 6);
+  RNA_def_property_ui_text(
+      prop, "Distance", "Radius around curves roots in which no new curves can be added");
+  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
 }
 
 void RNA_def_sculpt_paint(BlenderRNA *brna)