class AddOperation : public CurvesSculptStrokeOperation {
 private:
  KDTree_3d *curve_roots_kdtree_ = nullptr;

 public:
  ~AddOperation()
  {
    if (curve_roots_kdtree_ != nullptr) {
      BLI_kdtree_3d_free(curve_roots_kdtree_);
    }
  }

  /* Coordinate space naming convention:
   * 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.
   */

  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 float4x4 curves_to_world_mat = object.obmat;
    const float4x4 world_to_curves_mat = curves_to_world_mat.inverted();

    const Object &surface_ob = *curves_id.surface;
    const Mesh &surface = *static_cast<const Mesh *>(surface_ob.data);
    const float4x4 surface_to_world_mat = surface_ob.obmat;
    const float4x4 world_to_surface_mat = surface_to_world_mat.inverted();
    const float4x4 surface_to_curves_mat = world_to_curves_mat * surface_to_world_mat;

    CurvesSculpt &curves_sculpt = *scene.toolsettings->curves_sculpt;
    Brush &brush = *BKE_paint_brush(&curves_sculpt.paint);
    const float brush_radius_re = BKE_brush_size_get(&scene, &brush);
    const float2 brush_pos_re = stroke_extension.mouse_position;

    const bool use_front_face = brush.flag & BRUSH_FRONTFACE;
    const eBrushFalloffShape falloff_shape = static_cast<eBrushFalloffShape>(brush.falloff_shape);
    const int add_amount = std::max(0, brush.curves_sculpt_settings->add_amount);
    const bool interpolate_length = curves_sculpt.flag & CURVES_SCULPT_FLAG_INTERPOLATE_LENGTH;
    const bool interpolate_shape = curves_sculpt.flag & CURVES_SCULPT_FLAG_INTERPOLATE_SHAPE;
    const float new_curve_length = curves_sculpt.curve_length;
    const int points_per_curve = 8;

    if (add_amount == 0) {
      return;
    }

    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 Span<MLoopTri> surface_looptris{BKE_mesh_runtime_looptri_ensure(&surface),
                                          BKE_mesh_runtime_looptri_len(&surface)};

    Vector<float3> added_positions_cu;
    Vector<float3> added_bary_coords;
    Vector<int> added_looptri_indices;
    if (add_amount == 1) {
      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;

      added_positions_cu.append(brush_pos_cu);
      added_bary_coords.append(bary_coords);
      added_looptri_indices.append(looptri_index);
    }
    else if (falloff_shape == PAINT_FALLOFF_SHAPE_TUBE) {
      RandomNumberGenerator rng{(uint32_t)(PIL_check_seconds_timer() * 1000000.0f)};
      for (const int i : IndexRange(add_amount)) {
        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;

        added_positions_cu.append(pos_cu);
        added_bary_coords.append(bary_coords);
        added_looptri_indices.append(looptri_index);
      }
    }
    else if (falloff_shape == PAINT_FALLOFF_SHAPE_SPHERE) {
      /* Find 3d brush center and scatter n points around it. */
    }
    else {
      BLI_assert_unreachable();
    }

    if (added_bary_coords.is_empty()) {
      return;
    }

    if (interpolate_length || interpolate_shape) {
      if (curve_roots_kdtree_ == nullptr) {
        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(curve_roots_kdtree_, curve_i, root_pos_cu);
        }
        BLI_kdtree_3d_balance(curve_roots_kdtree_);
      }
    }

    const int tot_added_curves = added_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);

    MutableSpan<int> offsets = curves.offsets();
    MutableSpan<float3> positions_cu = curves.positions();

    for (const int i : IndexRange(tot_added_curves)) {
      const int curve_i = curves.curves_size() - tot_added_curves + i;
      const int first_point_i = offsets[curve_i];
      offsets[curve_i + 1] = offsets[curve_i] + points_per_curve;

      const float3 root_cu = added_positions_cu[i];
      const int looptri_index = added_looptri_indices[i];
      const MLoopTri &looptri = surface_looptris[looptri_index];

      if (interpolate_length || interpolate_shape) {
        KDTreeNearest_3d nearest;
        nearest.index = -1;
        nearest.dist = FLT_MAX;
        BLI_kdtree_3d_find_nearest(curve_roots_kdtree_, root_cu, &nearest);
        if (nearest.index >= 0) {
          const int nearest_curve_i = nearest.index;
          const IndexRange nearest_curve_points_range = curves.range_for_curve(nearest_curve_i);
          const float3 &nearest_curve_root_cu = positions_cu[nearest_curve_points_range[0]];
          const float3 offset_cu = root_cu - nearest_curve_root_cu;
          BLI_assert(nearest_curve_points_range.size() == points_per_curve);
          for (const int j : IndexRange(points_per_curve)) {
            positions_cu[first_point_i + j] = positions_cu[nearest_curve_points_range[j]] +
                                              offset_cu;
          }
          continue;
        }
      }
      const float3 v0_cu = surface_to_curves_mat *
                           float3(surface.mvert[surface.mloop[looptri.tri[0]].v].co);
      const float3 v1_cu = surface_to_curves_mat *
                           float3(surface.mvert[surface.mloop[looptri.tri[1]].v].co);
      const float3 v2_cu = surface_to_curves_mat *
                           float3(surface.mvert[surface.mloop[looptri.tri[2]].v].co);

      float3 normal_cu;
      normal_tri_v3(normal_cu, v0_cu, v1_cu, v2_cu);

      const float3 tip_cu = root_cu + normal_cu * new_curve_length;

      for (const int j : IndexRange(points_per_curve)) {
        positions_cu[first_point_i + j] = math::interpolate(
            root_cu, tip_cu, j / (float)(points_per_curve - 1));
      }
    }

    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;
  }
};