Modeling: use the central vertex when finding duplicates

Previously the method of picking the "target" duplicate wasn't
deterministic from a user perspective.

The behavior has been changed so:

- For a cluster of 3 or more vertices,
  use the vertex closest to the centroid.
- For a cluster of 2 use the lowest index.

This mitigates #78916, solving some cases where clusters have a
central vertex although can't be considered fixed as the 2 vertex
case doesn't work so well.

Added a BLI_kdtree_{N}d_calc_duplicates_cb function that lets the
caller choose the index to keep from a cluster of duplicates.

Refactored from !145851.

Ref !146492

Co-authored-by: tariqsulley <tariqsulley3c@gmail.com>

source/blender/blenlib/BLI_kdtree_impl.h

@@ -68,6 +68,29 @@ int BLI_kdtree_nd_(calc_duplicates_fast)(const KDTree *tree,
                                          bool use_index_order,
                                          int *duplicates);
 
+/**
+ * De-duplicate utility where the callback can evaluate duplicates and select the target
+ * which other indices are merged into.
+ *
+ * \param tree: A tree, all indices *must* be unique.
+ * \param deduplicate_cb: A function which receives duplicate indices,
+ * it must choose the the "target" index to keep which is returned.
+ * The return value is an index in the `cluster` array (a value from `0..cluster_num`).
+ * The last item in `cluster` is the index from which the search began.
+ *
+ * \note ~1.1x-1.5x slower than `calc_duplicates_fast` depending on the distribution of points.
+ *
+ * \note The duplicate search is performed in an order defined by the tree-nodes index,
+ * the index of the input (first to last) for predictability.
+ */
+int BLI_kdtree_nd_(calc_duplicates_cb)(const KDTree *tree,
+                                       const float range,
+                                       int *duplicates,
+                                       int (*deduplicate_cb)(void *user_data,
+                                                             const int *cluster,
+                                                             int cluster_num),
+                                       void *user_data);
+
 int BLI_kdtree_nd_(deduplicate)(KDTree *tree);
 
 /** Versions of find/range search that take a squared distance callback to support bias. */
@@ -124,6 +147,23 @@ inline int BLI_kdtree_nd_(find_nearest_cb_cpp)(const KDTree *tree,
       r_nearest);
 }
 
+template<typename Fn>
+inline int BLI_kdtree_nd_(calc_duplicates_cb_cpp)(const KDTree *tree,
+                                                  float distance,
+                                                  int *duplicates,
+                                                  const Fn &fn)
+{
+  return BLI_kdtree_nd_(calc_duplicates_cb)(
+      tree,
+      distance,
+      duplicates,
+      [](void *user_data, const int *cluster, int cluster_num) -> int {
+        const Fn &fn = *static_cast<const Fn *>(user_data);
+        return fn(cluster, cluster_num);
+      },
+      const_cast<Fn *>(&fn));
+}
+
 #undef _BLI_CONCAT_AUX
 #undef _BLI_CONCAT
 #undef BLI_kdtree_nd_

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"
@@ -934,6 +936,75 @@ int BLI_kdtree_nd_(calc_duplicates_fast)(const KDTree *tree,
 
 /** \} */
 
+/* -------------------------------------------------------------------- */
+/** \name BLI_kdtree_3d_calc_duplicates_cb
+ * \{ */
+
+int BLI_kdtree_nd_(calc_duplicates_cb)(const KDTree *tree,
+                                       const float range,
+                                       int *duplicates,
+                                       int (*duplicates_cb)(void *user_data,
+                                                            const int *cluster,
+                                                            int cluster_num),
+                                       void *user_data)
+{
+  BLI_assert(tree->is_balanced);
+  if (UNLIKELY(tree->root == KD_NODE_UNSET)) {
+    return 0;
+  }
+
+  /* Use `index_to_node_index` so coordinates are looked up in order first to last. */
+  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;
+    }
+
+    BLI_assert(cluster.is_empty());
+    const float *search_co = tree->nodes[index_to_node_index[node_index]].co;
+    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;
+    }
+    found += int(cluster.size());
+    cluster.append(node_index);
+
+    const int cluster_index = duplicates_cb(user_data, cluster.data(), int(cluster.size()));
+    BLI_assert(uint(cluster_index) < uint(cluster.size()));
+    const int target_index = cluster[cluster_index];
+    for (const int cluster_node_index : cluster) {
+      duplicates[cluster_node_index] = target_index;
+    }
+    cluster.clear();
+  }
+
+  return found;
+}
+
+/** \} */
+
 /* -------------------------------------------------------------------- */
 /** \name BLI_kdtree_3d_deduplicate
  * \{ */

source/blender/bmesh/operators/bmo_removedoubles.cc

@@ -648,7 +648,48 @@ 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;
+
+  /* Given a cluster of duplicates, pick the index to keep. */
+  auto deduplicate_target_calc_fn = [&verts](const int *cluster, const int cluster_num) -> int {
+    if (cluster_num == 2) {
+      /* Special case, no use in calculating centroid.
+       * Use the lowest index for stability. */
+      return (cluster[0] < cluster[1]) ? 0 : 1;
+    }
+    BLI_assert(cluster_num > 2);
+
+    blender::float3 centroid{0.0f};
+    for (int i = 0; i < cluster_num; i++) {
+      centroid += blender::float3(verts[cluster[i]]->co);
+    }
+    centroid /= float(cluster_num);
+
+    /* Now pick the most "central" index (with lowest index as a tie breaker). */
+    const int cluster_end = cluster_num - 1;
+    /* Assign `i_best` from the last index as this is the index where the search originated
+     * so it's most likely to be the best. */
+    int i_best = cluster_end;
+    float dist_sq_best = len_squared_v3v3(centroid, verts[cluster[i_best]]->co);
+    for (int i = 0; i < cluster_end; i++) {
+      const float dist_sq_test = len_squared_v3v3(centroid, verts[cluster[i]]->co);
+
+      if (dist_sq_test > dist_sq_best) {
+        continue;
+      }
+      if (dist_sq_test == dist_sq_best) {
+        if (cluster[i] > cluster[i_best]) {
+          continue;
+        }
+      }
+      i_best = i;
+      dist_sq_best = dist_sq_test;
+    }
+    return i_best;
+  };
+
+  found_duplicates = BLI_kdtree_3d_calc_duplicates_cb_cpp(
+                         tree, dist, duplicates, deduplicate_target_calc_fn) != 0;
+
   BLI_kdtree_3d_free(tree);
 
   if (!found_duplicates) {