Sculpt: Data oriented refactor for surface smooth brush

Part of #118145.

Pull Request: https://projects.blender.org/blender/blender/pulls/125394

source/blender/editors/sculpt_paint/CMakeLists.txt

@@ -151,6 +151,7 @@ set(SRC
   brushes/smooth.cc
   brushes/smooth_mask.cc
   brushes/snake_hook.cc
+  brushes/surface_smooth.cc
   brushes/thumb.cc
   brushes/topology_slide.cc
 

source/blender/editors/sculpt_paint/brushes/surface_smooth.cc

@@ -0,0 +1,450 @@
+/* SPDX-FileCopyrightText: 2024 Blender Authors
+ *
+ * SPDX-License-Identifier: GPL-2.0-or-later */
+
+#include "editors/sculpt_paint/brushes/types.hh"
+
+#include "DNA_brush_types.h"
+#include "DNA_mesh_types.h"
+#include "DNA_object_types.h"
+#include "DNA_scene_types.h"
+
+#include "BKE_mesh.hh"
+#include "BKE_paint.hh"
+#include "BKE_pbvh.hh"
+#include "BKE_subdiv_ccg.hh"
+
+#include "BLI_array.hh"
+#include "BLI_array_utils.hh"
+#include "BLI_enumerable_thread_specific.hh"
+#include "BLI_math_vector.hh"
+#include "BLI_task.hh"
+#include "BLI_virtual_array.hh"
+
+#include "editors/sculpt_paint/mesh_brush_common.hh"
+#include "editors/sculpt_paint/sculpt_intern.hh"
+
+namespace blender::ed::sculpt_paint {
+
+inline namespace surface_smooth_cc {
+
+struct LocalData {
+  Vector<float3> positions;
+  Vector<float> factors;
+  Vector<float> distances;
+  Vector<Vector<int>> vert_neighbors;
+  Vector<float3> laplacian_disp;
+  Vector<float3> average_positions;
+  Vector<float3> translations;
+};
+
+BLI_NOINLINE static void surface_smooth_laplacian_step(const Span<float3> positions,
+                                                       const Span<float3> orig_positions,
+                                                       const Span<float3> average_positions,
+                                                       const float alpha,
+                                                       MutableSpan<float3> laplacian_disp,
+                                                       MutableSpan<float3> translations)
+{
+  BLI_assert(positions.size() == orig_positions.size());
+  BLI_assert(positions.size() == average_positions.size());
+  BLI_assert(positions.size() == laplacian_disp.size());
+  BLI_assert(positions.size() == translations.size());
+
+  for (const int i : average_positions.index_range()) {
+    const float3 weighted_o = orig_positions[i] * alpha;
+    const float3 weighted_q = positions[i] * (1.0f - alpha);
+    const float3 d = weighted_o + weighted_q;
+    laplacian_disp[i] = average_positions[i] - d;
+    translations[i] = average_positions[i] - positions[i];
+  }
+}
+
+BLI_NOINLINE static void calc_displace_step(const Span<float3> laplacian_disp,
+                                            const Span<float3> average_laplacian_disp,
+                                            const float beta,
+                                            MutableSpan<float3> translations)
+{
+  BLI_assert(laplacian_disp.size() == average_laplacian_disp.size());
+  BLI_assert(laplacian_disp.size() == translations.size());
+
+  for (const int i : laplacian_disp.index_range()) {
+    float3 b_current_vert = average_laplacian_disp[i] * (1.0f - beta);
+    b_current_vert += laplacian_disp[i] * beta;
+    translations[i] = -b_current_vert;
+  }
+}
+
+BLI_NOINLINE static void clamp_factors(const MutableSpan<float> factors)
+{
+  for (float &factor : factors) {
+    factor = std::clamp(factor, 0.0f, 1.0f);
+  }
+}
+
+BLI_NOINLINE static void do_surface_smooth_brush_mesh(const Sculpt &sd,
+                                                      const Brush &brush,
+                                                      const Span<bke::pbvh::Node *> nodes,
+                                                      Object &object,
+                                                      const MutableSpan<float3> all_laplacian_disp)
+{
+  const SculptSession &ss = *object.sculpt;
+  const StrokeCache &cache = *ss.cache;
+  const float alpha = brush.surface_smooth_shape_preservation;
+  const float beta = brush.surface_smooth_current_vertex;
+
+  Mesh &mesh = *static_cast<Mesh *>(object.data);
+  const OffsetIndices faces = mesh.faces();
+  const Span<int> corner_verts = mesh.corner_verts();
+  const bke::AttributeAccessor attributes = mesh.attributes();
+  const VArraySpan hide_poly = *attributes.lookup<bool>(".hide_poly", bke::AttrDomain::Face);
+
+  const bke::pbvh::Tree &pbvh = *ss.pbvh;
+
+  const Span<float3> positions_eval = BKE_pbvh_get_vert_positions(pbvh);
+  const Span<float3> vert_normals = BKE_pbvh_get_vert_normals(pbvh);
+  MutableSpan<float3> positions_orig = mesh.vert_positions_for_write();
+
+  Array<int> node_offset_data;
+  const OffsetIndices node_offsets = create_node_vert_offsets(nodes, node_offset_data);
+  Array<float> all_factors(node_offsets.total_size());
+
+  threading::EnumerableThreadSpecific<LocalData> all_tls;
+  threading::parallel_for(nodes.index_range(), 1, [&](const IndexRange range) {
+    LocalData &tls = all_tls.local();
+    for (const int i : range) {
+      const Span<int> verts = bke::pbvh::node_unique_verts(*nodes[i]);
+
+      const MutableSpan<float> factors = all_factors.as_mutable_span().slice(node_offsets[i]);
+      fill_factor_from_hide_and_mask(mesh, verts, factors);
+      filter_region_clip_factors(ss, positions_eval, verts, factors);
+      if (brush.flag & BRUSH_FRONTFACE) {
+        calc_front_face(cache.view_normal, vert_normals, verts, factors);
+      }
+
+      tls.distances.reinitialize(verts.size());
+      const MutableSpan<float> distances = tls.distances;
+      calc_brush_distances(
+          ss, positions_eval, verts, eBrushFalloffShape(brush.falloff_shape), distances);
+      filter_distances_with_radius(cache.radius, distances, factors);
+      apply_hardness_to_distances(cache, distances);
+      calc_brush_strength_factors(cache, brush, distances, factors);
+
+      if (cache.automasking) {
+        auto_mask::calc_vert_factors(object, *cache.automasking, *nodes[i], verts, factors);
+      }
+
+      calc_brush_texture_factors(ss, brush, positions_eval, verts, factors);
+
+      scale_factors(factors, cache.bstrength);
+      clamp_factors(factors);
+    }
+  });
+
+  for ([[maybe_unused]] const int iteration : IndexRange(brush.surface_smooth_iterations)) {
+    threading::parallel_for(nodes.index_range(), 1, [&](const IndexRange range) {
+      LocalData &tls = all_tls.local();
+      for (const int i : range) {
+        const Span<int> verts = bke::pbvh::node_unique_verts(*nodes[i]);
+        const MutableSpan positions = gather_data_mesh(positions_eval, verts, tls.positions);
+        const OrigPositionData orig_data = orig_position_data_get_mesh(object, *nodes[i]);
+        const Span<float> factors = all_factors.as_span().slice(node_offsets[i]);
+
+        tls.vert_neighbors.reinitialize(verts.size());
+        calc_vert_neighbors(
+            faces, corner_verts, ss.vert_to_face_map, hide_poly, verts, tls.vert_neighbors);
+
+        tls.average_positions.reinitialize(verts.size());
+        const MutableSpan<float3> average_positions = tls.average_positions;
+        smooth::neighbor_data_average_mesh(positions_eval, tls.vert_neighbors, average_positions);
+
+        tls.laplacian_disp.reinitialize(verts.size());
+        const MutableSpan<float3> laplacian_disp = tls.laplacian_disp;
+        tls.translations.reinitialize(verts.size());
+        const MutableSpan<float3> translations = tls.translations;
+        surface_smooth_laplacian_step(positions,
+                                      orig_data.positions,
+                                      average_positions,
+                                      alpha,
+                                      laplacian_disp,
+                                      translations);
+        scale_translations(translations, factors);
+
+        scatter_data_mesh(laplacian_disp.as_span(), verts, all_laplacian_disp);
+
+        write_translations(sd, object, positions_eval, verts, translations, positions_orig);
+      }
+    });
+
+    threading::parallel_for(nodes.index_range(), 1, [&](const IndexRange range) {
+      LocalData &tls = all_tls.local();
+      for (const int i : range) {
+        const Span<int> verts = bke::pbvh::node_unique_verts(*nodes[i]);
+        const Span<float> factors = all_factors.as_span().slice(node_offsets[i]);
+
+        tls.laplacian_disp.reinitialize(verts.size());
+        const MutableSpan<float3> laplacian_disp = tls.laplacian_disp;
+        gather_data_mesh(all_laplacian_disp.as_span(), verts, laplacian_disp);
+
+        tls.vert_neighbors.reinitialize(verts.size());
+        calc_vert_neighbors(
+            faces, corner_verts, ss.vert_to_face_map, hide_poly, verts, tls.vert_neighbors);
+
+        tls.average_positions.reinitialize(verts.size());
+        const MutableSpan<float3> average_laplacian_disps = tls.average_positions;
+        smooth::neighbor_data_average_mesh(
+            all_laplacian_disp.as_span(), tls.vert_neighbors, average_laplacian_disps);
+
+        tls.translations.reinitialize(verts.size());
+        const MutableSpan<float3> translations = tls.translations;
+        calc_displace_step(laplacian_disp, average_laplacian_disps, beta, translations);
+        scale_translations(translations, factors);
+
+        write_translations(sd, object, positions_eval, verts, translations, positions_orig);
+      }
+    });
+  }
+}
+
+BLI_NOINLINE static void do_surface_smooth_brush_grids(
+    const Sculpt &sd,
+    const Brush &brush,
+    const Span<bke::pbvh::Node *> nodes,
+    Object &object,
+    const MutableSpan<float3> all_laplacian_disp)
+{
+  const SculptSession &ss = *object.sculpt;
+  const StrokeCache &cache = *ss.cache;
+  const float alpha = brush.surface_smooth_shape_preservation;
+  const float beta = brush.surface_smooth_current_vertex;
+
+  SubdivCCG &subdiv_ccg = *ss.subdiv_ccg;
+
+  Array<int> node_offset_data;
+  const OffsetIndices node_offsets = create_node_vert_offsets(
+      nodes, BKE_subdiv_ccg_key_top_level(subdiv_ccg), node_offset_data);
+  Array<float> all_factors(node_offsets.total_size());
+
+  threading::EnumerableThreadSpecific<LocalData> all_tls;
+  threading::parallel_for(nodes.index_range(), 1, [&](const IndexRange range) {
+    LocalData &tls = all_tls.local();
+    for (const int i : range) {
+      const Span<int> grids = bke::pbvh::node_grid_indices(*nodes[i]);
+      const MutableSpan positions = gather_grids_positions(subdiv_ccg, grids, tls.positions);
+
+      const MutableSpan<float> factors = all_factors.as_mutable_span().slice(node_offsets[i]);
+      fill_factor_from_hide_and_mask(subdiv_ccg, grids, factors);
+      filter_region_clip_factors(ss, positions, factors);
+      if (brush.flag & BRUSH_FRONTFACE) {
+        calc_front_face(cache.view_normal, subdiv_ccg, grids, factors);
+      }
+
+      tls.distances.reinitialize(positions.size());
+      const MutableSpan<float> distances = tls.distances;
+      calc_brush_distances(ss, positions, eBrushFalloffShape(brush.falloff_shape), distances);
+      filter_distances_with_radius(cache.radius, distances, factors);
+      apply_hardness_to_distances(cache, distances);
+      calc_brush_strength_factors(cache, brush, distances, factors);
+
+      if (cache.automasking) {
+        auto_mask::calc_grids_factors(object, *cache.automasking, *nodes[i], grids, factors);
+      }
+
+      calc_brush_texture_factors(ss, brush, positions, factors);
+
+      scale_factors(factors, cache.bstrength);
+      clamp_factors(factors);
+    }
+  });
+
+  for ([[maybe_unused]] const int iteration : IndexRange(brush.surface_smooth_iterations)) {
+    threading::parallel_for(nodes.index_range(), 1, [&](const IndexRange range) {
+      LocalData &tls = all_tls.local();
+      for (const int i : range) {
+        const Span<int> grids = bke::pbvh::node_grid_indices(*nodes[i]);
+        const MutableSpan positions = gather_grids_positions(subdiv_ccg, grids, tls.positions);
+        const OrigPositionData orig_data = orig_position_data_get_mesh(object, *nodes[i]);
+        const Span<float> factors = all_factors.as_span().slice(node_offsets[i]);
+
+        tls.average_positions.reinitialize(positions.size());
+        const MutableSpan<float3> average_positions = tls.average_positions;
+        smooth::neighbor_position_average_grids(subdiv_ccg, grids, average_positions);
+
+        tls.laplacian_disp.reinitialize(positions.size());
+        const MutableSpan<float3> laplacian_disp = tls.laplacian_disp;
+        tls.translations.reinitialize(positions.size());
+        const MutableSpan<float3> translations = tls.translations;
+        surface_smooth_laplacian_step(positions,
+                                      orig_data.positions,
+                                      average_positions,
+                                      alpha,
+                                      laplacian_disp,
+                                      translations);
+        scale_translations(translations, factors);
+
+        scatter_data_grids(subdiv_ccg, laplacian_disp.as_span(), grids, all_laplacian_disp);
+
+        clip_and_lock_translations(sd, ss, positions, translations);
+        apply_translations(translations, grids, subdiv_ccg);
+      }
+    });
+
+    threading::parallel_for(nodes.index_range(), 1, [&](const IndexRange range) {
+      LocalData &tls = all_tls.local();
+      for (const int i : range) {
+        const Span<int> grids = bke::pbvh::node_grid_indices(*nodes[i]);
+        const MutableSpan positions = gather_grids_positions(subdiv_ccg, grids, tls.positions);
+        const Span<float> factors = all_factors.as_span().slice(node_offsets[i]);
+
+        tls.laplacian_disp.reinitialize(positions.size());
+        const MutableSpan<float3> laplacian_disp = tls.laplacian_disp;
+        gather_data_grids(subdiv_ccg, all_laplacian_disp.as_span(), grids, laplacian_disp);
+
+        tls.average_positions.reinitialize(positions.size());
+        const MutableSpan<float3> average_laplacian_disps = tls.average_positions;
+        smooth::average_data_grids(
+            subdiv_ccg, all_laplacian_disp.as_span(), grids, average_laplacian_disps);
+
+        tls.translations.reinitialize(positions.size());
+        const MutableSpan<float3> translations = tls.translations;
+        calc_displace_step(laplacian_disp, average_laplacian_disps, beta, translations);
+        scale_translations(translations, factors);
+
+        clip_and_lock_translations(sd, ss, positions, translations);
+        apply_translations(translations, grids, subdiv_ccg);
+      }
+    });
+  }
+}
+
+BLI_NOINLINE static void do_surface_smooth_brush_bmesh(
+    const Sculpt &sd,
+    const Brush &brush,
+    const Span<bke::pbvh::Node *> nodes,
+    Object &object,
+    const MutableSpan<float3> all_laplacian_disp)
+{
+  const SculptSession &ss = *object.sculpt;
+  const StrokeCache &cache = *ss.cache;
+  const float alpha = brush.surface_smooth_shape_preservation;
+  const float beta = brush.surface_smooth_current_vertex;
+
+  Array<int> node_offset_data;
+  const OffsetIndices node_offsets = create_node_vert_offsets_bmesh(nodes, node_offset_data);
+  Array<float> all_factors(node_offsets.total_size());
+
+  threading::EnumerableThreadSpecific<LocalData> all_tls;
+  threading::parallel_for(nodes.index_range(), 1, [&](const IndexRange range) {
+    LocalData &tls = all_tls.local();
+    for (const int i : range) {
+      const Set<BMVert *, 0> &verts = BKE_pbvh_bmesh_node_unique_verts(nodes[i]);
+      const MutableSpan positions = gather_bmesh_positions(verts, tls.positions);
+
+      const MutableSpan<float> factors = all_factors.as_mutable_span().slice(node_offsets[i]);
+      fill_factor_from_hide_and_mask(*ss.bm, verts, factors);
+      filter_region_clip_factors(ss, positions, factors);
+      if (brush.flag & BRUSH_FRONTFACE) {
+        calc_front_face(cache.view_normal, verts, factors);
+      }
+
+      tls.distances.reinitialize(positions.size());
+      const MutableSpan<float> distances = tls.distances;
+      calc_brush_distances(ss, positions, eBrushFalloffShape(brush.falloff_shape), distances);
+      filter_distances_with_radius(cache.radius, distances, factors);
+      apply_hardness_to_distances(cache, distances);
+      calc_brush_strength_factors(cache, brush, distances, factors);
+
+      if (cache.automasking) {
+        auto_mask::calc_vert_factors(object, *cache.automasking, *nodes[i], verts, factors);
+      }
+
+      calc_brush_texture_factors(ss, brush, positions, factors);
+
+      scale_factors(factors, cache.bstrength);
+      clamp_factors(factors);
+    }
+  });
+
+  for ([[maybe_unused]] const int iteration : IndexRange(brush.surface_smooth_iterations)) {
+    threading::parallel_for(nodes.index_range(), 1, [&](const IndexRange range) {
+      LocalData &tls = all_tls.local();
+      for (const int i : range) {
+        const Set<BMVert *, 0> &verts = BKE_pbvh_bmesh_node_unique_verts(nodes[i]);
+        const MutableSpan positions = gather_bmesh_positions(verts, tls.positions);
+        Array<float3> orig_positions(verts.size());
+        Array<float3> orig_normals(verts.size());
+        orig_position_data_gather_bmesh(*ss.bm_log, verts, orig_positions, orig_normals);
+        const Span<float> factors = all_factors.as_span().slice(node_offsets[i]);
+
+        tls.average_positions.reinitialize(positions.size());
+        const MutableSpan<float3> average_positions = tls.average_positions;
+        smooth::neighbor_position_average_bmesh(verts, average_positions);
+
+        tls.laplacian_disp.reinitialize(positions.size());
+        const MutableSpan<float3> laplacian_disp = tls.laplacian_disp;
+        tls.translations.reinitialize(positions.size());
+        const MutableSpan<float3> translations = tls.translations;
+        surface_smooth_laplacian_step(
+            positions, orig_positions, average_positions, alpha, laplacian_disp, translations);
+        scale_translations(translations, factors);
+
+        scatter_data_vert_bmesh(laplacian_disp.as_span(), verts, all_laplacian_disp);
+
+        clip_and_lock_translations(sd, ss, positions, translations);
+        apply_translations(translations, verts);
+      }
+    });
+
+    threading::parallel_for(nodes.index_range(), 1, [&](const IndexRange range) {
+      LocalData &tls = all_tls.local();
+      for (const int i : range) {
+        const Set<BMVert *, 0> &verts = BKE_pbvh_bmesh_node_unique_verts(nodes[i]);
+        const MutableSpan positions = gather_bmesh_positions(verts, tls.positions);
+        const Span<float> factors = all_factors.as_span().slice(node_offsets[i]);
+
+        tls.laplacian_disp.reinitialize(positions.size());
+        const MutableSpan<float3> laplacian_disp = tls.laplacian_disp;
+        gather_data_vert_bmesh(all_laplacian_disp.as_span(), verts, laplacian_disp);
+
+        tls.average_positions.reinitialize(positions.size());
+        const MutableSpan<float3> average_laplacian_disps = tls.average_positions;
+        smooth::average_data_bmesh(all_laplacian_disp.as_span(), verts, average_laplacian_disps);
+
+        tls.translations.reinitialize(positions.size());
+        const MutableSpan<float3> translations = tls.translations;
+        calc_displace_step(laplacian_disp, average_laplacian_disps, beta, translations);
+        scale_translations(translations, factors);
+
+        clip_and_lock_translations(sd, ss, positions, translations);
+        apply_translations(translations, verts);
+      }
+    });
+  }
+}
+
+}  // namespace surface_smooth_cc
+
+void do_surface_smooth_brush(const Sculpt &sd, Object &object, const Span<bke::pbvh::Node *> nodes)
+{
+  SculptSession &ss = *object.sculpt;
+  const Brush &brush = *BKE_paint_brush_for_read(&sd.paint);
+
+  switch (ss.pbvh->type()) {
+    case bke::pbvh::Type::Mesh:
+      do_surface_smooth_brush_mesh(
+          sd, brush, nodes, object, ss.cache->surface_smooth_laplacian_disp);
+      break;
+    case bke::pbvh::Type::Grids: {
+      do_surface_smooth_brush_grids(
+          sd, brush, nodes, object, ss.cache->surface_smooth_laplacian_disp);
+      break;
+    }
+    case bke::pbvh::Type::BMesh: {
+      BM_mesh_elem_index_ensure(ss.bm, BM_VERT);
+      do_surface_smooth_brush_bmesh(
+          sd, brush, nodes, object, ss.cache->surface_smooth_laplacian_disp);
+      break;
+    }
+  }
+}
+
+}  // namespace blender::ed::sculpt_paint

source/blender/editors/sculpt_paint/brushes/types.hh

@@ -67,6 +67,7 @@ void do_smooth_mask_brush(const Sculpt &sd,
                           Span<bke::pbvh::Node *> nodes,
                           float brush_strength);
 void do_snake_hook_brush(const Sculpt &sd, Object &object, Span<bke::pbvh::Node *> nodes);
+void do_surface_smooth_brush(const Sculpt &sd, Object &object, Span<bke::pbvh::Node *> nodes);
 void do_thumb_brush(const Sculpt &sd, Object &object, Span<bke::pbvh::Node *> nodes);
 void do_topology_slide_brush(const Sculpt &sd, Object &object, Span<bke::pbvh::Node *> nodes);
 void do_topology_relax_brush(const Sculpt &sd, Object &object, Span<bke::pbvh::Node *> nodes);

source/blender/editors/sculpt_paint/sculpt.cc

@@ -3875,7 +3875,7 @@ static void do_brush_action(const Scene &scene,
         }
       }
       else if (brush.smooth_deform_type == BRUSH_SMOOTH_DEFORM_SURFACE) {
-        smooth::do_surface_smooth_brush(sd, ob, nodes);
+        do_surface_smooth_brush(sd, ob, nodes);
       }
       break;
     case SCULPT_TOOL_CREASE:
@@ -5963,11 +5963,7 @@ static void sculpt_stroke_update_step(bContext *C,
    *
    * For some brushes, flushing is done in the brush code itself.
    */
-  if ((ELEM(brush.sculpt_tool,
-            SCULPT_TOOL_BOUNDARY,
-            SCULPT_TOOL_CLOTH,
-            SCULPT_TOOL_POSE,
-            SCULPT_TOOL_SMOOTH) ||
+  if ((ELEM(brush.sculpt_tool, SCULPT_TOOL_BOUNDARY, SCULPT_TOOL_CLOTH, SCULPT_TOOL_POSE) ||
        ss.pbvh->type() != bke::pbvh::Type::Mesh))
   {
     if (ss.deform_modifiers_active) {

source/blender/editors/sculpt_paint/sculpt_intern.hh

@@ -1695,6 +1695,15 @@ void neighbor_data_average_mesh_check_loose(Span<T> src,
                                             Span<Vector<int>> vert_neighbors,
                                             MutableSpan<T> dst);
 
+template<typename T>
+void average_data_grids(const SubdivCCG &subdiv_ccg,
+                        Span<T> src,
+                        Span<int> grids,
+                        MutableSpan<T> dst);
+
+template<typename T>
+void average_data_bmesh(Span<T> src, const Set<BMVert *, 0> &verts, MutableSpan<T> dst);
+
 /* Surface Smooth Brush. */
 
 void surface_smooth_laplacian_step(SculptSession &ss,
@@ -1710,7 +1719,6 @@ void surface_smooth_displace_step(SculptSession &ss,
                                   PBVHVertRef vertex,
                                   float beta,
                                   float fade);
-void do_surface_smooth_brush(const Sculpt &sd, Object &ob, Span<blender::bke::pbvh::Node *> nodes);
 
 /* Slide/Relax */
 void relax_vertex(SculptSession &ss,

source/blender/editors/sculpt_paint/sculpt_smooth.cc

@@ -191,6 +191,73 @@ void neighbor_position_average_interior_grids(const OffsetIndices<int> faces,
   }
 }
 
+template<typename T>
+void average_data_grids(const SubdivCCG &subdiv_ccg,
+                        const Span<T> src,
+                        const Span<int> grids,
+                        const MutableSpan<T> dst)
+{
+  const CCGKey key = BKE_subdiv_ccg_key_top_level(subdiv_ccg);
+
+  BLI_assert(grids.size() * key.grid_area == src.size());
+
+  for (const int i : grids.index_range()) {
+    const int grid = grids[i];
+    const int node_verts_start = i * key.grid_area;
+
+    /* TODO: This loop could be optimized in the future by skipping unnecessary logic for
+     * non-boundary grid vertices. */
+    for (const int y : IndexRange(key.grid_size)) {
+      for (const int x : IndexRange(key.grid_size)) {
+        const int offset = CCG_grid_xy_to_index(key.grid_size, x, y);
+        const int node_vert_index = node_verts_start + offset;
+
+        SubdivCCGCoord coord{};
+        coord.grid_index = grid;
+        coord.x = x;
+        coord.y = y;
+
+        SubdivCCGNeighbors neighbors;
+        BKE_subdiv_ccg_neighbor_coords_get(subdiv_ccg, coord, false, neighbors);
+
+        T sum{};
+        for (const SubdivCCGCoord neighbor : neighbors.coords) {
+          const int index = neighbor.grid_index * key.grid_area +
+                            CCG_grid_xy_to_index(key.grid_size, neighbor.x, neighbor.y);
+          sum += src[index];
+        }
+        dst[node_vert_index] = sum / neighbors.coords.size();
+      }
+    }
+  }
+}
+
+template<typename T>
+void average_data_bmesh(const Span<T> src, const Set<BMVert *, 0> &verts, const MutableSpan<T> dst)
+{
+  Vector<BMVert *, 64> neighbor_data;
+
+  int i = 0;
+  for (BMVert *vert : verts) {
+    T sum{};
+    neighbor_data.clear();
+    const Span<BMVert *> neighbors = vert_neighbors_get_bmesh(*vert, neighbor_data);
+    for (const BMVert *neighbor : neighbors) {
+      sum += src[BM_elem_index_get(neighbor)];
+    }
+    dst[i] = sum / neighbors.size();
+    i++;
+  }
+}
+
+template void average_data_grids<float3>(const SubdivCCG &,
+                                         Span<float3>,
+                                         Span<int>,
+                                         MutableSpan<float3>);
+template void average_data_bmesh<float3>(Span<float3> src,
+                                         const Set<BMVert *, 0> &,
+                                         MutableSpan<float3>);
+
 static float3 average_positions(const Span<const BMVert *> verts)
 {
   const float factor = math::rcp(float(verts.size()));
@@ -457,108 +524,4 @@ void surface_smooth_displace_step(SculptSession &ss,
   }
 }
 
-static void do_surface_smooth_brush_laplacian_task(Object &ob,
-                                                   const Brush &brush,
-                                                   bke::pbvh::Node *node)
-{
-  SculptSession &ss = *ob.sculpt;
-  const float bstrength = ss.cache->bstrength;
-  float alpha = brush.surface_smooth_shape_preservation;
-
-  PBVHVertexIter vd;
-
-  SculptBrushTest test;
-  SculptBrushTestFn sculpt_brush_test_sq_fn = SCULPT_brush_test_init_with_falloff_shape(
-      ss, test, brush.falloff_shape);
-  const int thread_id = BLI_task_parallel_thread_id(nullptr);
-
-  SculptOrigVertData orig_data = SCULPT_orig_vert_data_init(ob, *node, undo::Type::Position);
-  auto_mask::NodeData automask_data = auto_mask::node_begin(
-      ob, ss.cache->automasking.get(), *node);
-
-  BKE_pbvh_vertex_iter_begin (*ss.pbvh, node, vd, PBVH_ITER_UNIQUE) {
-    SCULPT_orig_vert_data_update(orig_data, vd);
-    if (!sculpt_brush_test_sq_fn(test, vd.co)) {
-      continue;
-    }
-
-    auto_mask::node_update(automask_data, vd);
-
-    const float fade = bstrength * SCULPT_brush_strength_factor(ss,
-                                                                brush,
-                                                                vd.co,
-                                                                sqrtf(test.dist),
-                                                                vd.no,
-                                                                vd.fno,
-                                                                vd.mask,
-                                                                vd.vertex,
-                                                                thread_id,
-                                                                &automask_data);
-
-    float disp[3];
-    surface_smooth_laplacian_step(
-        ss, disp, vd.co, ss.cache->surface_smooth_laplacian_disp, vd.vertex, orig_data.co, alpha);
-    madd_v3_v3fl(vd.co, disp, clamp_f(fade, 0.0f, 1.0f));
-  }
-  BKE_pbvh_vertex_iter_end;
-}
-
-static void do_surface_smooth_brush_displace_task(Object &ob,
-                                                  const Brush &brush,
-                                                  bke::pbvh::Node *node)
-{
-  SculptSession &ss = *ob.sculpt;
-  const float bstrength = ss.cache->bstrength;
-  const float beta = brush.surface_smooth_current_vertex;
-
-  PBVHVertexIter vd;
-
-  SculptBrushTest test;
-  SculptBrushTestFn sculpt_brush_test_sq_fn = SCULPT_brush_test_init_with_falloff_shape(
-      ss, test, brush.falloff_shape);
-  const int thread_id = BLI_task_parallel_thread_id(nullptr);
-  auto_mask::NodeData automask_data = auto_mask::node_begin(
-      ob, ss.cache->automasking.get(), *node);
-
-  BKE_pbvh_vertex_iter_begin (*ss.pbvh, node, vd, PBVH_ITER_UNIQUE) {
-    if (!sculpt_brush_test_sq_fn(test, vd.co)) {
-      continue;
-    }
-
-    auto_mask::node_update(automask_data, vd);
-
-    const float fade = bstrength * SCULPT_brush_strength_factor(ss,
-                                                                brush,
-                                                                vd.co,
-                                                                sqrtf(test.dist),
-                                                                vd.no,
-                                                                vd.fno,
-                                                                vd.mask,
-                                                                vd.vertex,
-                                                                thread_id,
-                                                                &automask_data);
-    surface_smooth_displace_step(
-        ss, vd.co, ss.cache->surface_smooth_laplacian_disp, vd.vertex, beta, fade);
-  }
-  BKE_pbvh_vertex_iter_end;
-}
-
-void do_surface_smooth_brush(const Sculpt &sd, Object &ob, Span<bke::pbvh::Node *> nodes)
-{
-  const Brush &brush = *BKE_paint_brush_for_read(&sd.paint);
-
-  for (int i = 0; i < brush.surface_smooth_iterations; i++) {
-    threading::parallel_for(nodes.index_range(), 1, [&](const IndexRange range) {
-      for (const int i : range) {
-        do_surface_smooth_brush_laplacian_task(ob, brush, nodes[i]);
-      }
-    });
-    threading::parallel_for(nodes.index_range(), 1, [&](const IndexRange range) {
-      for (const int i : range) {
-        do_surface_smooth_brush_displace_task(ob, brush, nodes[i]);
-      }
-    });
-  }
-}
-
 }  // namespace blender::ed::sculpt_paint::smooth