diff --git a/source/blender/blenlib/BLI_kdtree_impl.h b/source/blender/blenlib/BLI_kdtree_impl.h
index 775e0d0f87e..d83b0e0dba3 100644
--- a/source/blender/blenlib/BLI_kdtree_impl.h
+++ b/source/blender/blenlib/BLI_kdtree_impl.h
@@ -68,6 +68,18 @@ int BLI_kdtree_nd_(calc_duplicates_fast)(const KDTree *tree,
                                          bool use_index_order,
                                          int *duplicates);
 
+/**
+ * Stable clustering and centroid computation to ensure consistent survivor selection.
+ *
+ * \param tree: A tree, all indices *must* be unique.
+ *
+ * \note ~1.1x-1.5x slower than `calc_duplicates_fast` depending on the distribution of points.
+ */
+int BLI_kdtree_nd_(calc_duplicates_and_center)(const KDTree *tree,
+                                               const float range,
+                                               int *duplicates,
+                                               float (*r_cluster_center)[KD_DIMS]);
+
 int BLI_kdtree_nd_(deduplicate)(KDTree *tree);
 
 /** Versions of find/range search that take a squared distance callback to support bias. */
diff --git a/source/blender/blenlib/intern/kdtree_impl.h b/source/blender/blenlib/intern/kdtree_impl.h
index 6e25a808f9c..d3c9b5313dd 100644
--- a/source/blender/blenlib/intern/kdtree_impl.h
+++ b/source/blender/blenlib/intern/kdtree_impl.h
@@ -8,6 +8,8 @@
 
 #include "MEM_guardedalloc.h"
 
+#include "BLI_array.hh"
+#include "BLI_bit_vector.hh"
 #include "BLI_kdtree_impl.h"
 #include "BLI_math_base.h"
 #include "BLI_utildefines.h"
@@ -66,6 +68,20 @@ struct KDTree {
 /** \name Local Math API
  * \{ */
 
+static void add_vn_vn(float v0[KD_DIMS], const float v1[KD_DIMS])
+{
+  for (uint j = 0; j < KD_DIMS; j++) {
+    v0[j] += v1[j];
+  }
+}
+
+static void div_vn_fl(float v[KD_DIMS], const float f)
+{
+  for (uint j = 0; j < KD_DIMS; j++) {
+    v[j] /= f;
+  }
+}
+
 static void copy_vn_vn(float v0[KD_DIMS], const float v1[KD_DIMS])
 {
   for (uint j = 0; j < KD_DIMS; j++) {
@@ -934,6 +950,83 @@ int BLI_kdtree_nd_(calc_duplicates_fast)(const KDTree *tree,
 
 /** \} */
 
+/* -------------------------------------------------------------------- */
+/** \name BLI_kdtree_3d_calc_duplicates_and_center
+ * \{ */
+
+int BLI_kdtree_nd_(calc_duplicates_and_center)(const KDTree *tree,
+                                               const float range,
+                                               int *duplicates,
+                                               float (*r_cluster_center)[KD_DIMS])
+{
+  BLI_assert(tree->is_balanced);
+  if (UNLIKELY(tree->root == KD_NODE_UNSET)) {
+    return 0;
+  }
+
+  const uint nodes_len = tree->nodes_len;
+  blender::Array<int> index_to_node_index(tree->max_node_index + 1);
+  for (uint i = 0; i < nodes_len; i++) {
+    index_to_node_index[tree->nodes[i].index] = int(i);
+  }
+
+  blender::BitVector<> visited(tree->max_node_index + 1, false);
+  /* Could be inline, declare here to avoid re-allocation. */
+  blender::Vector<int> cluster;
+
+  int found = 0;
+
+  for (uint i = 0; i < nodes_len; i++) {
+    const int node_index = tree->nodes[i].index;
+    if ((duplicates[node_index] != -1) || visited[node_index]) {
+      continue;
+    }
+
+    const float *search_co = tree->nodes[index_to_node_index[node_index]].co;
+    BLI_assert(search_co != nullptr);
+    BLI_assert(cluster.is_empty());
+
+    visited[node_index].set();
+    auto accumulate_neighbors_fn = [&duplicates, &visited, &cluster](int neighbor_index,
+                                                                     const float * /*co*/,
+                                                                     float /*dist_sq*/) -> bool {
+      if ((duplicates[neighbor_index] == -1) && !visited[neighbor_index]) {
+        cluster.append(neighbor_index);
+        visited[neighbor_index].set();
+      }
+      return true;
+    };
+
+    BLI_kdtree_nd_(range_search_cb_cpp)(tree, search_co, range, accumulate_neighbors_fn);
+    if (cluster.is_empty()) {
+      continue;
+    }
+    cluster.append(node_index);
+
+    /* Compute centroid and choose survivor in one pass. */
+    float centroid[KD_DIMS] = {};
+    int survivor_index = node_index; /* Same as `cluster[0]`. */
+
+    for (int cluster_node_index : cluster) {
+      add_vn_vn(centroid, tree->nodes[index_to_node_index[cluster_node_index]].co);
+      survivor_index = std::min(cluster_node_index, survivor_index);
+    }
+
+    div_vn_fl(centroid, float(cluster.size()));
+    copy_vn_vn(r_cluster_center[survivor_index], centroid);
+
+    for (int cluster_node_index : cluster) {
+      duplicates[cluster_node_index] = survivor_index;
+    }
+    found += int(cluster.size()) - 1;
+    cluster.clear();
+  }
+
+  return found;
+}
+
+/** \} */
+
 /* -------------------------------------------------------------------- */
 /** \name BLI_kdtree_3d_deduplicate
  * \{ */
diff --git a/source/blender/bmesh/operators/bmo_removedoubles.cc b/source/blender/bmesh/operators/bmo_removedoubles.cc
index e75ae876fbf..620c3637387 100644
--- a/source/blender/bmesh/operators/bmo_removedoubles.cc
+++ b/source/blender/bmesh/operators/bmo_removedoubles.cc
@@ -633,8 +633,13 @@ static int *bmesh_find_doubles_by_distance_impl(BMesh *bm,
                                                 const float dist,
                                                 const bool has_keep_vert)
 {
+  if (verts_len == 0) {
+    return nullptr;
+  }
+
   int *duplicates = MEM_malloc_arrayN<int>(verts_len, __func__);
   bool found_duplicates = false;
+  blender::Vector<blender::float3> survivor_cos(verts_len);
 
   KDTree_3d *tree = BLI_kdtree_3d_new(verts_len);
   for (int i = 0; i < verts_len; i++) {
@@ -648,13 +653,23 @@ static int *bmesh_find_doubles_by_distance_impl(BMesh *bm,
   }
 
   BLI_kdtree_3d_balance(tree);
-  found_duplicates = BLI_kdtree_3d_calc_duplicates_fast(tree, dist, false, duplicates) != 0;
+  found_duplicates = BLI_kdtree_3d_calc_duplicates_and_center(
+                         tree, dist, duplicates, (float(*)[3])survivor_cos.data()) != 0;
+
   BLI_kdtree_3d_free(tree);
 
   if (!found_duplicates) {
     MEM_freeN(duplicates);
     duplicates = nullptr;
   }
+  else {
+    for (int i = 0; i < verts_len; i++) {
+      if (duplicates[i] == i) {
+        copy_v3_v3(verts[i]->co, survivor_cos[i]);
+      }
+    }
+  }
+
   return duplicates;
 }