diff --git a/source/blender/editors/curves/intern/curves_ops.cc b/source/blender/editors/curves/intern/curves_ops.cc
index 1d2b1264477..5673c6876ed 100644
--- a/source/blender/editors/curves/intern/curves_ops.cc
+++ b/source/blender/editors/curves/intern/curves_ops.cc
@@ -635,7 +635,7 @@ static void snap_curves_to_surface_exec_object(Object &curves_ob,
             continue;
           }
 
-          const MLoopTri &looptri = *lookup_result.looptri;
+          const MLoopTri &looptri = surface_looptris[lookup_result.looptri_index];
           const float3 &bary_coords = lookup_result.bary_weights;
 
           const float3 &p0_su = verts[loops[looptri.tri[0]].v].co;
diff --git a/source/blender/editors/sculpt_paint/curves_sculpt_add.cc b/source/blender/editors/sculpt_paint/curves_sculpt_add.cc
index 0d2c2d3f0c9..fa9cf3d3a04 100644
--- a/source/blender/editors/sculpt_paint/curves_sculpt_add.cc
+++ b/source/blender/editors/sculpt_paint/curves_sculpt_add.cc
@@ -228,6 +228,7 @@ struct AddOperationExecutor {
     add_inputs.fallback_curve_length = brush_settings_->curve_length;
     add_inputs.fallback_point_count = std::max(2, brush_settings_->points_per_curve);
     add_inputs.transforms = &transforms_;
+    add_inputs.surface_looptris = surface_looptris_orig;
     add_inputs.reverse_uv_sampler = &reverse_uv_sampler;
     add_inputs.surface = &surface_orig;
     add_inputs.corner_normals_su = corner_normals_su;
diff --git a/source/blender/editors/sculpt_paint/curves_sculpt_density.cc b/source/blender/editors/sculpt_paint/curves_sculpt_density.cc
index 78e2d55e6b9..78eea553737 100644
--- a/source/blender/editors/sculpt_paint/curves_sculpt_density.cc
+++ b/source/blender/editors/sculpt_paint/curves_sculpt_density.cc
@@ -280,6 +280,7 @@ struct DensityAddOperationExecutor {
     add_inputs.transforms = &transforms_;
     add_inputs.surface = surface_orig_;
     add_inputs.corner_normals_su = corner_normals_su;
+    add_inputs.surface_looptris = surface_looptris_orig;
     add_inputs.reverse_uv_sampler = &reverse_uv_sampler;
     add_inputs.old_roots_kdtree = self_->original_curve_roots_kdtree_;
 
diff --git a/source/blender/editors/sculpt_paint/curves_sculpt_slide.cc b/source/blender/editors/sculpt_paint/curves_sculpt_slide.cc
index ae89bc1c58b..89110f72e4b 100644
--- a/source/blender/editors/sculpt_paint/curves_sculpt_slide.cc
+++ b/source/blender/editors/sculpt_paint/curves_sculpt_slide.cc
@@ -292,8 +292,9 @@ struct SlideOperationExecutor {
       /* Compute the normal at the initial surface position. */
       const float3 normal_cu = math::normalize(
           transforms_.surface_to_curves_normal *
-          geometry::compute_surface_point_normal(
-              *result.looptri, result.bary_weights, corner_normals_orig_su_));
+          geometry::compute_surface_point_normal(surface_looptris_orig_[result.looptri_index],
+                                                 result.bary_weights,
+                                                 corner_normals_orig_su_));
 
       r_curves_to_slide.append({curve_i, radius_falloff, normal_cu});
     }
@@ -389,7 +390,7 @@ struct SlideOperationExecutor {
           found_invalid_uv_mapping_.store(true);
           continue;
         }
-        const MLoopTri &looptri_orig = *result.looptri;
+        const MLoopTri &looptri_orig = surface_looptris_orig_[result.looptri_index];
         const float3 &bary_weights_orig = result.bary_weights;
 
         /* Gather old and new surface normal. */
@@ -397,7 +398,7 @@ struct SlideOperationExecutor {
         const float3 new_normal_cu = math::normalize(
             transforms_.surface_to_curves_normal *
             geometry::compute_surface_point_normal(
-                *result.looptri, result.bary_weights, corner_normals_orig_su_));
+                looptri_orig, result.bary_weights, corner_normals_orig_su_));
 
         /* Gather old and new surface position. */
         const float3 old_first_pos_orig_cu = self_->initial_positions_cu_[first_point_i];
diff --git a/source/blender/geometry/GEO_add_curves_on_mesh.hh b/source/blender/geometry/GEO_add_curves_on_mesh.hh
index 2d2ba89a18f..d0f72deeb46 100644
--- a/source/blender/geometry/GEO_add_curves_on_mesh.hh
+++ b/source/blender/geometry/GEO_add_curves_on_mesh.hh
@@ -30,6 +30,7 @@ struct AddCurvesOnMeshInputs {
 
   /** Information about the surface that the new curves are attached to. */
   const Mesh *surface = nullptr;
+  Span<MLoopTri> surface_looptris;
   const ReverseUVSampler *reverse_uv_sampler = nullptr;
   Span<float3> corner_normals_su;
 
diff --git a/source/blender/geometry/GEO_reverse_uv_sampler.hh b/source/blender/geometry/GEO_reverse_uv_sampler.hh
index ee91e0b0731..62c6382c9d8 100644
--- a/source/blender/geometry/GEO_reverse_uv_sampler.hh
+++ b/source/blender/geometry/GEO_reverse_uv_sampler.hh
@@ -35,7 +35,7 @@ class ReverseUVSampler {
 
   struct Result {
     ResultType type = ResultType::None;
-    const MLoopTri *looptri = nullptr;
+    int looptri_index = -1;
     float3 bary_weights;
   };
 
diff --git a/source/blender/geometry/intern/add_curves_on_mesh.cc b/source/blender/geometry/intern/add_curves_on_mesh.cc
index a03c9b994a9..1cc13a40fc4 100644
--- a/source/blender/geometry/intern/add_curves_on_mesh.cc
+++ b/source/blender/geometry/intern/add_curves_on_mesh.cc
@@ -140,6 +140,7 @@ static void interpolate_position_with_interpolation(CurvesGeometry &curves,
                                                     const Span<float> new_lengths_cu,
                                                     const Span<float3> new_normals_su,
                                                     const bke::CurvesSurfaceTransforms &transforms,
+                                                    const Span<MLoopTri> looptris,
                                                     const ReverseUVSampler &reverse_uv_sampler,
                                                     const Span<float3> corner_normals_su)
 {
@@ -178,7 +179,7 @@ static void interpolate_position_with_interpolation(CurvesGeometry &curves,
         }
 
         const float3 neighbor_normal_su = compute_surface_point_normal(
-            *result.looptri, result.bary_weights, corner_normals_su);
+            looptris[result.looptri_index], result.bary_weights, corner_normals_su);
         const float3 neighbor_normal_cu = math::normalize(transforms.surface_to_curves_normal *
                                                           neighbor_normal_su);
 
@@ -254,7 +255,7 @@ AddCurvesOnMeshOutputs add_curves_on_mesh(CurvesGeometry &curves,
       outputs.uv_error = true;
       continue;
     }
-    const MLoopTri &looptri = *result.looptri;
+    const MLoopTri &looptri = inputs.surface_looptris[result.looptri_index];
     bary_coords.append(result.bary_weights);
     looptris.append(&looptri);
     const float3 root_position_su = attribute_math::mix3<float3>(
@@ -358,6 +359,7 @@ AddCurvesOnMeshOutputs add_curves_on_mesh(CurvesGeometry &curves,
                                             new_lengths_cu,
                                             new_normals_su,
                                             *inputs.transforms,
+                                            inputs.surface_looptris,
                                             *inputs.reverse_uv_sampler,
                                             inputs.corner_normals_su);
   }
diff --git a/source/blender/geometry/intern/reverse_uv_sampler.cc b/source/blender/geometry/intern/reverse_uv_sampler.cc
index f66e4a3ac2e..d69935a8fe3 100644
--- a/source/blender/geometry/intern/reverse_uv_sampler.cc
+++ b/source/blender/geometry/intern/reverse_uv_sampler.cc
@@ -48,7 +48,7 @@ ReverseUVSampler::Result ReverseUVSampler::sample(const float2 &query_uv) const
 
   float best_dist = FLT_MAX;
   float3 best_bary_weights;
-  const MLoopTri *best_looptri;
+  int best_looptri;
 
   /* The distance to an edge that is allowed to be inside or outside the triangle. Without this,
    * the lookup can fail for floating point accuracy reasons when the uv is almost exact on an
@@ -84,7 +84,7 @@ ReverseUVSampler::Result ReverseUVSampler::sample(const float2 &query_uv) const
     if (dist < best_dist) {
       best_dist = dist;
       best_bary_weights = bary_weights;
-      best_looptri = &looptri;
+      best_looptri = looptri_index;
     }
   }
 
diff --git a/source/blender/nodes/geometry/nodes/node_geo_deform_curves_on_surface.cc b/source/blender/nodes/geometry/nodes/node_geo_deform_curves_on_surface.cc
index dabd2a1a9f2..d926e7e4e30 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_deform_curves_on_surface.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_deform_curves_on_surface.cc
@@ -68,9 +68,11 @@ static void deform_curves(const CurvesGeometry &curves,
 
   const Span<MVert> surface_verts_old = surface_mesh_old.verts();
   const Span<MLoop> surface_loops_old = surface_mesh_old.loops();
+  const Span<MLoopTri> surface_looptris_old = surface_mesh_old.looptris();
 
   const Span<MVert> surface_verts_new = surface_mesh_new.verts();
   const Span<MLoop> surface_loops_new = surface_mesh_new.loops();
+  const Span<MLoopTri> surface_looptris_new = surface_mesh_new.looptris();
 
   threading::parallel_for(curves.curves_range(), 256, [&](const IndexRange range) {
     for (const int curve_i : range) {
@@ -85,8 +87,8 @@ static void deform_curves(const CurvesGeometry &curves,
         continue;
       }
 
-      const MLoopTri &looptri_old = *surface_sample_old.looptri;
-      const MLoopTri &looptri_new = *surface_sample_new.looptri;
+      const MLoopTri &looptri_old = surface_looptris_old[surface_sample_old.looptri_index];
+      const MLoopTri &looptri_new = surface_looptris_new[surface_sample_new.looptri_index];
       const float3 &bary_weights_old = surface_sample_old.bary_weights;
       const float3 &bary_weights_new = surface_sample_new.bary_weights;
 
diff --git a/source/blender/nodes/geometry/nodes/node_geo_sample_uv_surface.cc b/source/blender/nodes/geometry/nodes/node_geo_sample_uv_surface.cc
index 0a6dd569395..0e6e68036d9 100644
--- a/source/blender/nodes/geometry/nodes/node_geo_sample_uv_surface.cc
+++ b/source/blender/nodes/geometry/nodes/node_geo_sample_uv_surface.cc
@@ -105,9 +105,8 @@ static void node_gather_link_searches(GatherLinkSearchOpParams &params)
   }
 }
 
-class SampleUVSurfaceFunction : public fn::MultiFunction {
+class SampleMeshBarycentricFunction : public fn::MultiFunction {
   GeometrySet source_;
-  Field<float2> src_uv_map_field_;
   GField src_field_;
 
   /**
@@ -122,60 +121,50 @@ class SampleUVSurfaceFunction : public fn::MultiFunction {
   std::optional<bke::MeshFieldContext> source_context_;
   std::unique_ptr<FieldEvaluator> source_evaluator_;
   const GVArray *source_data_;
-  VArraySpan<float2> source_uv_map_;
 
-  std::optional<ReverseUVSampler> reverse_uv_sampler_;
+  Span<MLoopTri> looptris_;
 
  public:
-  SampleUVSurfaceFunction(GeometrySet geometry, Field<float2> src_uv_map_field, GField src_field)
-      : source_(std::move(geometry)),
-        src_uv_map_field_(std::move(src_uv_map_field)),
-        src_field_(std::move(src_field))
+  SampleMeshBarycentricFunction(GeometrySet geometry, GField src_field)
+      : source_(std::move(geometry)), src_field_(std::move(src_field))
   {
     source_.ensure_owns_direct_data();
-    signature_ = this->create_signature();
     this->set_signature(&signature_);
+    signature_ = this->create_signature();
     this->evaluate_source();
   }
 
   fn::MFSignature create_signature()
   {
-    blender::fn::MFSignatureBuilder signature{"Sample UV Surface"};
-    signature.single_input<float2>("Sample UV");
+    fn::MFSignatureBuilder signature{"Sample Barycentric Triangles"};
+    signature.single_input<int>("Triangle Index");
+    signature.single_input<float3>("Barycentric Weight");
     signature.single_output("Value", src_field_.cpp_type());
-    signature.single_output<bool>("Is Valid");
     return signature.build();
   }
 
   void call(IndexMask mask, fn::MFParams params, fn::MFContext /*context*/) const override
   {
-    const VArray<float2> &sample_uvs = params.readonly_single_input<float2>(0, "Sample UV");
-    GMutableSpan dst = params.uninitialized_single_output_if_required(1, "Value");
-    MutableSpan<bool> valid_dst = params.uninitialized_single_output_if_required<bool>(2,
-                                                                                       "Is Valid");
+    const VArraySpan<int> triangle_indices = params.readonly_single_input<int>(0,
+                                                                               "Triangle Index");
+    const VArraySpan<float3> bary_weights = params.readonly_single_input<float3>(
+        1, "Barycentric Weight");
+    GMutableSpan dst = params.uninitialized_single_output(2, "Value");
 
-    const CPPType &type = src_field_.cpp_type();
-    attribute_math::convert_to_static_type(type, [&](auto dummy) {
+    attribute_math::convert_to_static_type(src_field_.cpp_type(), [&](auto dummy) {
       using T = decltype(dummy);
       const VArray<T> src_typed = source_data_->typed<T>();
       MutableSpan<T> dst_typed = dst.typed<T>();
       for (const int i : mask) {
-        const float2 sample_uv = sample_uvs[i];
-        const ReverseUVSampler::Result result = reverse_uv_sampler_->sample(sample_uv);
-        const bool valid = result.type == ReverseUVSampler::ResultType::Ok;
-        if (!dst_typed.is_empty()) {
-          if (valid) {
-            dst_typed[i] = attribute_math::mix3(result.bary_weights,
-                                                src_typed[result.looptri->tri[0]],
-                                                src_typed[result.looptri->tri[1]],
-                                                src_typed[result.looptri->tri[2]]);
-          }
-          else {
-            dst_typed[i] = {};
-          }
+        const int triangle_index = triangle_indices[i];
+        if (triangle_indices[i] != -1) {
+          dst_typed[i] = attribute_math::mix3(bary_weights[i],
+                                              src_typed[looptris_[triangle_index].tri[0]],
+                                              src_typed[looptris_[triangle_index].tri[1]],
+                                              src_typed[looptris_[triangle_index].tri[2]]);
         }
-        if (!valid_dst.is_empty()) {
-          valid_dst[i] = valid;
+        else {
+          dst_typed[i] = {};
         }
       }
     });
@@ -185,14 +174,91 @@ class SampleUVSurfaceFunction : public fn::MultiFunction {
   void evaluate_source()
   {
     const Mesh &mesh = *source_.get_mesh_for_read();
+    looptris_ = mesh.looptris();
     source_context_.emplace(bke::MeshFieldContext{mesh, domain_});
     const int domain_size = mesh.attributes().domain_size(domain_);
     source_evaluator_ = std::make_unique<FieldEvaluator>(*source_context_, domain_size);
-    source_evaluator_->add(src_uv_map_field_);
     source_evaluator_->add(src_field_);
     source_evaluator_->evaluate();
+    source_data_ = &source_evaluator_->get_evaluated(0);
+  }
+};
+
+class ReverseUVSampleFunction : public fn::MultiFunction {
+  GeometrySet source_;
+  Field<float2> src_uv_map_field_;
+
+  std::optional<bke::MeshFieldContext> source_context_;
+  std::unique_ptr<FieldEvaluator> source_evaluator_;
+  VArraySpan<float2> source_uv_map_;
+
+  std::optional<ReverseUVSampler> reverse_uv_sampler_;
+
+ public:
+  ReverseUVSampleFunction(GeometrySet geometry, Field<float2> src_uv_map_field)
+      : source_(std::move(geometry)), src_uv_map_field_(std::move(src_uv_map_field))
+  {
+    source_.ensure_owns_direct_data();
+    static fn::MFSignature signature = create_signature();
+    this->set_signature(&signature);
+    this->evaluate_source();
+  }
+
+  static fn::MFSignature create_signature()
+  {
+    fn::MFSignatureBuilder signature{"Sample UV Surface"};
+    signature.single_input<float2>("Sample UV");
+    signature.single_output<bool>("Is Valid");
+    signature.single_output<int>("Triangle Index");
+    signature.single_output<float3>("Barycentric Weights");
+    return signature.build();
+  }
+
+  void call(IndexMask mask, fn::MFParams params, fn::MFContext /*context*/) const override
+  {
+    const VArraySpan<float2> sample_uvs = params.readonly_single_input<float2>(0, "Sample UV");
+    MutableSpan<bool> is_valid = params.uninitialized_single_output_if_required<bool>(1,
+                                                                                      "Is Valid");
+    MutableSpan<int> tri_index = params.uninitialized_single_output_if_required<int>(
+        2, "Triangle Index");
+    MutableSpan<float3> bary_weights = params.uninitialized_single_output_if_required<float3>(
+        3, "Barycentric Weights");
+
+    Array<ReverseUVSampler::Result> results(mask.min_array_size());
+    reverse_uv_sampler_->sample_many(sample_uvs, results);
+
+    if (!is_valid.is_empty()) {
+      std::transform(results.begin(),
+                     results.end(),
+                     is_valid.begin(),
+                     [](const ReverseUVSampler::Result &result) {
+                       return result.type == ReverseUVSampler::ResultType::Ok;
+                     });
+    }
+    if (!tri_index.is_empty()) {
+      std::transform(results.begin(),
+                     results.end(),
+                     tri_index.begin(),
+                     [](const ReverseUVSampler::Result &result) { return result.looptri_index; });
+    }
+
+    if (!bary_weights.is_empty()) {
+      std::transform(results.begin(),
+                     results.end(),
+                     bary_weights.begin(),
+                     [](const ReverseUVSampler::Result &result) { return result.bary_weights; });
+    }
+  }
+
+ private:
+  void evaluate_source()
+  {
+    const Mesh &mesh = *source_.get_mesh_for_read();
+    source_context_.emplace(bke::MeshFieldContext{mesh, ATTR_DOMAIN_CORNER});
+    source_evaluator_ = std::make_unique<FieldEvaluator>(*source_context_, mesh.totloop);
+    source_evaluator_->add(src_uv_map_field_);
+    source_evaluator_->evaluate();
     source_uv_map_ = source_evaluator_->get_evaluated<float2>(0);
-    source_data_ = &source_evaluator_->get_evaluated(1);
 
     reverse_uv_sampler_.emplace(source_uv_map_, mesh.looptris());
   }
@@ -260,19 +326,24 @@ static void node_geo_exec(GeoNodeExecParams params)
     return;
   }
 
-  const CPPType &float2_type = CPPType::get<float2>();
-
+  /* Do reverse sampling of the UV map first. */
   const bke::DataTypeConversions &conversions = bke::get_implicit_type_conversions();
+  const CPPType &float2_type = CPPType::get<float2>();
   Field<float2> source_uv_map = conversions.try_convert(
       params.extract_input<Field<float3>>("Source UV Map"), float2_type);
-  GField field = get_input_attribute_field(params, data_type);
   Field<float2> sample_uvs = conversions.try_convert(
       params.extract_input<Field<float3>>("Sample UV"), float2_type);
-  auto fn = std::make_shared<SampleUVSurfaceFunction>(
-      std::move(geometry), std::move(source_uv_map), std::move(field));
-  auto op = FieldOperation::Create(std::move(fn), {std::move(sample_uvs)});
-  output_attribute_field(params, GField(op, 0));
-  params.set_output("Is Valid", Field<bool>(op, 1));
+  auto uv_op = FieldOperation::Create(
+      std::make_shared<ReverseUVSampleFunction>(geometry, std::move(source_uv_map)),
+      {std::move(sample_uvs)});
+  params.set_output("Is Valid", Field<bool>(uv_op, 0));
+
+  /* Use the output of the UV sampling to interpolate the mesh attribute. */
+  GField field = get_input_attribute_field(params, data_type);
+  auto sample_op = FieldOperation::Create(
+      std::make_shared<SampleMeshBarycentricFunction>(std::move(geometry), std::move(field)),
+      {Field<int>(uv_op, 1), Field<float3>(uv_op, 2)});
+  output_attribute_field(params, GField(sample_op, 0));
 }
 
 }  // namespace blender::nodes::node_geo_sample_uv_surface_cc